LOBPCG

The method performs an iterative maximization (or minimization) of the generalized Rayleigh quotient

${\displaystyle \rho (x):=\rho (A,B;x):={\frac {x^{T}Ax}{x^{T}Bx}},}$

which results in finding largest (or smallest) eigenpairs of ${\displaystyle Ax=\lambda Bx.}$

The direction of the steepest ascent, which is the gradient, of the generalized Rayleigh quotient is positively proportional to the vector

${\displaystyle r:=Ax-\rho (x)Bx,}$

called the eigenvector residual. If a preconditioner ${\displaystyle T}$ is available, it is applied to the residual and gives the vector

${\displaystyle w:=Tr,}$

called the preconditioned residual. Without preconditioning, we set ${\displaystyle T:=I}$ and so ${\displaystyle w:=r}$. An iterative method

${\displaystyle x^{i+1}:=x^{i}+\alpha ^{i}T(Ax^{i}-\rho (x^{i})Bx^{i}),}$

or, in short,

${\displaystyle x^{i+1}:=x^{i}+\alpha ^{i}w^{i},\,}$

${\displaystyle w^{i}:=Tr^{i},\,}$

${\displaystyle r^{i}:=Ax^{i}-\rho (x^{i})Bx^{i},}$

is known as preconditioned steepest ascent (or descent), where the scalar ${\displaystyle \alpha ^{i}}$ is called the step size. The optimal step size can be determined by maximizing the Rayleigh quotient, i.e.,

${\displaystyle x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i}\}}\rho (y)}$

(or ${\displaystyle \arg \min }$ in case of minimizing), in which case the method is called locally optimal.

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:

${\displaystyle x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i},x^{i-1}\}}\rho (y)}$

(use ${\displaystyle \arg \min }$ 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[8] but increases computation costs, so is not generally recommended.

As the iterations converge, the vectors ${\displaystyle x^{i}}$ and ${\displaystyle x^{i-1}}$ 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 ${\displaystyle x^{i-1}}$ with a vector ${\displaystyle p^{i}}$, that may be further away from ${\displaystyle x^{i}}$, in the basis of the three-dimensional subspace ${\displaystyle span\{x^{i},w^{i},x^{i-1}\}}$, while keeping the subspace unchanged and avoiding orthogonalization or any other extra operations.[8] Furthermore, orthogonalizing the basis of the three-dimensional subspace may be needed for ill-conditioned eigenvalue problems to improve stability and attainable accuracy.

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.[8] Even in the trivial case ${\displaystyle T=I}$ and ${\displaystyle B=I}$ the resulting approximation with ${\displaystyle i>3}$ will be different from that obtained by the Lanczos algorithm, although both approximations will belong to the same Krylov subspace.

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.

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.[8] 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.

The block approach in LOBPCG replaces single-vectors ${\displaystyle x^{i},\,w^{i},}$ and ${\displaystyle p^{i}}$ with block-vectors, i.e. matrices ${\displaystyle X^{i},\,W^{i},}$ and ${\displaystyle P^{i}}$, where, e.g., every column of ${\displaystyle X^{i}}$ approximates one of the eigenvectors. All columns are iterated simultaneously, and the next matrix of approximate eigenvectors ${\displaystyle X^{i+1}}$ is determined by the Rayleigh–Ritz method on the subspace spanned by all columns of matrices ${\displaystyle X^{i},\,W^{i},}$ and ${\displaystyle P^{i}}$. Each column of ${\displaystyle W^{i}}$ is computed simply as the preconditioned residual for every column of ${\displaystyle X^{i}.}$ The matrix ${\displaystyle P^{i}}$ is determined such that the subspaces spanned by the columns of ${\displaystyle [X^{i},\,P^{i}]}$ and of ${\displaystyle [X^{i},\,X^{i-1}]}$ are the same.

The outcome of the Rayleigh–Ritz method is determined by the subspace spanned by all columns of matrices ${\displaystyle X^{i},\,W^{i},}$ and ${\displaystyle P^{i}}$, 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,[9][10] utilize unstable but efficient Cholesky decomposition of the normal matrix, which is performed only on individual matrices ${\displaystyle W^{i}}$ and ${\displaystyle P^{i}}$, rather than on the whole subspace. The constantly increasing amount of computer memory allows typical block sizes nowadays in the ${\displaystyle 10^{3}-10^{4}}$ range, where the percentage of compute time spend on orthogonalizations and the Rayleigh-Ritz method starts dominating.

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,[9][10] follow,[8][11] 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.

LOBPCG is implemented in ABINIT[29] (including CUDA version) and Octopus.[30] It has been used for multi-billion size matrices by Gordon Bell Prize finalists, on the Earth Simulator supercomputer in Japan.[31][32] 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.[33] There are MATLAB [34] and Julia[35][36] [37] versions of LOBPCG for Kohn-Sham equations and density functional theory (DFT) using the plain-wave basis. Recent implementations include TTPY,[38] Platypus‐QM,[39] MFDn,[40] ACE-Molecule,[41] LACONIC.[42]

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[43] in LS-DYNA.

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.

Iterative LOBPCG-based approximate low-pass filter can be used for denoising; see,[44] e.g., to accelerate total variation denoising.

Image segmentation via spectral clustering performs a low-dimension embedding using an 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 [45] via spectral graph partitioning using the graph Laplacian for the bilateral filter. Scikit-learn uses LOBPCG from SciPy with algebraic multigrid preconditioning for solving the eigenvalue problem.[46]

Software packages scikit-learn and Megaman[47] use LOBPCG to scale spectral clustering[48] and manifold learning[49] via Laplacian eigenmaps to large data sets. NVIDIA has implemented[50] LOBPCG in its nvGRAPH library introduced in CUDA 8.