Compute pseudoinverse using eigendecomposition

vene · amueller · commit 0946f5945968 · 2012-09-04T20:53:58.000+01:00
diff --git a/sklearn/covariance/empirical_covariance_.py b/sklearn/covariance/empirical_covariance_.py
@@ -17,7 +17,7 @@
 
 from ..base import BaseEstimator
 from ..utils import array2d
-from ..utils.extmath import fast_logdet
+from ..utils.extmath import fast_logdet, symmetric_pinv
 
 
 def log_likelihood(emp_cov, precision):
@@ -113,7 +113,7 @@ def _set_covariance(self, covariance):
         self.covariance_ = covariance
         # set precision
         if self.store_precision:
-            self.precision_ = linalg.pinv(covariance)
+            self.precision_ = symmetric_pinv(covariance)
         else:
             self.precision_ = None
 
@@ -129,7 +129,7 @@ def get_precision(self):
         if self.store_precision:
             precision = self.precision_
         else:
-            precision = linalg.pinv(self.covariance_)
+            precision = symmetric_pinv(self.covariance_)
         return precision
 
     def fit(self, X):
diff --git a/sklearn/covariance/graph_lasso_.py b/sklearn/covariance/graph_lasso_.py
@@ -17,6 +17,7 @@
                 EmpiricalCovariance, log_likelihood
 
 from ..utils import ConvergenceWarning
+from ..utils.extmath import symmetric_pinv
 from ..linear_model import lars_path
 from ..linear_model import cd_fast
 from ..cross_validation import check_cv, cross_val_score
@@ -143,7 +144,7 @@ def graph_lasso(emp_cov, alpha, cov_init=None, mode='cd', tol=1e-4,
     covariance_ *= 0.95
     diagonal = emp_cov.flat[::n_features + 1]
     covariance_.flat[::n_features + 1] = diagonal
-    precision_ = linalg.pinv(covariance_)
+    precision_ = symmetric_pinv(covariance_)
 
     indices = np.arange(n_features)
     costs = list()
diff --git a/sklearn/covariance/robust_covariance.py b/sklearn/covariance/robust_covariance.py
@@ -13,7 +13,7 @@
 from scipy.stats import chi2
 
 from . import empirical_covariance, EmpiricalCovariance
-from ..utils.extmath import fast_logdet
+from ..utils.extmath import fast_logdet, symmetric_pinv
 from ..utils import check_random_state
 
 
@@ -85,7 +85,7 @@ def c_step(X, n_support, remaining_iterations=30, initial_estimates=None,
         location = initial_estimates[0]
         covariance = initial_estimates[1]
         # run a special iteration for that case (to get an initial support)
-        precision = linalg.pinv(covariance)
+        precision = symmetric_pinv(covariance)
         X_centered = X - location
         dist = (np.dot(X_centered, precision) * X_centered).sum(1)
         # compute new estimates
@@ -98,15 +98,15 @@ def c_step(X, n_support, remaining_iterations=30, initial_estimates=None,
     # Iterative procedure for Minimum Covariance Determinant computation
     det = fast_logdet(covariance)
     while (det < previous_det) and (remaining_iterations > 0):
-        # compute a new support from the full data set mahalanobis distances
-        precision = linalg.pinv(covariance)
-        X_centered = X - location
-        dist = (np.dot(X_centered, precision) * X_centered).sum(axis=1)
         # save old estimates values
         previous_location = location
         previous_covariance = covariance
         previous_det = det
         previous_support = support
+        # compute a new support from the full data set mahalanobis distances
+        precision = symmetric_pinv(covariance)
+        X_centered = X - location
+        dist = (np.dot(X_centered, precision) * X_centered).sum(axis=1)
         # compute new estimates
         support = np.zeros(n_samples).astype(bool)
         support[np.argsort(dist)[:n_support]] = True
@@ -343,7 +343,8 @@ def fast_mcd(X, support_fraction=None,
         support[np.argsort(np.abs(X - location), axis=0)[:n_support]] = True
         covariance = np.asarray([[np.var(X[support])]])
         location = np.array([location])
-        dist = (np.dot(X_centered, linalg.pinv(covariance)) \
+        precision = symmetric_pinv(covariance)
+        dist = (np.dot(X_centered, precision) \
                     * (X_centered)).sum(axis=1)
 
     ### Starting FastMCD algorithm for p-dimensional case
@@ -542,7 +543,8 @@ def fit(self, X):
             raw_location = np.zeros(n_features)
             raw_covariance = self._nonrobust_covariance(
                     X[raw_support], assume_centered=True)
-            raw_dist = np.sum(np.dot(X, linalg.pinv(raw_covariance)) * X, 1)
+            precision = symmetric_pinv(raw_covariance)
+            raw_dist = np.sum(np.dot(X, precision) * X, 1)
         self.raw_location_ = raw_location
         self.raw_covariance_ = raw_covariance
         self.raw_support_ = raw_support
diff --git a/sklearn/utils/extmath.py b/sklearn/utils/extmath.py
@@ -297,3 +297,65 @@ def weighted_mode(a, w, axis=0):
         oldcounts = np.maximum(counts, oldcounts)
         oldmostfreq = mostfrequent
     return mostfrequent, oldcounts
+
+
+def symmetric_pinv(a, cond=None, rcond=None):
+    """Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+    Calculate a generalized inverse of a symmetric matrix using its
+    eigenvalue decomposition and including all 'large' eigenvalues.
+
+    Inspired by ``scipy.linalg.pinv2``, credited to Pearu Peterson and Travis
+    Oliphant.
+
+    Parameters
+    ----------
+    a : array, shape (N, N)
+        Symmetric 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.
+
+    Returns
+    -------
+    B : array, shape (N, N)
+
+    Raises LinAlgError if eigenvalue does not converge
+
+    Examples
+    --------
+    >>> from numpy import *
+    >>> a = random.randn(9, 6)
+    >>> a = np.dot(a, a.T)
+    >>> B = symmetric_pinv(a)
+    >>> allclose(a, dot(a, dot(B, a)))
+    True
+    >>> allclose(B, dot(B, dot(a, B)))
+    True
+
+    """
+    a = np.asarray_chkfinite(a)
+    s, u = linalg.eigh(a)
+    # eigh returns eigvals in reverse order, but this doesn't affect anything.
+
+    t = u.dtype.char
+    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]]
+    n = a.shape[0]
+    cutoff = cond * np.maximum.reduce(s)
+    psigma = np.zeros((n, 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
+    #XXX: above comment is from scipy, but I (@vene)'ll take a look
+    return np.transpose(np.conjugate(np.dot(np.dot(u, psigma),
+                                     u.T.conjugate())))
diff --git a/sklearn/utils/tests/test_utils.py b/sklearn/utils/tests/test_utils.py
@@ -1,13 +1,17 @@
-from nose.tools import assert_equal, assert_raises, assert_true
 import warnings
 
 import numpy as np
 import scipy.sparse as sp
+from scipy.linalg import pinv2
+
+from nose.tools import assert_equal, assert_raises, assert_true
+from numpy.testing import assert_almost_equal
 
 from sklearn.utils import check_random_state
 from sklearn.utils import deprecated
 from sklearn.utils import resample
 from sklearn.utils import safe_mask
+from sklearn.utils.extmath import symmetric_pinv
 
 
 def test_make_rng():
@@ -87,3 +91,9 @@ def test_safe_mask():
 
     mask = safe_mask(X_csr, mask)
     assert_equal(X_csr[mask].shape[0], 3)
+
+
+def test_symmetric_pinv():
+    a = np.random.randn(5, 3)
+    a = np.dot(a, a.T)  # symmetric singular matrix
+    assert_almost_equal(pinv2(a), symmetric_pinv(a))
