Permalink
Browse files

BUG: pinvh: account for non-positive hermitian matrices

  • Loading branch information...
1 parent 8ffd440 commit 4c45c08d16434a0eb44efcb45e069dc81c6bd482 @jakevdp jakevdp committed with pv Aug 16, 2012
Showing with 15 additions and 4 deletions.
  1. +2 −1 scipy/linalg/basic.py
  2. +13 −3 scipy/linalg/tests/test_basic.py
View
3 scipy/linalg/basic.py
@@ -648,7 +648,8 @@ def pinvh(a, cond=None, rcond=None):
factor = {'f': 1E3, 'd': 1E6}
cond = factor[t] * np.finfo(t).eps
- above_cutoff = (s > cond * np.max(s))
+ # unlike svd case, eigh can lead to negative eigenvalues
+ above_cutoff = (abs(s) > cond * np.max(abs(s)))
psigma_diag = np.zeros_like(s)
psigma_diag[above_cutoff] = 1.0 / s[above_cutoff]
View
16 scipy/linalg/tests/test_basic.py
@@ -541,8 +541,8 @@ def test_simple_real(self):
assert_array_almost_equal(dot(a,a_pinv), np.eye(3))
def test_simple_complex(self):
- a = (array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
- + 1j * array([[9, 8, 7], [6, 5, 4], [3, 2, 1]]))
+ a = (array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
+ + 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]]))
a_pinv = pinv(a)
assert_array_almost_equal(dot(a, a_pinv), np.eye(3))
a_pinv = pinv2(a)
@@ -567,14 +567,24 @@ def test_simple_rows(self):
assert_array_almost_equal(a_pinv,a_pinv2)
-class TestPinv(TestCase):
+class TestPinvSymmetric(TestCase):
def test_simple_real(self):
a = array([[1,2,3], [4,5,6.], [7,8,10]])
a = np.dot(a, a.T)
a_pinv = pinvh(a)
assert_array_almost_equal(np.dot(a, a_pinv), np.eye(3))
+ def test_nonpositive(self):
+ a = array([[1,2,3], [4,5,6.], [7,8,9]])
+ a = np.dot(a, a.T)
+ u, s, vt = np.linalg.svd(a)
+ s[0] *= -1
+ a = np.dot(u * s, vt) # a is now symmetric non-positive and singular
+ a_pinv = pinv2(a)
+ a_pinvh = pinvh(a)
+ assert_array_almost_equal(a_pinv, a_pinvh)
+
def test_simple_complex(self):
a = (array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
+ 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]]))

0 comments on commit 4c45c08

Please sign in to comment.