[Scipy-svn] r6273 - trunk/scipy/sparse/linalg/eigen/arpack

scipy-svn@scip... scipy-svn@scip...
Fri Mar 26 00:35:35 CDT 2010

```Author: cdavid
Date: 2010-03-26 00:35:35 -0500 (Fri, 26 Mar 2010)
New Revision: 6273

Modified:
trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
Log:
ENH: use implicit A'A computation for sparse SVD.

Modified: trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-03-26 05:35:25 UTC (rev 6272)
+++ trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-03-26 05:35:35 UTC (rev 6273)
@@ -519,22 +519,35 @@
else:
op = lambda x: x.T

-    def _left(x):
+    tp = A.dtype.char
+    linear_at = aslinearoperator(op(A))
+    linear_a = aslinearoperator(A)
+
+    def _left(x, sz):
x = csc_matrix(x)
-        m = op(x) * x

-        eigvals, eigvec = eigen_symmetric(m, k)
+        matvec = lambda x: linear_at.matvec(linear_a.matvec(x))
+        params = _SymmetricArpackParams(sz, k, tp, matvec)
+
+        while not params.converged:
+            params.iterate()
+        eigvals, eigvec = params.extract(True)
s = np.sqrt(eigvals)

v = eigvec
u = (x * v) / s
return u, s, op(v)

-    def _right(x):
+    def _right(x, sz):
x = csr_matrix(x)
-        m = x * op(x)

-        eigvals, eigvec = eigen_symmetric(m, k)
+        matvec = lambda x: linear_a.matvec(linear_at.matvec(x))
+        params = _SymmetricArpackParams(sz, k, tp, matvec)
+
+        while not params.converged:
+            params.iterate()
+        eigvals, eigvec = params.extract(True)
+
s = np.sqrt(eigvals)

u = eigvec
@@ -542,6 +555,6 @@
return u, s, vh

if n > m:
-        return _left(A)
+        return _left(A, m)
else:
-        return _right(A)
+        return _right(A, n)

```