scipy.sparse.linalg.lobpcg#
- scipy.sparse.linalg.lobpcg(A, X, B=None, M=None, Y=None, tol=None, maxiter=None, largest=True, verbosityLevel=0, retLambdaHistory=False, retResidualNormsHistory=False)[source]#
Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG).
LOBPCG is a preconditioned eigensolver for large symmetric positive definite (SPD) generalized eigenproblems.
- Parameters:
- A{sparse matrix, dense matrix, LinearOperator}
The symmetric linear operator of the problem, usually a sparse matrix. Often called the “stiffness matrix”.
- Xndarray, float32 or float64
Initial approximation to the
k
eigenvectors (non-sparse). If A hasshape=(n,n)
then X should have shapeshape=(n,k)
.- B{dense matrix, sparse matrix, LinearOperator}, optional
The right hand side operator in a generalized eigenproblem. By default,
B = Identity
. Often called the “mass matrix”.- M{dense matrix, sparse matrix, LinearOperator}, optional
Preconditioner to A; by default
M = Identity
. M should approximate the inverse of A.- Yndarray, float32 or float64, optional
An n-by-sizeY matrix of constraints (non-sparse), sizeY < n. The iterations will be performed in the B-orthogonal complement of the column-space of Y. Y must be full rank.
- tolscalar, optional
Solver tolerance (stopping criterion). The default is
tol=n*sqrt(eps)
.- maxiterint, optional
Maximum number of iterations. The default is
maxiter = 20
.- largestbool, optional
When True, solve for the largest eigenvalues, otherwise the smallest.
- verbosityLevelint, optional
Controls solver output. The default is
verbosityLevel=0
.- retLambdaHistorybool, optional
Whether to return eigenvalue history. Default is False.
- retResidualNormsHistorybool, optional
Whether to return history of residual norms. Default is False.
- Returns:
- wndarray
Array of
k
eigenvalues.- vndarray
An array of
k
eigenvectors. v has the same shape as X.- lambdaslist of ndarray, optional
The eigenvalue history, if retLambdaHistory is True.
- rnormslist of ndarray, optional
The history of residual norms, if retResidualNormsHistory is True.
Notes
If both
retLambdaHistory
andretResidualNormsHistory
are True, the return tuple has the following format(lambda, V, lambda history, residual norms history)
.In the following
n
denotes the matrix size andk
the number of required eigenvalues (smallest or largest).The LOBPCG code internally solves eigenproblems of the size
3k
on every iteration by calling the “standard” dense eigensolver, so ifk
is not small enough compared ton
, it does not make sense to call the LOBPCG code, but rather one should use the “standard” eigensolver, e.g. numpy or scipy function in this case. If one calls the LOBPCG algorithm for5k > n
, it will most likely break internally, so the code tries to call the standard function instead.It is not that
n
should be large for the LOBPCG to work, but rather the ration / k
should be large. It you call LOBPCG withk=1
andn=10
, it works thoughn
is small. The method is intended for extremely largen / k
.The convergence speed depends basically on two factors:
Relative separation of the seeking eigenvalues from the rest of the eigenvalues. One can vary
k
to improve the absolute separation and use proper preconditioning to shrink the spectral spread. For example, a rod vibration test problem (under tests directory) is ill-conditioned for largen
, so convergence will be slow, unless efficient preconditioning is used. For this specific problem, a good simple preconditioner function would be a linear solve for A, which is easy to code since A is tridiagonal.Quality of the initial approximations X to the seeking eigenvectors. Randomly distributed around the origin vectors work well if no better choice is known.
References
[1]A. V. Knyazev (2001), Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific Computing 23, no. 2, pp. 517-541. DOI:10.1137/S1064827500366124
[2]A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov (2007), Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc. arXiv:0705.2626
[3]A. V. Knyazev’s C and MATLAB implementations: https://github.com/lobpcg/blopex
Examples
Solve
A x = lambda x
with constraints and preconditioning.>>> import numpy as np >>> from scipy.sparse import spdiags, issparse >>> from scipy.sparse.linalg import lobpcg, LinearOperator
The square matrix size:
>>> n = 100 >>> vals = np.arange(1, n + 1)
The first mandatory input parameter, in this test a sparse 2D array representing the square matrix of the eigenvalue problem to solve:
>>> A = spdiags(vals, 0, n, n) >>> A.toarray() array([[ 1, 0, 0, ..., 0, 0, 0], [ 0, 2, 0, ..., 0, 0, 0], [ 0, 0, 3, ..., 0, 0, 0], ..., [ 0, 0, 0, ..., 98, 0, 0], [ 0, 0, 0, ..., 0, 99, 0], [ 0, 0, 0, ..., 0, 0, 100]])
Initial guess for eigenvectors, should have linearly independent columns. The second mandatory input parameter, a 2D array with the row dimension determining the number of requested eigenvalues. If no initial approximations available, randomly oriented vectors commonly work best, e.g., with components normally disrtibuted around zero or uniformly distributed on the interval [-1 1].
>>> rng = np.random.default_rng() >>> X = rng.normal(size=(n, 3))
Constraints - an optional input parameter is a 2D array comprising of column vectors that the eigenvectors must be orthogonal to:
>>> Y = np.eye(n, 3)
Preconditioner in the inverse of A in this example:
>>> invA = spdiags([1./vals], 0, n, n)
The preconditiner must be defined by a function:
>>> def precond( x ): ... return invA @ x
The argument x of the preconditioner function is a matrix inside
lobpcg
, thus the use of matrix-matrix product@
.The preconditioner function is passed to lobpcg as a
LinearOperator
:>>> M = LinearOperator(matvec=precond, matmat=precond, ... shape=(n, n), dtype=np.float64)
Let us now solve the eigenvalue problem for the matrix A:
>>> eigenvalues, _ = lobpcg(A, X, Y=Y, M=M, largest=False) >>> eigenvalues array([4., 5., 6.])
Note that the vectors passed in Y are the eigenvectors of the 3 smallest eigenvalues. The results returned are orthogonal to those.