LOBPCG

Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) is a matrix-free method for finding the largest (or smallest) eigenvalues and the corresponding eigenvectors of a symmetric generalized eigenvalue problem
for a given pair formula_2 of complex Hermitian or real symmetric matrices, where
the matrix formula_3 is also assumed positive-definite.
Background.
Kantorovich in 1948 proposed calculating the smallest eigenvalue formula_4 of a symmetric matrix formula_5 by steepest descent using a direction formula_6 of a scaled gradient of a Rayleigh quotient formula_7 in a scalar product formula_8, with the step size computed by minimizing the Rayleigh quotient in the linear span of the vectors formula_9 and formula_10, i.e. in a locally optimal manner. Samokish proposed applying a preconditioner formula_11 to the residual vector formula_12 to generate the preconditioned direction formula_13 and derived asymptotic, as formula_9 approaches the eigenvector, convergence rate bounds. D'yakonov suggested spectrally equivalent preconditioning and derived non-asymptotic convergence rate bounds. Block locally optimal multi-step steepest descent for eigenvalue problems was described in. Local minimization of the Rayleigh quotient on the subspace spanned by the current approximation, the current residual and the previous approximation, as well as its block version, appeared in. The preconditioned version was analyzed in and.
Algorithm.
Single-vector version.
Preliminaries: Gradient descent for eigenvalue problems.
The method performs an iterative maximization (or minimization) of the generalized Rayleigh quotient
which results in finding largest (or smallest) eigenpairs of formula_16
The direction of the steepest ascent, which is the gradient, of the generalized Rayleigh quotient is positively proportional to the vector
called the eigenvector residual. If a preconditioner formula_11 is available, it is applied to the residual and gives the vector
called the preconditioned residual. Without preconditioning, we set 
formula_20 and so formula_21. An iterative method
or, in short,
is known as preconditioned steepest ascent (or descent), where the scalar 
formula_26 is called the step size. The optimal step size can be determined by maximizing the Rayleigh quotient, i.e.,
(or formula_28 in case of minimizing), 
in which case the method is called locally optimal.
Three-term recurrence.
To dramatically accelerate the convergence of the locally optimal preconditioned steepest ascent (or descent), one extra vector can be added to the two-term recurrence relation to make it three-term:
(use formula_28 in case of minimizing). The maximization/minimization of the Rayleigh quotient in a 3-dimensional subspace can be performed numerically by the Rayleigh–Ritz method. Adding more vectors, see, e.g., Richardson extrapolation, does not result in significant acceleration but increases computation costs, so is not generally recommended.
Numerical stability improvements.
As the iterations converge, the vectors formula_31 and formula_32 become nearly linearly dependent, resulting in a precision loss and making the Rayleigh–Ritz method numerically unstable in the presence of round-off errors. The loss of precision may be avoided by substituting the vector formula_32 with a vector formula_34, that may be further away from formula_35, in the basis of the three-dimensional subspace formula_36, while keeping the subspace unchanged and avoiding orthogonalization or any other extra operations. Furthermore, orthogonalizing the basis of the three-dimensional subspace may be needed for ill-conditioned eigenvalue problems to improve stability and attainable accuracy.
Krylov subspace analogs.
This is a single-vector version of the LOBPCG method—one of possible generalization of the preconditioned conjugate gradient linear solvers to the case of symmetric eigenvalue problems. Even in the trivial case formula_37 and formula_38 the resulting approximation with formula_39 will be different from that obtained by the Lanczos algorithm, although both approximations will belong to the same Krylov subspace.
Practical use scenarios.
Extreme simplicity and high efficiency of the single-vector version of LOBPCG make it attractive for eigenvalue-related applications under sever hardware limitations, ranging from spectral clustering based real-time anomaly detection via graph partitioning on embedded ASIC or FPGA to modelling physical phenomena of record computing complexity on exascale TOP500 supercomputers.
Block version.
Summary.
Subsequent eigenpairs can be computed one-by-one via single-vector LOBPCG supplemented with an orthogonal deflation or simultaneously as a block. In the former approach, imprecisions in already computed approximate eigenvectors additively affect the accuracy of the subsequently computed eigenvectors, thus increasing the error with every new computation. Iterating several approximate eigenvectors together in a block in a locally optimal fashion in the block version of the LOBPCG. allows fast, accurate, and robust computation of eigenvectors, including those corresponding to nearly-multiple eigenvalues where the single-vector LOBPCG suffers from slow convergence. The block size can be tuned to balance numerical stability vs. convergence speed vs. computer costs of orthogonalizations and the Rayleigh-Ritz method on every iteration.
Core design.
The block approach in LOBPCG replaces single-vectors formula_40 and formula_41 with block-vectors, i.e. matrices formula_42 and formula_43, where, e.g., every column of formula_44 approximates one of the eigenvectors. All columns are iterated simultaneously, and the next matrix of approximate eigenvectors formula_45 is determined by the Rayleigh–Ritz method on the subspace spanned by all columns of matrices formula_42 and formula_43. Each column of formula_48 is computed simply as the preconditioned residual for every column of formula_49 The matrix formula_43 is determined such that the subspaces spanned by the columns of formula_51 and of formula_52 are the same.
Numerical stability vs. efficiency.
The outcome of the Rayleigh–Ritz method is determined by the subspace spanned by all columns of matrices formula_42 and formula_43, where a basis of the subspace can theoretically be arbitrary. However, in inexact computer arithmetic the Rayleigh–Ritz method becomes numerically unstable if some of the basis vectors are approximately linearly depended. Numerical instabilities typically occur, e.g., if some of the eigenvectors in the iterative block already reach attainable accuracy for a given computer precision and are especially prominent in low precision, e.g., single precision.
The art of multiple different implementation of LOBPCG is to ensure numerical stability of the Rayleigh–Ritz method at minimal computing costs by choosing a good basis of the subspace. The arguably most stable approach of making the basis vectors orthogonal, e.g., by the Gram–Schmidt process, is also the most computational expensive. For example, LOBPCG implementations, utilize unstable but efficient Cholesky decomposition of the normal matrix, which is performed only on individual matrices formula_48 and formula_43, rather than on the whole subspace. The constantly increasing amount of computer memory allows typical block sizes nowadays in the formula_57 range, where the percentage of compute time spend on orthogonalizations and the Rayleigh-Ritz method starts dominating.
Locking of previously converged eigenvectors.
Block methods for eigenvalue problems that iterate subspaces commonly have some of the iterative eigenvectors converged faster than others that motivates locking the already converged eigenvectors, i.e. removing them from the iterative loop, in order to eliminate unnecessary computations and improve numerical stability. A simple removal of an eigenvector may likely result in forming its duplicate in still iterating vectors. The fact that the eigenvectors of symmetric eigenvalue problems are pair-wise orthogonal suggest keeping all iterative vectors orthogonal to the locked vectors.
Locking can be implemented differently maintaining numerical accuracy and stability while minimizing the compute costs. For example, LOBPCG implementations, follow, separating hard locking, i.e. a deflation by restriction, where the locked eigenvectors serve as a code input and do not change, from soft locking, where the locked vectors do not participate in the typically most expensive iterative step of computing the residuals, however, fully participate in the Rayleigh—Ritz method and thus are allowed to be changed by the Rayleigh—Ritz method.
Convergence theory and practice.
LOBPCG by construction is guaranteed to minimize the Rayleigh quotient not slower than the block steepest gradient descent, which has a comprehensive convergence theory. Every eigenvector is a stationary point of the Rayleigh quotient, where the gradient vanishes. Thus, the gradient descent may slow down in a vicinity of any eigenvector, however, it is guaranteed to either converge to the eigenvector with a linear convergence rate or, if this eigenvector is a saddle point, the iterative Rayleigh quotient is more likely to drop down below the corresponding eigenvalue and start converging linearly to the next eigenvalue below. The worst value of the linear linear convergence rate has been determined and depends on the relative gap between the eigenvalue and the rest of the matrix spectrum and the quality of the preconditioner, if present.
For a general matrix, there is evidently no way to predict the eigenvectors and thus generate the initial approximations that always work well. The iterative solution by LOBPCG may be sensitive to the initial eigenvectors approximations, e.g., taking longer to converge slowing down as passing intermediate eigenpairs. Moreover, in theory, one cannot guarantee convergence necessarily to the smallest eigenpair, although the probability of the miss is zero. A good quality random Gaussian function with the zero mean is commonly the default in LOBPCG to generate the initial approximations. To fix the initial approximations, one can select a fixed seed for the random number generator.
In contrast to the Lanczos method, LOBPCG rarely exhibits asymptotic superlinear convergence in practice.
Partial Principal component analysis (PCA) and Singular Value Decomposition (SVD).
LOBPCG can be trivially adopted for computing several largest singular values and the corresponding singular vectors (partial SVD), e.g., for iterative computation of PCA, for a data matrix with zero mean, without explicitly computing the covariance matrix , i.e. in matrix-free fashion. The main calculation is evaluation of a function of the product of the covariance matrix and the block-vector that iteratively approximates the desired singular vectors. PCA needs the largest eigenvalues of the covariance matrix, while LOBPCG is typically implemented to calculate the smallest ones. A simple work-around is to negate the function, substituting for and thus reversing the order of the eigenvalues, since LOBPCG does not care if the matrix of the eigenvalue problem is positive definite or not.
LOBPCG for PCA and SVD is implemented in SciPy since revision 1.4.0
General software implementations.
LOBPCG's inventor, Andrew Knyazev, published a reference implementation called Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) with interfaces to PETSc, hypre, and Parallel Hierarchical Adaptive MultiLevel method (PHAML). Other implementations are available in, e.g., GNU Octave, MATLAB (including for distributed or tiling arrays), Java, Anasazi (Trilinos), SLEPc, SciPy
, Julia, MAGMA, Pytorch, Rust, OpenMP and OpenACC, CuPy (A NumPy-compatible array library accelerated by CUDA), and NVIDIA AMGX. LOBPCG is implemented, but not included, in TensorFlow.
Applications.
Material sciences.
LOBPCG is implemented in ABINIT (including CUDA version) and Octopus. It has been used for multi-billion size matrices by Gordon Bell Prize finalists, on the Earth Simulator supercomputer in Japan. 
Hubbard model for strongly-correlated electron systems to understand the mechanism behind the superconductivity uses LOBPCG to calculate the ground state of the Hamiltonian on the K computer. 
There are MATLAB and Julia
versions of LOBPCG for [[Kohn-Sham]] equations and [[density functional theory]] (DFT) using the plain-wave basis.
Recent implementations include TTPY, Platypus‐QM, MFDn, ACE-Molecule, LACONIC.
[[Mechanics]] and [[fluid]]s.
LOBPCG from BLOPEX is used for [[preconditioner]] setup in Multilevel [[Balancing Domain Decomposition by Constraints]] (BDDC) solver library BDDCML, which is incorporated into OpenFTL (Open [[Finite element]] Template Library) and Flow123d simulator of underground water flow, solute and [[heat transport]] in fractured [[porous media]]. LOBPCG has been implemented in [[LS-DYNA]].
[[Maxwell's equations]].
LOBPCG is one of core eigenvalue solvers in PYFEMax and high performance multiphysics [[finite element]] software Netgen/NGSolve. LOBPCG from [[hypre]] is incorporated into [[open source]] lightweight scalable [[C++]] library for [[finite element]] methods [[MFEM]], which is used in many projects, including [[BLAST]], XBraid, [[VisIt]], xSDK, the FASTMath institute in [[SciDAC]], and the co-design Center for Efficient Exascale Discretizations (CEED) in the [[Exascale computing]] Project.
[[Denoising]].
Iterative LOBPCG-based approximate [[low-pass filter]] can be used for [[denoising]]; see, e.g., to accelerate [[total variation denoising]].
[[Image segmentation]].
[[Image segmentation]] via [[spectral clustering]] performs a low-dimension [[embedding]] using an [[Affine transformation|affinity]] matrix between pixels, followed by clustering of the components of the eigenvectors in the low dimensional space. LOBPCG with [[multigrid]] [[preconditioning]] has been first applied to [[image segmentation]] in via [[spectral]] [[graph partitioning]] using the [[graph Laplacian]] for the [[bilateral filter]]. 
[[Scikit-learn]] uses LOBPCG from [[SciPy]] with [[Multigrid method#Algebraic multigrid (AMG)|algebraic multigrid preconditioning]] for solving the eigenvalue problem.
[[Data mining]].
Software packages [[scikit-learn]] and Megaman use LOBPCG to scale [[spectral clustering]] and [[manifold learning]] via [[Nonlinear dimensionality reduction#Laplacian eigenmaps|Laplacian eigenmaps]] to large data sets. [[NVIDIA]] has implemented LOBPCG in its nvGRAPH library introduced in [[CUDA]] 8.
External links.
[[Category:Numerical linear algebra]]
[[Category:Scientific simulation software]]