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 positive definite generalized eigenvalue problem

$Ax=\lambda Bx,$

for a given pair $(A,B)$ of complex Hermitian or real symmetric matrices, where the matrix $B$ is also assumed positive-definite.

Background

Kantorovich in 1948 proposed calculating the smallest eigenvalue $\lambda _{1}$ of a symmetric matrix $A$ by steepest descent using a direction $r=Ax-\lambda (x)$ of a scaled gradient of a Rayleigh quotient $\lambda (x)=(x,Ax)/(x,x)$ in a scalar product $(x,y)=x'y$, with the step size computed by minimizing the Rayleigh quotient in the linear span of the vectors $x$ and $w$, i.e. in a locally optimal manner. Samokish proposed applying a preconditioner $T$ to the residual vector $r$ to generate the preconditioned direction $w=Tr$ and derived asymptotic, as $x$ 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.

Main features

• LOBPCG is matrix-free, i.e. does not require storing the coefficient matrix explicitly, but can accesses the matrix by evaluating matrix-vector products.
• The costs per iteration and the memory use in LOBPCG are competitive with those of the Lanczos method, computing a single extreme eigenpair.
• Linear convergence is theoretically guaranteed and practically observed, since local optimality implies that LOBPCG converges at least as fast as the gradient descent method. In numerical tests, LOBPCG typically shows no super-linear convergence.
• LOBPCG blocking allows utilizing highly efficient matrix-matrix operations, e.g., BLAS 3.
• LOBPCG can directly take advantage of preconditioning, in contrast to the Lanczos method. LOBPCG allows variable and non-symmetric as well as fixed and positive definite preconditioning.
• LOBPCG allows warm starts and computes an approximation to the eigenvector on every iteration. It has no numerical stability issues similar to those of the Lanczos method.
• LOBPCG is reasonably easy to implement, so many implementations have appeared.
• LOBPCG general technology can also be viewed as a particular case of generalized block Davidson diagonalization methods with thick restart, or accelerated block gradient descent with plane-search.
• Very large block sizes in LOBPCG become expensive to deal with due to orthogonalizations and the use of the Rayleigh-Ritz method on every iteration.

Algorithm

Single-vector version

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

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

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

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

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

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

$w:=Tr,$

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

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

or, in short,

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

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

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

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

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

(or $\arg \min$ in case of minimizing), in which case the method is called locally optimal. To further accelerate the convergence of the locally optimal preconditioned steepest ascent (or descent), one can add one extra vector to the two-term recurrence relation to make it three-term:

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

(use $\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. As the iterations converge, the vectors $x^{i}$ and $x^{i-1}$ become nearly linearly dependent, making the Rayleigh–Ritz method numerically unstable in the presence of round-off errors. It is possible to substitute the vector $x^{i-1}$ with an explicitly computed difference $p^{i}=x^{i-1}-x^{i}$ making the Rayleigh–Ritz method more stable; see.

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 $T=I$ and $B=I$ the resulting approximation with $i>3$ will be different from that obtained by the Lanczos algorithm, although both approximations will belong to the same Krylov subspace.

Block version

Iterating several approximate eigenvectors together in a block in a similar locally optimal fashion, gives the full block version of the LOBPCG. It allows robust computation of eigenvectors corresponding to nearly-multiple eigenvalues.

General software implementations

LOBPCG's inventor, Andrew Knyazev, published an implementation called Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) with interfaces to PETSc and hypre. Other implementations are available in, e.g., Octave, MATLAB, Java, Anasazi (Trilinos), SLEPc, SciPy, MAGMA and NVIDIA AMGX.

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. Recent implementations include TTPY, Platypus‐QM, and MFDn. 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.

Maxwell's equations

LOBPCG is one of core eigenvalue solvers in PYFEMax, NGSolve, and MFEM.

Denoising

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

Image Segmentation

Hypre implementation of LOBPCG with multigrid preconditioning has been applied to image segmentation in via spectral graph partitioning using the graph Laplacian for the bilateral filter.

Data Mining

Software package Megaman uses LOBPCG to scale manifold learning algorithms to large data sets. NVIDIA has implemented LOBPCG in its nvGRAPH library introduced in CUDA 8.

External links

• LOBPCG code in MATLAB
• LOBPCG code in Octave
• LOBPCG code in SciPy
• LOBPCG code in Java at Google Code
• LOBPCG in Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) at Bitbucket
• LOBPCG in Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) at Google Code (retiring)