From a4fd9e6f6b09f08f50fae3df6fb04f8e298b0503 Mon Sep 17 00:00:00 2001
From: Lijun Zhu <ljzhu@gps.caltech.edu>
Date: Thu, 25 Jul 2019 15:44:53 -0700
Subject: [PATCH] add positive definteness check for covariance matrix of the
 Gaussian proposal

---
 altar/bayesian/COV.py |  12 +++-
 ext/altar.cc          |   4 ++
 ext/condition.cc      | 127 ++++++++++++++++++++++++++++++++++++++++++
 ext/condition.h       |  27 +++++++++
 4 files changed, 167 insertions(+), 3 deletions(-)
 create mode 100644 ext/condition.cc
 create mode 100644 ext/condition.h

diff --git a/altar/bayesian/COV.py b/altar/bayesian/COV.py
index 0370dea1..c817a5b0 100644
--- a/altar/bayesian/COV.py
+++ b/altar/bayesian/COV.py
@@ -43,6 +43,11 @@ class COV(altar.component, family="altar.schedulers.cov", implements=scheduler):
     maxiter = altar.properties.int(default=10**3)
     maxiter.doc = 'the maximum number of iterations while looking for a δβ'
 
+    check_positive_definiteness = altar.properties.bool(default=True)
+    check_positive_definiteness.doc = 'whether to check the positive definiteness of Σ matrix and condition it accordingly'
+
+    min_eigenvalue_ratio = altar.properties.float(default=0.001)
+    min_eigenvalue_ratio.doc = 'the desired minimal eigenvalue of Σ matrix, as a ratio to the max eigenvalue'
 
     # public data
     w = None # the vector of re-sampling weights
@@ -158,7 +163,8 @@ def computeCovariance(self, step):
                 Σ[j,i] = Σ[i,j]
 
         # condition the covariance matrix
-        self.conditionCovariance(Σ=Σ)
+        if self.check_positive_definiteness:
+            self.conditionCovariance(Σ=Σ)
 
         # all done
         return Σ
@@ -221,8 +227,8 @@ def conditionCovariance(self, Σ):
         """
         Make sure the covariance matrix Σ is symmetric and positive definite
         """
-        # no need to symmetrize it since it is symmetric by construction
-        # NYI: check the eigenvalues to verify positive definiteness
+        # replaces negative or small eigenvalues with min_eigenvalue_ratio*max_eigenvalue
+        altar.libaltar.matrix_condition(Σ.data, self.min_eigenvalue_ratio)
         # all done
         return Σ
 
diff --git a/ext/altar.cc b/ext/altar.cc
index 7a62bddf..c29683a9 100644
--- a/ext/altar.cc
+++ b/ext/altar.cc
@@ -17,6 +17,7 @@
 #include "exceptions.h"
 #include "metadata.h"
 #include "dbeta.h"
+#include "condition.h"
 
 
 // put everything in my private namespace
@@ -36,6 +37,9 @@ namespace altar {
             { cov__name__, cov, METH_VARARGS, cov__doc__},
             { dbeta__name__, dbeta, METH_VARARGS, dbeta__doc__},
 
+            // matrix_condition for positive definite
+            { matrix_condition__name__, matrix_condition, METH_VARARGS, matrix_condition__doc__},
+
             // sentinel
             {0, 0, 0, 0}
         };
diff --git a/ext/condition.cc b/ext/condition.cc
new file mode 100644
index 00000000..27911654
--- /dev/null
+++ b/ext/condition.cc
@@ -0,0 +1,127 @@
+// -*- C++ -*-
+// -*- coding: utf-8 -*-
+//
+// (c) 2013-2019 parasim inc
+// (c) 2010-2019 california institute of technology
+// all rights reserved
+//
+// Author(s): AlTar-1 team, rearranged by Lijun Zhu
+
+#include <portinfo>
+#include <Python.h>
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+
+#include <gsl/gsl_sys.h>
+#include <gsl/gsl_min.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_roots.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_statistics.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_sort_vector.h>
+#include <gsl/gsl_eigen.h>
+
+#include <pyre/journal.h>
+#include <pyre/gsl/capsules.h>
+
+// local includes
+#include "condition.h"
+#include "capsules.h"
+
+
+// condition a matrix to be positive definite
+// replace negative or small eigen-values with ratio*max_eigenvalue
+const char * const altar::extensions::matrix_condition__name__ = "matrix_condition";
+const char * const altar::extensions::matrix_condition__doc__ =
+    "condition a matrix to be positive definite";
+
+PyObject *
+altar::extensions::matrix_condition(PyObject *, PyObject * args)
+{
+    // the arguments
+    PyObject * matrixCapsule;
+
+     // set minimum eigenvalue ratio
+    double eval_ratio_min;
+
+    // build my debugging channel
+    pyre::journal::debug_t debug("altar.matrix_condition");
+
+    // unpack the argument tuple
+    int status = PyArg_ParseTuple(
+                                  args, "O!d:matrix_condition",
+                                  &PyCapsule_Type, &matrixCapsule,
+                                  &eval_ratio_min
+                                  );
+    // if something went wrong
+    if (!status) return 0;
+    // bail out if the {matrix} capsule is not valid
+    if (!PyCapsule_IsValid(matrixCapsule, gsl::matrix::capsule_t)) {
+        PyErr_SetString(PyExc_TypeError, "invalid matrix capsule");
+        return 0;
+    }
+
+    // get the {sigma} matrix
+    gsl_matrix * sigma =
+        static_cast<gsl_matrix *>(PyCapsule_GetPointer(matrixCapsule, gsl::matrix::capsule_t));
+
+    // get matrix size
+    size_t m = sigma->size1;
+
+    // solve the eigen value problem
+    gsl_vector *eval = gsl_vector_alloc (m);
+    gsl_matrix *evec = gsl_matrix_alloc (m, m);
+    gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (m);
+
+    gsl_eigen_symmv (sigma, eval, evec, w);
+    gsl_eigen_symmv_free (w);
+
+    // sort the eigen values in ascending order (magnitude)
+    gsl_eigen_symmv_sort (eval, evec,  GSL_EIGEN_SORT_ABS_ASC);
+
+    // make a transpose of the eigen vector matrix
+    gsl_matrix *evecT = gsl_matrix_alloc(m,m);
+    gsl_matrix_transpose_memcpy(evecT, evec);
+
+    // allocate a matrix for conditioned eigen values
+    gsl_matrix *diagM = gsl_matrix_calloc(m,m);
+
+    // set the minimum eigen value as the max * ratio
+    double eval_min = eval_ratio_min*gsl_vector_get(eval,m-1);
+    double eval_i;
+    // copy the eigenvalues, set it to eval_min if smaller
+    for (size_t i = 0; i < m; i++)
+    {
+        eval_i  = gsl_vector_get (eval, i);
+        if (eval_i<eval_min) gsl_matrix_set(diagM, i, i, eval_min);
+        else gsl_matrix_set(diagM, i, i, eval_i);
+    }
+
+    // reconstruct sigma from the conditioned eigen values
+    gsl_matrix *tmp = gsl_matrix_alloc(m,m);
+    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, diagM, evecT, 0.0, tmp);
+    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, evec, tmp, 0.0, sigma);
+
+    // make sigma symmetric
+    gsl_matrix_transpose_memcpy(tmp, sigma);
+    gsl_matrix_add(sigma,tmp);
+    gsl_matrix_scale(sigma, 0.5);
+
+    // free temporary data
+    gsl_vector_free (eval);
+    gsl_matrix_free (evec);
+    gsl_matrix_free (evecT);
+    gsl_matrix_free (diagM);
+    gsl_matrix_free (tmp);
+
+    // all done
+    // return None
+    Py_INCREF(Py_None);
+    return Py_None;
+}
+
+// end of file
diff --git a/ext/condition.h b/ext/condition.h
new file mode 100644
index 00000000..71eb5c40
--- /dev/null
+++ b/ext/condition.h
@@ -0,0 +1,27 @@
+// -*- C++ -*-
+// -*- coding: utf-8 -*-
+//
+// (c) 2013-2019 parasim inc
+// (c) 2010-2019 california institute of technology
+// all rights reserved
+//
+// Author(s): Lijun Zhu
+
+#if !defined(altar_extensions_condition_h)
+#define altar_extensions_condition_h
+
+
+// place everything in my private namespace
+namespace altar {
+    namespace extensions {
+
+            extern const char * const matrix_condition__name__;
+            extern const char * const matrix_condition__doc__;
+            PyObject * matrix_condition(PyObject *, PyObject *);
+
+    } // of namespace extensions
+} // of namespace altar
+
+#endif //altar_extensions_condition_h
+
+// end of file