# Pinv enhancements #289

Closed
wants to merge 7 commits into
from
Commits
Failed to load files and symbols.
+152 −28
 @@ -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
@@ -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
@@ -500,7 +501,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.
@@ -515,11 +516,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
------
@@ -540,10 +545,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.
@@ -560,11 +571,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
------
@@ -582,21 +597,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}

#### charris Aug 16, 2012

Member

I wonder if the updated version in np.matrix_rank would be a better choice for the default condition? See the discussion at http://thread.gmane.org/gmane.comp.python.numeric.general/50396/focus=50912

#### jakevdp Aug 16, 2012

Member

Interesting discussion: here I just duplicated the previous behavior. Any change would be a potential compatibility issue: before making that decision we should probably discuss it on the mailing list.

+ 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
 @@ -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 @@ -534,30 +534,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): @@ -608,6 +644,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()