Merge pull request #110 from QInfer/fix-lw-cov-esti

Fix for #108.

Author: cgranade

src/qinfer/resamplers.py

@@ -42,7 +42,7 @@
 import scipy.linalg as la
 import warnings
 
-from .utils import outer_product, particle_meanfn, particle_covariance_mtx
+from .utils import outer_product, particle_meanfn, particle_covariance_mtx, sqrtm_psd
 
 from abc import ABCMeta, abstractmethod, abstractproperty
 from future.utils import with_metaclass
@@ -260,7 +260,7 @@ def __call__(self, model, particle_weights, particle_locations,
                 ResamplerWarning
             )
             cov = self._zero_cov_comp * np.eye(cov.shape[0])
-        S, S_err = la.sqrtm(cov, disp=False)
+        S, S_err = sqrtm_psd(cov)
         if not np.isfinite(S_err):
             raise ResamplerError(
                 "Infinite error in computing the square root of the "

src/qinfer/tests/test_utils.py

@@ -36,9 +36,11 @@
 import numpy as np
 from numpy.testing import assert_equal, assert_almost_equal
 
+from scipy.linalg import sqrtm
+
 from qinfer.tests.base_test import DerandomizedTestCase, MockModel, assert_warns
 
-from qinfer.utils import in_ellipsoid, assert_sigfigs_equal
+from qinfer.utils import in_ellipsoid, assert_sigfigs_equal, sqrtm_psd
 
 ## TESTS #####################################################################
 
@@ -80,7 +82,7 @@ def test_assert_sigfigs_equal(self):
             np.array([1728]),
             4
         )
-        
+
 
 class TestEllipsoids(DerandomizedTestCase):
 
@@ -116,4 +118,26 @@ def test_in_ellipsoid(self):
         assert_equal(
             in_ellipsoid(x, A, c), 
             np.array([1,0,1,0], dtype=bool)
-        )
+        )
+
+class TestLinearAlgebra(DerandomizedTestCase):
+    def test_sqrtm_psd(self):
+        # Construct Y = XX^T as a PSD matrix.
+        X = np.random.random((5, 5))
+        Y = np.dot(X, X.T)
+        sqrt_Y = sqrtm_psd(Y, est_error=False)
+
+        np.testing.assert_allclose(
+            np.dot(sqrt_Y, sqrt_Y),
+            Y
+        )
+
+        # Try again, but with a singular matrix.
+        Y_singular = np.zeros((6, 6))
+        Y_singular[:5, :5] = Y
+        sqrt_Y_singular = sqrtm_psd(Y_singular, est_error=False)
+
+        np.testing.assert_allclose(
+            np.dot(sqrt_Y_singular, sqrt_Y_singular),
+            Y_singular
+        )

src/qinfer/utils.py

@@ -38,6 +38,8 @@
 import numpy as np
 import numpy.linalg as la
 
+from scipy.linalg import eigh
+
 from scipy.stats import logistic, binom
 from scipy.special import gammaln, gamma
 from scipy.linalg import sqrtm
@@ -455,6 +457,21 @@ def safe_shape(arr, idx=0, default=1):
     shape = np.shape(arr)
     return shape[idx] if idx < len(shape) else default
 
+def sqrtm_psd(A, est_error=True, check_finite=True):
+    """
+    Returns the matrix square root of a positive semidefinite matrix,
+    truncating negative eigenvalues.
+    """
+    w, v = eigh(A, check_finite=check_finite)
+    mask = w <= 0
+    w[mask] = 0
+    np.sqrt(w, out=w)
+    A_sqrt = (v * w).dot(v.conj().T)
+
+    if est_error:
+        return A_sqrt, np.linalg.norm(np.dot(A_sqrt, A_sqrt) - A, 'fro')
+    else:
+        return A_sqrt
 
 #==============================================================================
 #Test Code