[Scipy-svn] r6992 - in trunk/scipy/sparse/linalg/eigen/arpack: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Dec 4 15:25:55 EST 2010
Author: ptvirtan
Date: 2010-12-04 14:25:55 -0600 (Sat, 04 Dec 2010)
New Revision: 6992
Modified:
trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py
Log:
ENH: sparse/arpack: rewrite the sparse svd routine to handle complex matrices
The ARPACK user guide implies that the only way to compute eigenvalues
of Hermitian matrices is to use zneupd, as there is no special support
for hermitian matrices.
Modified: trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/arpack.py 2010-12-04 20:25:44 UTC (rev 6991)
+++ trunk/scipy/sparse/linalg/eigen/arpack/arpack.py 2010-12-04 20:25:55 UTC (rev 6992)
@@ -43,8 +43,8 @@
import _arpack
import numpy as np
-from scipy.sparse.linalg.interface import aslinearoperator
-from scipy.sparse import csc_matrix, csr_matrix
+from scipy.sparse.linalg.interface import aslinearoperator, LinearOperator
+from scipy.sparse import csc_matrix, csr_matrix, isspmatrix
_type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
_ndigits = {'f':5, 'd':12, 'F':5, 'D':12}
@@ -656,70 +656,62 @@
return params.extract(return_eigenvectors)
-def svds(A, k=6):
+def svds(A, k=6, ncv=None, tol=0):
"""Compute a few singular values/vectors for a sparse matrix using ARPACK.
Parameters
----------
- A: sparse matrix
- Array to compute the SVD on.
- k: int
+ A : sparse matrix
+ Array to compute the SVD on
+ k : int, optional
Number of singular values and vectors to compute.
+ ncv : integer
+ The number of Lanczos vectors generated
+ ncv must be greater than k+1 and smaller than n;
+ it is recommended that ncv > 2*k
+ tol : float, optional
+ Tolerance for singular values. Zero (default) means machine precision.
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.
+ This is a naive implementation using an eigensolver on A.H * A or
+ A * A.H, depending on which one is more efficient.
- Complex support is not implemented yet
"""
- # TODO: implement complex support once ARPACK-based eigen_hermitian is
- # available
+ if not (isinstance(A, np.ndarray) or isspmatrix(A)):
+ A = np.asarray(A)
+
n, m = A.shape
- if np.iscomplexobj(A):
- raise NotImplementedError("Complex support for sparse SVD not "
- "implemented yet")
- op = lambda x: x.T.conjugate()
+ if np.issubdtype(A.dtype, np.complexfloating):
+ herm = lambda x: x.T.conjugate()
+ eigensolver = eigs
else:
- op = lambda x: x.T
+ herm = lambda x: x.T
+ eigensolver = eigsh
- tp = A.dtype.char
- linear_at = aslinearoperator(op(A))
- linear_a = aslinearoperator(A)
+ if n > m:
+ X = A
+ XH = herm(A)
+ else:
+ XH = A
+ X = herm(A)
- def _left(x, sz):
- x = csc_matrix(x)
+ def matvec_XH_X(x):
+ return XH.dot(X.dot(x))
- matvec = lambda x: linear_at.matvec(linear_a.matvec(x))
- params = _SymmetricArpackParams(sz, k, tp, matvec)
+ XH_X = LinearOperator(matvec=matvec_XH_X, dtype=X.dtype,
+ shape=(X.shape[1], X.shape[1]))
- while not params.converged:
- params.iterate()
- eigvals, eigvec = params.extract(True)
- s = np.sqrt(eigvals)
+ eigvals, eigvec = eigensolver(XH_X, k=k, tol=tol**2)
+ s = np.sqrt(eigvals)
+ if n > m:
v = eigvec
- u = (x * v) / s
- return u, s, op(v)
-
- def _right(x, sz):
- x = csr_matrix(x)
-
- 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 = X.dot(v) / s
+ vh = herm(v)
+ else:
u = eigvec
- vh = (op(u) * x) / s[:, None]
- return u, s, vh
+ vh = herm(X.dot(u) / s)
- if n > m:
- return _left(A, m)
- else:
- return _right(A, n)
+ return u, s, vh
Modified: trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py 2010-12-04 20:25:44 UTC (rev 6991)
+++ trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py 2010-12-04 20:25:55 UTC (rev 6992)
@@ -356,6 +356,8 @@
def sorted_svd(m, k):
"""Compute svd of a dense matrix m, and return singular vectors/values
sorted."""
+ if isspmatrix(m):
+ m = m.todense()
u, s, vh = svd(m)
ii = np.argsort(s)[-k:]
@@ -370,9 +372,14 @@
[3, 4, 3],
[1, 0, 2],
[0, 0, 1]], np.float)
+ y = np.array([[1, 2, 3, 8],
+ [3, 4, 3, 5],
+ [1, 0, 2, 3],
+ [0, 0, 1, 0]], np.float)
+ z = csc_matrix(x)
- for m in [x.T, x]:
- for k in range(1, 3):
+ for m in [x.T, x, y, z, z.T]:
+ for k in range(1, min(m.shape)):
u, s, vh = sorted_svd(m, k)
su, ss, svh = svds(m, k)
@@ -381,16 +388,19 @@
assert_array_almost_equal_nulp(m_hat, sm_hat, nulp=1000)
- @dec.knownfailureif(True, "Complex sparse SVD not implemented (depends on "
- "Hermitian support in eigsh")
def test_simple_complex(self):
x = np.array([[1, 2, 3],
[3, 4, 3],
[1+1j, 0, 2],
[0, 0, 1]], np.complex)
+ y = np.array([[1, 2, 3, 8+5j],
+ [3-2j, 4, 3, 5],
+ [1, 0, 2, 3],
+ [0, 0, 1, 0]], np.complex)
+ z = csc_matrix(x)
- for m in [x, x.T.conjugate()]:
- for k in range(1, 3):
+ for m in [x, x.T.conjugate(), x.T, y, y.conjugate(), z, z.T]:
+ for k in range(1, min(m.shape)-1):
u, s, vh = sorted_svd(m, k)
su, ss, svh = svds(m, k)
More information about the Scipy-svn
mailing list