Python Interface

primme.eigsh(A, k=6, M=None, sigma=None, which='LM', v0=None, ncv=None, maxiter=None, tol=0, return_eigenvectors=True, Minv=None, OPinv=None, mode='normal', lock=None, return_stats=False, maxBlockSize=0, minRestartSize=0, maxPrevRetain=0, method=None, return_history=False, convtest=None, **kargs)

Find k eigenvalues and eigenvectors of the real symmetric square matrix or complex Hermitian matrix A.

Solves A * x[i] = w[i] * x[i], the standard eigenvalue problem for w[i] eigenvalues with corresponding eigenvectors x[i].

If M is specified, solves A * x[i] = w[i] * M * x[i], the generalized eigenvalue problem for w[i] eigenvalues with corresponding eigenvectors x[i]

Parameters

  • A (An N x N matrix, array, sparse matrix, or LinearOperator) – the operation A * x, where A is a real symmetric matrix or complex Hermitian.

  • k (int, optional) – The number of eigenvalues and eigenvectors to be computed. Must be 1 <= k < min(A.shape).

  • M (An N x N matrix, array, sparse matrix, or LinearOperator) –

    the operation M * x for the generalized eigenvalue problem

    A * x = w * M * x.

    M must represent a real, symmetric matrix if A is real, and must represent a complex, Hermitian matrix if A is complex. For best results, the data type of M should be the same as that of A.

  • sigma (real, optional) – Find eigenvalues near sigma.

  • v0 (N x i, ndarray, optional) – Initial guesses to the eigenvectors.

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

  • which (str ['LM' | 'SM' | 'LA' | 'SA' | number]) –

    Which k eigenvectors and eigenvalues to find:

    ’LM’ : Largest in magnitude eigenvalues; the farthest from sigma

    ’SM’ : Smallest in magnitude eigenvalues; the closest to sigma

    ’LA’ : Largest algebraic eigenvalues

    ’SA’ : Smallest algebraic eigenvalues

    ’CLT’ : closest but left to sigma

    ’CGT’ : closest but greater than sigma

    number : the closest to which

    When sigma == None, ‘LM’, ‘SM’, ‘CLT’, and ‘CGT’ treat sigma as zero.

  • maxiter (int, optional) – Maximum number of iterations.

  • tol (float) –

    Tolerance for eigenpairs (stopping criterion). The default value is sqrt of machine precision.

    An eigenpair (lamba,v) is marked as converged when ||A*v - lambda*B*v|| < max(|eig(A,B)|)*tol.

    The value is ignored if convtest is provided.

  • Minv ((not supported yet)) – The inverse of M in the generalized eigenproblem.

  • OPinv (N x N matrix, array, sparse matrix, or LinearOperator, optional) – Preconditioner to accelerate the convergence. Usually it is an approximation of the inverse of (A - sigma*M).

  • return_eigenvectors (bool, optional) – Return eigenvectors (True) in addition to eigenvalues

  • mode (string ['normal' | 'buckling' | 'cayley']) – Only ‘normal’ mode is supported.

  • lock (N x i, ndarray, optional) – Seek the eigenvectors orthogonal to these ones. The provided vectors should be orthonormal. Useful to avoid converging to previously computed solutions.

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

  • minRestartSize (int, optional) – Number of approximate eigenvectors kept during restart.

  • maxPrevRetain (int, optional) – Number of approximate eigenvectors kept from previous iteration in restart. Also referred as +k vectors in GD+k.

  • method (int, optional) –

    Preset method, one of:

    • DEFAULT_MIN_TIME : a variant of JDQMR,

    • DEFAULT_MIN_MATVECS : GD+k

    • DYNAMIC : choose dynamically between these previous methods.

    See a detailed description of the methods and other possible values in 2.

  • convtest (callable) –

    User-defined function to mark an approximate eigenpair as converged.

    The function is called as convtest(eval, evec, resNorm) and returns True if the eigenpair with value eval, vector evec 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 (see hist in Returns).

Returns

  • w (array) – Array of k eigenvalues ordered to best satisfy “which”.

  • v (array) – An array representing the k eigenvectors. The column v[:, i] is the eigenvector corresponding to the eigenvalue w[i].

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

    • ”numOuterIterations”: number of outer iterations

    • ”numRestarts”: number of restarts

    • ”numMatvecs”: number of A*v

    • ”numPreconds”: number of OPinv*v

    • ”elapsedTime”: time that took

    • ”estimateMinEVal”: the leftmost Ritz value seen

    • ”estimateMaxEVal”: the rightmost Ritz value seen

    • ”estimateLargestSVal”: the largest singular value seen

    • ”rnorms” : ||A*x[i] - x[i]*w[i]||

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

      • ”elapsedTime”: time spent up to now

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

      • ”nconv”: number of converged pair

      • ”eval”: eigenvalue of the first unconverged pair

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

Raises

PrimmeError – When the requested convergence is not obtained. The PRIMME error code can be found as err attribute of the exception object.

See also

scipy.sparse.linalg.eigs()

eigenvalues and eigenvectors for a general (nonsymmetric) matrix A

primme.svds()

singular value decomposition for a matrix A

Notes

This function is a wrapper to PRIMME functions to find the eigenvalues and eigenvectors 1.

References

1

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

2

Preset Methods, http://www.cs.wm.edu/~andreas/software/doc/readme.html#preset-methods

3

A. Stathopoulos and J. R. McCombs PRIMME: PReconditioned Iterative MultiMethod Eigensolver: Methods and software description, ACM Transaction on Mathematical Software Vol. 37, No. 2, (2010), 21:1-21:30.

Examples

>>> import primme, scipy.sparse
>>> A = scipy.sparse.spdiags(range(100), [0], 100, 100) # sparse diag. matrix
>>> evals, evecs = primme.eigsh(A, 3, tol=1e-6, which='LA')
>>> evals # the three largest eigenvalues of A
array([99., 98., 97.])
>>> new_evals, new_evecs = primme.eigsh(A, 3, tol=1e-6, which='LA', lock=evecs)
>>> new_evals # the next three largest eigenvalues
array([96., 95., 94.])
>>> evals, evecs = primme.eigsh(A, 3, tol=1e-6, which=50.1)
>>> evals # the three closest eigenvalues to 50.1
array([50.,  51.,  49.])
>>> M = scipy.sparse.spdiags(np.asarray(range(99,-1,-1)), [0], 100, 100)
>>> # the smallest eigenvalues of the eigenproblem (A,M)
>>> evals, evecs = primme.eigsh(A, 3, M=M, tol=1e-6, which='SA')
>>> evals 
array([1.0035e-07, 1.0204e-02, 2.0618e-02])

>>> # Giving the matvec as a function
>>> import primme, scipy.sparse, numpy as np
>>> Adiag = np.arange(0, 100).reshape((100,1))
>>> def Amatmat(x):
...    if len(x.shape) == 1: x = x.reshape((100,1))
...    return Adiag * x   # equivalent to diag(Adiag).dot(x)
...
>>> A = scipy.sparse.linalg.LinearOperator((100,100), matvec=Amatmat, matmat=Amatmat)
>>> evals, evecs = primme.eigsh(A, 3, tol=1e-6, which='LA')
>>> evals
array([99., 98., 97.])