[6] Source:[7] The method performs an iterative maximization (or minimization) of the generalized Rayleigh quotient which results in finding largest (or smallest) eigenpairs of
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.
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
Adding more vectors, see, e.g., Richardson extrapolation, does not result in significant acceleration[8] but increases computation costs, so is not generally recommended.
become nearly linearly dependent, resulting in a precision loss and making the Rayleigh–Ritz method numerically unstable in the presence of round-off errors.
[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.
Extreme simplicity and high efficiency of the single-vector version of LOBPCG make it attractive for eigenvalue-related applications under severe 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.
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.
However, in inexact computer arithmetic the Rayleigh–Ritz method becomes numerically unstable if some of the basis vectors are approximately linearly dependent.
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.
For example, LOBPCG implementations,[9][10] utilize unstable but efficient Cholesky decomposition of the normal matrix, which is performed only on individual matrices
The global Rayleigh-Ritz procedure for all desired eigenpairs is only applied periodically at the end of a fixed number of unblocked LOBPCG iterations.
Individually running branches of the single-vector LOBPCG may not follow continuous iterative paths flipping instead and creating duplicated approximations to the same eigenvector.
LOBPCG by construction is guaranteed[8] to minimize the Rayleigh quotient not slower than the block steepest gradient descent, which has a comprehensive convergence theory.
The worst value of the linear convergence rate has been determined[8] and depends on the relative gap between the eigenvalue and the rest of the matrix spectrum and the quality of the preconditioner, if present.
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.
A good quality random Gaussian function with the zero mean is commonly the default in LOBPCG to generate the initial approximations.
The main calculation is evaluation of a function of the product DT(D X) of the covariance matrix DTD and the block-vector X 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.
[9] LOBPCG for PCA and SVD is implemented in SciPy since revision 1.4.0[13] LOBPCG's inventor, Andrew Knyazev, published a reference implementation called Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX)[14][15] with interfaces to PETSc, hypre, and Parallel Hierarchical Adaptive MultiLevel method (PHAML).
[16] Other implementations are available in, e.g., GNU Octave,[17] MATLAB (including for distributed or tiling arrays),[9] Java,[18] Anasazi (Trilinos),[19] SLEPc,[20][21] SciPy ,[10] Julia,[22] MAGMA,[23] Pytorch,[24] Rust,[25] OpenMP and OpenACC,[26] CuPy (A NumPy-compatible array library accelerated by CUDA),[27] Google JAX,[28] and NVIDIA AMGX.
Software packages scikit-learn and Megaman[31] use LOBPCG to scale spectral clustering[32] and manifold learning[33] via Laplacian eigenmaps to large data sets.
[37] It has been used for multi-billion size matrices by Gordon Bell Prize finalists, on the Earth Simulator supercomputer in Japan.
[41] There are MATLAB[42] and Julia[43][44] versions of LOBPCG for Kohn-Sham equations and density functional theory (DFT) using the plane wave basis.
[49] 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.
[51] 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.
[55] The latter approach has been later implemented in Python scikit-learn[56] that uses LOBPCG from SciPy with algebraic multigrid preconditioning for solving the eigenvalue problem for the graph Laplacian.