LOBPCG finds the lowest states of a particle in a box

The problem

A particle in a square box of side 1, with $\hslash =m=1$, has the Hamiltonian $H=-\frac{1}{2}{\nabla }^2$ and $\psi =0$ on the walls. Put it on a grid of 24 by 24 interior points, spacing $h=1/25$, and the Laplacian becomes the five-point stencil: at each point, four times the value there minus the four neighbors, over $h^2$. The grid Hamiltonian is then a symmetric matrix with $576$ rows. Its eigenvalues are known exactly, so every number the solver produces can be checked. The mode with $j$ half-waves across and $k$ down has energy

a little below the continuous box's ${\pi }^2(j^2+k^2)/2$. We want the lowest six: $(1,1)$, then $(1,2)$ and $(2,1)$, which have the same energy, then $(2,2)$, then $(1,3)$ and $(3,1)$, again a pair.

Why not diagonalize

eig_sym on the 576 by 576 matrix would return all 576 eigenpairs, and its cost grows as the cube of the size. We want six of them. Worse, the matrix has to exist: $576^2$ numbers, of which at most five per row are not zero. On a 3-D grid of 100 points a side it would have $10^{12}$ entries, 8 TB of doubles, while applying it to a vector touches seven numbers per point.

How do you find the lowest eigenvectors of $H$ when the one thing you can do with it is compute $Hx$? In this program that is apply_box, a few whole-array shifts. In the plane-wave DFT code on this site it is two FFTs and a projector sum, and the difference is measured: silicon with the matrix built and diagonalized at every k-point takes several times as long as with this solver (the timings, native on one core, are in dft/README.md).

A minimum instead of a factorization

The Rayleigh quotient of a vector is $\rho (x)=x^{H}Hx/x^{H}x$. Write $x$ in the eigenvectors of $H$, $x={\sum }_i{c}_i{v}_i$, and it becomes ${\sum }_i|c_i{|}^2{\lambda }_i/{\sum }_i|c_i{|}^2$: a weighted average of the eigenvalues. An average is smallest when all the weight sits on the smallest value, so the minimum of $\rho$ over all vectors is ${\lambda }_1$, reached at the ground state. Finding the lowest eigenvector is therefore a minimization.

The gradient of $\rho$ is $2(Hx-\rho x)/x^{H}x$. The vector $r=Hx-\rho x$ is the residual: it is zero exactly when $x$ is an eigenvector, and its norm bounds how far $\rho$ is from the nearest eigenvalue. So the direction downhill costs one application of $H$, and the residual's norm tells you when to stop.

For the lowest several states, carry several vectors at once, kept orthonormal. The lowest $k$ eigenvalues are what the $k$ Rayleigh quotients add up to at the minimum, and the $k$ vectors that reach it span the lowest $k$ eigenvectors.

Three blocks

Steepest descent would move each vector along its residual, $x-\alpha r$, with $\alpha$ chosen to make $\rho$ as small as possible. Choosing that $\alpha$ is itself a small eigenvalue problem. For any set of orthonormal vectors $S$ (as rows), the best combinations of them are found by forming the small matrix $S^{\ast }H{S}^T$ and diagonalizing it; its eigenvalues are the lowest Rayleigh quotients available in the span of $S$ (Ritz values) and its eigenvectors say which combinations reach them (Ritz vectors). This is called Rayleigh-Ritz. With $S=[x;r]$ the best $\alpha$ comes out of a 2 by 2.

LOBPCG searches a bigger span. Each step it takes the lowest Ritz vectors of three blocks stacked together:

• x, the current approximations,
• w, their residuals, made orthogonal to x and scaled to unit length, and
• p, the step each approximation took the last time.

The third block is the conjugate-gradient idea. Steepest descent on a long, narrow valley zigzags from wall to wall, because each new step undoes part of the last one. Keeping the last step in the search space lets the method combine the two, and the combination points along the valley. "Locally optimal" means that instead of a formula for the step it takes the best point in the span of all three.

Here each block has eight rows, six wanted and two spare, so the eigenproblems a step solves are at most 24 by 24, cheap next to anything done to the full vectors. The spares are there because the highest vectors in a block converge slowest: nothing above them in the block keeps them separated from the rest of the spectrum. Degenerate pairs need no special treatment: any combination of $(1,2)$ and $(2,1)$ is an eigenvector with the same energy, and the block converges to some orthonormal pair spanning that plane.

The P in the name is for preconditioner: replace $r$ by an approximation to $(H-\lambda {)}^{-1}r$, which points closer to the correction needed. The DFT engine uses one built from the kinetic energy. This run uses none, so the figure is about the three blocks alone.

Locking

A state whose residual norm is under the tolerance ($10^{-8}$ here) is locked: its residual is left out of w, so no direction is spent pushing it further. It stays in the Rayleigh-Ritz step, though, so it still improves from the other states' directions. By then its residual is mostly rounding error, and adding it to the search space only makes the space more nearly dependent, the usual way LOBPCG runs into trouble (the notes in lib/lobpcg have the details).

What the run shows

The ground state is the first under $10^{-8}$, at step 75, and it keeps going down while locked, to $10^{-9.2}$ by the end. The two highest wanted states, $(1,3)$ and $(3,1)$, are the last, at steps 102 and 103, and the run stops there. Every level then agrees with the formula above to $3\times 10^{-13}$ (LEVEL_ERROR).

The eigenvalues get good much sooner than the vectors do. At step 50 the worst residual norm is $10^{-2.2}$ and the worst eigenvalue error is $10^{-6.9}$. That is because $\rho$ is flat at an eigenvector: an error $ϵ$ in the vector moves the Rayleigh quotient by about ${ϵ}^2$.

Run the same iteration without p, which converge does by telling every step it is the first, and it is block steepest descent. From the same starting vectors it needs 797 steps to meet the same tolerance.

The levels also sit where they should relative to the continuous box. The five-point stencil drops the $h^2$ term of the second derivative's Taylor series, which lowers each level by about $h^2{\pi }^4(j^4+k^4)/24$; the measured shifts match that to within half a percent.