Python Interface

primme.svds(A, k=6, ncv=None, tol=0, which='LM', v0=None, maxiter=None, return_singular_vectors=True, precAHA=None, precAAH=None, precAug=None, u0=None, orthou0=None, orthov0=None, return_stats=False, maxBlockSize=0, method=None, methodStage1=None, methodStage2=None, return_history=False, convtest=None, **kargs)

Compute k singular values and vectors of the matrix A.

Parameters
  • A ({sparse matrix, LinearOperator}) – Array to compute the SVD on, of shape (M, N)

  • k (int, optional) – Number of singular values and vectors to compute. Must be 1 <= k < min(A.shape).

  • ncv (int, optional) – The maximum size of the basis

  • tol (float, optional) –

    Tolerance for singular values. Zero (default) means 10**4 times the machine precision.

    A triplet (u,sigma,v) is marked as converged when (||A*v - sigma*u||**2 + ||A.H*u - sigma*v||**2)**.5 is less than “tol” * ||A||, or close to the minimum tolerance that the method can achieve. See the note.

    The value is ignored if convtest is provided.

  • which (str ['LM' | 'SM'] or number, optional) –

    Which k singular values to find:

    • ’LM’ : largest singular values

    • ’SM’ : smallest singular values

    • number : closest singular values to (referred as sigma later)

  • u0 (ndarray, optional) –

    Initial guesses for the left singular vectors.

    If only u0 or v0 is provided, the other is computed. If both are provided, u0 and v0 should have the same number of columns.

  • v0 (ndarray, optional) – Initial guesses for the right singular vectors.

  • maxiter (int, optional) – Maximum number of matvecs with A and A.H.

  • precAHA ({N x N matrix, array, sparse matrix, LinearOperator}, optional) – Approximate inverse of (A.H*A - sigma**2*I). If provided and M>=N, it usually accelerates the convergence.

  • precAAH ({M x M matrix, array, sparse matrix, LinearOperator}, optional) – Approximate inverse of (A*A.H - sigma**2*I). If provided and M<N, it usually accelerates the convergence.

  • precAug ({(M+N) x (M+N) matrix, array, sparse matrix, LinearOperator}, optional) – Approximate inverse of ([zeros() A.H; zeros() A] - sigma*I).

  • orthou0 (ndarray, optional) –

    Left orthogonal vector constrain.

    Seek singular triplets orthogonal to orthou0 and orthov0. The provided vectors should be orthonormal. If only orthou0 or orthov0 is provided, the other is computed. Useful to avoid converging to previously computed solutions.

  • orthov0 (ndarray, optional) – Right orthogonal vector constrain. See orthou0.

  • maxBlockSize (int, optional) – Maximum number of vectors added at every iteration.

  • convtest (callable) –

    User-defined function to mark an approximate singular triplet as converged.

    The function is called as convtest(sval, svecleft, svecright, resNorm) and returns True if the triplet with value sval, left vector svecleft, right vector svecright, and residual norm resNorm is considered converged.

  • return_stats (bool, optional) – If True, the function returns extra information (see stats in Returns).

  • return_history (bool, optional) – If True, the function returns performance information at every iteration

Returns

  • u (ndarray, shape=(M, k), optional) – Unitary matrix having left singular vectors as columns. Returned if return_singular_vectors is True.

  • s (ndarray, shape=(k,)) – The singular values.

  • vt (ndarray, shape=(k, N), optional) – Unitary matrix having right singular vectors as rows. Returned if return_singular_vectors is True.

  • stats (dict, optional (if return_stats)) – Extra information reported by PRIMME:

    • ”numOuterIterations”: number of outer iterations

    • ”numRestarts”: number of restarts

    • ”numMatvecs”: number of matvecs with A and A.H

    • ”numPreconds”: cumulative number of applications of precAHA, precAAH and precAug

    • ”elapsedTime”: time that took

    • ”rnorms” : (||A*v[:,i] - sigma[i]*u[:,i]||**2 + ||A.H*u[:,i] - sigma[i]*v[:,i]||**2)**.5

    • ”hist” : (if return_history) report at every outer iteration of:

      • ”elapsedTime”: time spent up to now

      • ”numMatvecs”: number of A*v and A.H*v spent up to now

      • ”nconv”: number of converged triplets

      • ”sval”: singular value of the first unconverged triplet

      • ”resNorm”: residual norm of the first unconverged triplet

Notes

The default method used is the hybrid method, which first solves the equivalent eigenvalue problem A.H*A or A*A.H (normal equations) and then refines the solution solving the augmented problem. The minimum tolerance that this method can achieve is ||A||*epsilon, where epsilon is the machine precision. However it may not return triplets with singular values smaller than ||A||*epsilon if “tol” is smaller than ||A||*epsilon/sigma.

This function is a wrapper to PRIMME functions to find singular values and vectors 1.

References

1

PRIMME Software, https://github.com/primme/primme

2

L. Wu, E. Romero and A. Stathopoulos, PRIMME_SVDS: A High- Performance Preconditioned SVD Solver for Accurate Large-Scale Computations. https://arxiv.org/abs/1607.01404

See also

primme.eigsh()

eigenvalue decomposition for a sparse symmetrix/complex Hermitian matrix A

scipy.sparse.linalg.eigs()

eigenvalues and eigenvectors for a general (nonsymmetric) matrix A

Examples

>>> import primme, scipy.sparse
>>> A = scipy.sparse.spdiags(range(1, 11), [0], 100, 10) # sparse diag. rect. matrix
>>> svecs_left, svals, svecs_right = primme.svds(A, 3, tol=1e-6, which='LM')
>>> svals # the three largest singular values of A
array([10., 9., 8.])
>>> import primme, scipy.sparse, numpy as np
>>> A = scipy.sparse.rand(10000, 100, random_state=10)
>>> prec = scipy.sparse.spdiags(np.reciprocal(A.multiply(A).sum(axis=0)),
...           [0], 100, 100) # square diag. preconditioner
>>> # the three smallest singular values of A, using preconditioning
>>> svecs_left, svals, svecs_right = primme.svds(A, 3, which='SM', tol=1e-6, precAHA=prec)
>>> ["%.5f" % x for x in svals.flat] 
['4.57263', '4.78752', '4.82229']
>>> # Giving the matvecs as functions
>>> import primme, scipy.sparse, numpy as np
>>> Bdiag = np.arange(0, 100).reshape((100,1))
>>> Bdiagr = np.concatenate((np.arange(0, 100).reshape((100,1)).astype(np.float32), np.zeros((100,1), dtype=np.float32)), axis=None).reshape((200,1))
>>> def Bmatmat(x):
...    if len(x.shape) == 1: x = x.reshape((100,1))
...    return np.vstack((Bdiag * x, np.zeros((100, x.shape[1]), dtype=np.float32)))
...
>>> def Brmatmat(x):
...    if len(x.shape) == 1: x = x.reshape((200,1))
...    return (Bdiagr * x)[0:100,:]
...
>>> B = scipy.sparse.linalg.LinearOperator((200,100), matvec=Bmatmat, matmat=Bmatmat, rmatvec=Brmatmat, dtype=np.float32)
>>> svecs_left, svals, svecs_right = primme.svds(B, 5, which='LM', tol=1e-6)
>>> svals 
array([99., 98., 97., 96., 95.])