[Scipy-svn] r6245 - in trunk/scipy/sparse/linalg/eigen/arpack: . tests

[scipy-svn at scipy.org]

Author: cdavid
Date: 2010-02-22 03:02:18 -0600 (Mon, 22 Feb 2010)
New Revision: 6245

Modified:
   trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
   trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py
Log:
ENH: implement naive SVD for sparse real matrices.

Modified: trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-02-21 21:04:04 UTC (rev 6244)
+++ trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-02-22 09:02:18 UTC (rev 6245)
@@ -39,13 +39,14 @@
 
 __docformat__ = "restructuredtext en"
 
-__all___=['eigen','eigen_symmetric']
+__all___=['eigen','eigen_symmetric', 'svd']
 
 import warnings
 
 import _arpack
 import numpy as np
 from scipy.sparse.linalg.interface import aslinearoperator
+from scipy.sparse import csc_matrix, csr_matrix
 
 _type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
 _ndigits = {'f':5, 'd':12, 'F':5, 'D':12}
@@ -488,3 +489,58 @@
     if return_eigenvectors:
         return d,z
     return d
+
+def svd(A, k=6):
+    """Compute a few singular values/vectors for a sparse matrix using ARPACK.
+
+    Parameters
+    ----------
+    A: sparse matrix
+        Array to compute the SVD on.
+    k: int
+        Number of singular values and vectors to compute.
+
+    Note
+    ----
+    This is a naive implementation using the symmetric eigensolver on A.T * A
+    or A * A.T, depending on which one is more efficient.
+
+    Complex support is not implemented yet
+    """
+    # TODO: implement complex support once ARPACK-based eigen_hermitian is
+    # available
+    n, m = A.shape
+
+    if np.iscomplexobj(A):
+        raise NotImplementedError("Complex support for sparse SVD not " \
+                                  "implemented yet")
+        op = lambda x: x.T.conjugate()
+    else:
+        op = lambda x: x.T
+
+    def _left(x):
+        x = csc_matrix(x)
+        m = op(x) * x
+
+        eigvals, eigvec = eigen_symmetric(m, k)
+        s = np.sqrt(eigvals)
+
+        v = eigvec
+        u = (x * v) / s
+        return u, s, op(v)
+
+    def _right(x):
+        x = csr_matrix(x)
+        m = x * op(x)
+
+        eigvals, eigvec = eigen_symmetric(m, k)
+        s = np.sqrt(eigvals)
+
+        u = eigvec
+        vh = (op(u) * x) / s[:, None]
+        return u, s, vh
+
+    if n > m:
+        return _left(A)
+    else:
+        return _right(A)

Modified: trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py	2010-02-21 21:04:04 UTC (rev 6244)
+++ trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py	2010-02-22 09:02:18 UTC (rev 6245)
@@ -4,12 +4,15 @@
   python tests/test_arpack.py [-l<int>] [-v<int>]
 
 """
+import numpy as np
 
 from numpy.testing import *
 
 from numpy import array, finfo, argsort, dot, round, conj, random
-from scipy.sparse.linalg.eigen.arpack import eigen_symmetric, eigen
+from scipy.sparse.linalg.eigen.arpack import eigen_symmetric, eigen, svd
 
+from scipy.linalg import svd as dsvd
+
 def assert_almost_equal_cc(actual,desired,decimal=7,err_msg='',verbose=True):
     # almost equal or complex conjugates almost equal
     try:
@@ -263,5 +266,51 @@
                 for m in self.nonsymmetric:
                     self.eval_evec(m,typ,k,which)
 
+def sorted_svd(m, k):
+    """Compute svd of a dense matrix m, and return singular vectors/values
+    sorted."""
+    u, s, vh = dsvd(m)
+    ii = np.argsort(s)[-k:]
+
+    return u[:, ii], s[ii], vh[ii]
+
+def svd_estimate(u, s, vh):
+    return np.dot(u, np.dot(np.diag(s), vh))
+
+class TestSparseSvd(TestCase):
+    def test_simple_real(self):
+        x = np.array([[1, 2, 3],
+                      [3, 4, 3],
+                      [1, 0, 2],
+                      [0, 0, 1]], np.float)
+
+        for m in [x.T, x]:
+            for k in range(1, 3):
+                u, s, vh = sorted_svd(m, k)
+                su, ss, svh = svd(m, k)
+
+                m_hat = svd_estimate(u, s, vh)
+                sm_hat = svd_estimate(su, ss, svh)
+
+                assert_array_almost_equal_nulp(m_hat, sm_hat, nulp=1000)
+
+    @dec.knownfailureif(True, "Complex sparse SVD not implemented (depends on "
+                              "Hermitian support in eigen_symmetric")
+    def test_simple_complex(self):
+        x = np.array([[1, 2, 3],
+                      [3, 4, 3],
+                      [1+1j, 0, 2],
+                      [0, 0, 1]], np.complex)
+
+        for m in [x, x.T.conjugate()]:
+            for k in range(1, 3):
+                u, s, vh = sorted_svd(m, k)
+                su, ss, svh = svd(m, k)
+
+                m_hat = svd_estimate(u, s, vh)
+                sm_hat = svd_estimate(su, ss, svh)
+
+                assert_array_almost_equal_nulp(m_hat, sm_hat, nulp=1000)
+
 if __name__ == "__main__":
     run_module_suite()