Browse files

ENH: linalg: Merge PR #289 Pinv enhancements

This is a PR based on some of @vene's work in scikit-learn: see
scikit-learn/scikit-learn#1015

There are three enhancements:

- scipy.linalg.pinv2 is sped up by only computing necessary singular
  values, and not allocating the entire psigma matrix

- scipy.linalg.pinvh added: this uses eigh to significantly speed up the
  pseudo-inverse computation in the symmetric/hermitian case

- add return_rank keyword to obtain the effective rank from all pinv*
  routines.
  • Loading branch information...
2 parents 05a1fbe + 6cc12fe commit 5adeb123e4e6ec064dd4877870fed9b6f7d032be @pv pv committed Oct 21, 2012
Showing with 152 additions and 28 deletions.
  1. +1 −0 scipy/linalg/__init__.py
  2. +104 −19 scipy/linalg/basic.py
  3. +47 −9 scipy/linalg/tests/test_basic.py
View
1 scipy/linalg/__init__.py
@@ -31,6 +31,7 @@
lstsq - Solve a linear least-squares problem
pinv - Pseudo-inverse (Moore-Penrose) using lstsq
pinv2 - Pseudo-inverse using svd
+ pinvh - Pseudo-inverse of hermitian matrix
kron - Kronecker product of two arrays
tril - Construct a lower-triangular matrix from a given matrix
triu - Construct an upper-triangular matrix from a given matrix
View
123 scipy/linalg/basic.py
@@ -2,9 +2,10 @@
# Author: Pearu Peterson, March 2002
#
# w/ additions by Travis Oliphant, March 2002
+# and Jake Vanderplas, August 2012
__all__ = ['solve', 'solve_triangular', 'solveh_banded', 'solve_banded',
- 'inv', 'det', 'lstsq', 'pinv', 'pinv2']
+ 'inv', 'det', 'lstsq', 'pinv', 'pinv2', 'pinvh']
import numpy as np
@@ -13,7 +14,7 @@
from misc import LinAlgError, _datacopied
from scipy.linalg import calc_lwork
from decomp_schur import schur
-import decomp_svd
+import decomp, decomp_svd
# Linear equations
@@ -502,7 +503,7 @@ def lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False):
return x, resids, rank, s
-def pinv(a, cond=None, rcond=None):
+def pinv(a, cond=None, rcond=None, return_rank=False):
"""
Compute the (Moore-Penrose) pseudo-inverse of a matrix.
@@ -517,11 +518,15 @@ def pinv(a, cond=None, rcond=None):
Cutoff for 'small' singular values in the least-squares solver.
Singular values smaller than ``rcond * largest_singular_value``
are considered zero.
+ return_rank : bool, optional
+ if True, return the effective rank of the matrix
Returns
-------
B : array, shape (N, M)
The pseudo-inverse of matrix `a`.
+ rank : int
+ The effective rank of the matrix. Returned if return_rank == True
Raises
------
@@ -542,10 +547,16 @@ def pinv(a, cond=None, rcond=None):
b = np.identity(a.shape[0], dtype=a.dtype)
if rcond is not None:
cond = rcond
- return lstsq(a, b, cond=cond)[0]
+
+ x, resids, rank, s = lstsq(a, b, cond=cond)
+
+ if return_rank:
+ return x, rank
+ else:
+ return x
-def pinv2(a, cond=None, rcond=None):
+def pinv2(a, cond=None, rcond=None, return_rank=False):
"""
Compute the (Moore-Penrose) pseudo-inverse of a matrix.
@@ -562,11 +573,15 @@ def pinv2(a, cond=None, rcond=None):
Singular values smaller than ``rcond*largest_singular_value``
are considered zero.
If None or -1, suitable machine precision is used.
+ return_rank : bool, optional
+ if True, return the effective rank of the matrix
Returns
-------
B : array, shape (N, M)
The pseudo-inverse of matrix `a`.
+ rank : int
+ The effective rank of the matrix. Returned if return_rank == True
Raises
------
@@ -584,21 +599,91 @@ def pinv2(a, cond=None, rcond=None):
"""
a = np.asarray_chkfinite(a)
- u, s, vh = decomp_svd.svd(a)
- t = u.dtype.char
+ u, s, vh = decomp_svd.svd(a, full_matrices=False)
+
if rcond is not None:
cond = rcond
if cond in [None,-1]:
- eps = np.finfo(np.float).eps
- feps = np.finfo(np.single).eps
- _array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
- cond = {0: feps*1e3, 1: eps*1e6}[_array_precision[t]]
- m, n = a.shape
- cutoff = cond*np.maximum.reduce(s)
- psigma = np.zeros((m, n), t)
- for i in range(len(s)):
- if s[i] > cutoff:
- psigma[i,i] = 1.0/np.conjugate(s[i])
- #XXX: use lapack/blas routines for dot
- return np.transpose(np.conjugate(np.dot(np.dot(u,psigma),vh)))
+ t = u.dtype.char.lower()
+ factor = {'f': 1E3, 'd': 1E6}
+ cond = factor[t] * np.finfo(t).eps
+
+ above_cutoff = (s > cond * np.max(s))
+ psigma_diag = np.zeros_like(s)
+ psigma_diag[above_cutoff] = 1.0 / s[above_cutoff]
+
+ B = np.transpose(np.conjugate(np.dot(u * psigma_diag, vh)))
+
+ if return_rank:
+ return B, np.sum(above_cutoff)
+ else:
+ return B
+
+
+def pinvh(a, cond=None, rcond=None, lower=True, return_rank=False):
+ """Compute the (Moore-Penrose) pseudo-inverse of a hermetian matrix.
+
+ Calculate a generalized inverse of a symmetric matrix using its
+ eigenvalue decomposition and including all 'large' eigenvalues.
+
+ Parameters
+ ----------
+ a : array, shape (N, N)
+ Real symmetric or complex hermetian matrix to be pseudo-inverted
+ cond, rcond : float or None
+ Cutoff for 'small' eigenvalues.
+ Singular values smaller than rcond * largest_eigenvalue are considered
+ zero.
+
+ If None or -1, suitable machine precision is used.
+ lower : boolean
+ Whether the pertinent array data is taken from the lower or upper
+ triangle of a. (Default: lower)
+ return_rank : bool, optional
+ if True, return the effective rank of the matrix
+
+ Returns
+ -------
+ B : array, shape (N, N)
+ The pseudo-inverse of matrix `a`.
+ rank : int
+ The effective rank of the matrix. Returned if return_rank == True
+
+ Raises
+ ------
+ LinAlgError
+ If eigenvalue does not converge
+
+ Examples
+ --------
+ >>> from numpy import *
+ >>> a = random.randn(9, 6)
+ >>> a = np.dot(a, a.T)
+ >>> B = pinvh(a)
+ >>> allclose(a, dot(a, dot(B, a)))
+ True
+ >>> allclose(B, dot(B, dot(a, B)))
+ True
+ """
+ a = np.asarray_chkfinite(a)
+ s, u = decomp.eigh(a, lower=lower)
+
+ if rcond is not None:
+ cond = rcond
+ if cond in [None, -1]:
+ t = u.dtype.char.lower()
+ factor = {'f': 1E3, 'd': 1E6}
+ cond = factor[t] * np.finfo(t).eps
+
+ # For hermitian matrices, singular values equal abs(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]
+
+ B = np.dot(u * psigma_diag, np.conjugate(u).T)
+
+ if return_rank:
+ return B, np.sum(above_cutoff)
+ else:
+ return B
View
56 scipy/linalg/tests/test_basic.py
@@ -28,7 +28,7 @@
assert_equal, assert_almost_equal, assert_array_almost_equal, assert_, \
assert_allclose
-from scipy.linalg import solve, inv, det, lstsq, pinv, pinv2, norm,\
+from scipy.linalg import solve, inv, det, lstsq, pinv, pinv2, pinvh, norm,\
solve_banded, solveh_banded, solve_triangular
from scipy.linalg._testutils import assert_no_overwrite
@@ -533,30 +533,66 @@ def test_random_complex_overdet(self):
class TestPinv(TestCase):
- def test_simple(self):
- a=array([[1,2,3],[4,5,6.],[7,8,10]])
+ def test_simple_real(self):
+ a = array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float)
a_pinv = pinv(a)
- assert_array_almost_equal(dot(a,a_pinv),[[1,0,0],[0,1,0],[0,0,1]])
+ assert_array_almost_equal(dot(a,a_pinv), np.eye(3))
a_pinv = pinv2(a)
- assert_array_almost_equal(dot(a,a_pinv),[[1,0,0],[0,1,0],[0,0,1]])
- def test_simple_0det(self):
- a=array([[1,2,3],[4,5,6.],[7,8,9]])
+ 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, 10]], dtype=float)
+ + 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]], dtype=float))
+ a_pinv = pinv(a)
+ assert_array_almost_equal(dot(a, a_pinv), np.eye(3))
+ a_pinv = pinv2(a)
+ assert_array_almost_equal(dot(a, a_pinv), np.eye(3))
+
+ def test_simple_singular(self):
+ a = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)
a_pinv = pinv(a)
a_pinv2 = pinv2(a)
assert_array_almost_equal(a_pinv,a_pinv2)
def test_simple_cols(self):
- a=array([[1,2,3],[4,5,6.]])
+ a = array([[1, 2, 3], [4, 5, 6]], dtype=float)
a_pinv = pinv(a)
a_pinv2 = pinv2(a)
assert_array_almost_equal(a_pinv,a_pinv2)
def test_simple_rows(self):
- a=array([[1,2],[3,4],[5,6]])
+ a = array([[1, 2], [3, 4], [5, 6]], dtype=float)
a_pinv = pinv(a)
a_pinv2 = pinv2(a)
assert_array_almost_equal(a_pinv,a_pinv2)
+
+class TestPinvSymmetric(TestCase):
+
+ def test_simple_real(self):
+ a = array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float)
+ 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]], dtype=float)
+ 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]], dtype=float)
+ + 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]], dtype=float))
+ a = np.dot(a, a.conj().T)
+ a_pinv = pinvh(a)
+ assert_array_almost_equal(np.dot(a, a_pinv), np.eye(3))
+
+
class TestNorm(object):
def test_types(self):
@@ -614,6 +650,8 @@ def test_pinv(self):
assert_no_overwrite(pinv, [(3,3)])
def test_pinv2(self):
assert_no_overwrite(pinv2, [(3,3)])
+ def test_pinvh(self):
+ assert_no_overwrite(pinvh, [(3,3)])
if __name__ == "__main__":
run_module_suite()

0 comments on commit 5adeb12

Please sign in to comment.