From 92ff36b916b4ad96e6ff9cb5f2f6672d3b42991b Mon Sep 17 00:00:00 2001
From: Kevin Murray <spam@kdmurray.id.au>
Date: Sun, 10 Jan 2016 17:16:10 +1100
Subject: [PATCH] Warn users on non-PSD kernel matrix

---
 src/kernel.cc     | 8 ++++++++
 src/kwip-utils.cc | 9 +++++++++
 src/kwip-utils.hh | 5 +++++
 3 files changed, 22 insertions(+)

diff --git a/src/kernel.cc b/src/kernel.cc
index 0994a46..74fe953 100644
--- a/src/kernel.cc
+++ b/src/kernel.cc
@@ -113,6 +113,14 @@ calculate_pairwise(std::vector<std::string> &hash_fnames)
     if (verbosity > 0) {
         *outstream << "Done all!" << std::endl;
     }
+
+    if (!matrix_is_pos_semidef(_kernel_m)) {
+        *outstream << "WARNING: The kernel matrix is not positive semidefinite."
+                   << std::endl;
+    } else {
+        *outstream << "The kernel matrix is positive semidefinite."
+                   << std::endl;
+    }
 }
 
 void
diff --git a/src/kwip-utils.cc b/src/kwip-utils.cc
index b44cdb1..1d921b3 100644
--- a/src/kwip-utils.cc
+++ b/src/kwip-utils.cc
@@ -11,6 +11,8 @@
 
 
 #include "kwip-utils.hh"
+#include <Eigen/Eigenvalues>
+
 
 namespace kwip
 {
@@ -120,4 +122,11 @@ kernel_to_distance(MatrixXd &dist, MatrixXd &kernel, bool normalise)
     }
 }
 
+bool
+matrix_is_pos_semidef(MatrixXd &mat)
+{
+    VectorXd eigenvalues = mat.eigenvalues().real();
+    return eigenvalues.minCoeff() > -1e-5;
+}
+
 } // end namespace kwip
diff --git a/src/kwip-utils.hh b/src/kwip-utils.hh
index c3ca2c1..109439f 100644
--- a/src/kwip-utils.hh
+++ b/src/kwip-utils.hh
@@ -22,6 +22,7 @@
 #include <Eigen/Core>
 
 using Eigen::MatrixXd;
+using Eigen::VectorXd;
 
 #include <kwip-config.hh>
 
@@ -37,6 +38,10 @@ void print_lsmat(MatrixXd &mat, std::ostream &outstream,
 void normalise_matrix(MatrixXd &norm, MatrixXd &input);
 void kernel_to_distance(MatrixXd &dist, MatrixXd &kernel, bool normalise=true);
 
+// Checks if a matrix is postitive semi-definite. Specifically, that all
+// eigenvalues are > -1e-5.
+bool matrix_is_pos_semidef(MatrixXd &mat);
+
 template <typename eltype>
 eltype
 vec_min(std::vector<eltype> &vec)