Sparse Eigensolvers

Sparse eigensolvers compute selected eigenvalues and eigenvectors of large matrices without forming or diagonalizing a dense matrix. In quantum mechanics, the matrix is often a sparse Hamiltonian, and the desired output is usually the ground state, a few low-lying excited states, or a small set of eigenvalues near a chosen energy.

The key idea is to learn spectral information from repeated matrix-vector products,

$$v \mapsto Hv,$$

rather than from dense factorization of the whole matrix. This is why the storage and matrix-vector-product language of Sparse Matrices is prerequisite.

Problem Solved

Let $H$ be a large finite Hamiltonian matrix. A sparse eigensolver typically seeks a small number $k$ of eigenpairs

$$Hv_j = E_jv_j, \qquad j=1,\dots,k,$$

where $k\ll n$ and $n$ is the matrix dimension.

For closed-system bound-state calculations, $H$ is usually Hermitian. Then the eigenvalues are real, the eigenvectors can be chosen orthonormal, and the most common target is the lowest part of the spectrum:

$$E_0\le E_1\le E_2\le\cdots.$$

Dense Matrix Diagonalization is still the right tool for modest matrices when the full spectrum is needed. Sparse eigensolvers are for the regime where storing and diagonalizing a dense matrix is the wrong computational model.

Krylov Subspaces

Most standard sparse eigensolvers are Krylov methods. Starting from a nonzero vector $q$, define

$$\mathcal K_m(H,q) = \operatorname{span} \{q,Hq,H^2q,\dots,H^{m-1}q\}.$$

This subspace contains vectors of the form $p(H)q$, where $p$ is a polynomial of degree less than $m$. If $q$ has overlap with the desired eigenvectors, repeated application of $H$ builds a subspace that can approximate them.

The method then projects $H$ onto this smaller subspace, solves the small projected eigenvalue problem, and lifts the approximate eigenvectors back to the full space. The approximate eigenvalues are called Ritz values, and the approximate eigenvectors are called Ritz vectors. Their ordered upper-bound and interlacing properties are established in Upper Bounds and the Min–Max Principle.

The same Krylov-subspace idea can approximate $e^A v$ for time evolution; see Matrix Exponentials Numerically.

Projection Picture

Let the columns of $Q_m$ be an orthonormal basis for the Krylov subspace:

$$Q_m^\dagger Q_m = I_m.$$

The projected matrix is

$$H_m = Q_m^\dagger H Q_m.$$

Solving

$$H_m y = \theta y$$

gives a Ritz pair for the original problem:

$$\lambda\approx\theta, \qquad v\approx Q_m y.$$

The whole art is to build a useful $Q_m$ without losing orthogonality, wasting memory, or confusing convergence of Ritz values with convergence of eigenvectors.

Lanczos Method

For Hermitian matrices, the Lanczos method builds an orthonormal Krylov basis using a three-term recurrence. With $q_0=0$ and $\beta_1=0$, the recurrence has the schematic form

$$\beta_{j+1}q_{j+1} = Hq_j - \alpha_j q_j - \beta_j q_{j-1},$$

where

$$\alpha_j = q_j^\dagger Hq_j.$$

In exact arithmetic, the projected matrix is tridiagonal:

$$T_m = \begin{pmatrix} \alpha_1 & \beta_2 & 0 & \cdots & 0\\ \beta_2 & \alpha_2 & \beta_3 & \cdots & 0\\ 0 & \beta_3 & \alpha_3 & \ddots & \vdots\\ \vdots & \vdots & \ddots & \ddots & \beta_m\\ 0 & 0 & \cdots & \beta_m & \alpha_m \end{pmatrix}.$$

The eigenvalues of $T_m$ approximate selected eigenvalues of $H$. For ground-state computations, the lowest Ritz value often converges rapidly if the starting vector has nonzero overlap with the ground state.

Lanczos is attractive because each iteration needs one matrix-vector product, a few vector operations, and only short recurrence data in exact arithmetic. In finite precision, however, loss of orthogonality can produce repeated or spurious Ritz values. Practical implementations use reorthogonalization, restarts, locking, or thick-restart variants.

Lanczos Method Preview applies this machinery to symmetry-resolved many-body ground states, residual evidence, and continued-fraction response spectra.

Arnoldi Method

Arnoldi iteration is the corresponding Krylov method for general, possibly non-Hermitian matrices. It builds an orthonormal basis satisfying

$$HQ_m = Q_m A_m + h_{m+1,m}q_{m+1}e_m^T,$$

where $A_m$ is upper Hessenberg rather than tridiagonal.

Arnoldi is more expensive than Lanczos because orthogonalization against all previous basis vectors is usually required. It is the natural method when the operator is non-Hermitian, as in absorbing-boundary approximations, effective non-Hermitian Hamiltonians, Liouvillian problems, or some scattering discretizations.

For Hermitian Hamiltonians, Lanczos or a Hermitian-specialized variant is usually preferred.

Ground-State Computation

A common quantum workflow is:

1. Choose a basis, symmetry sector, grid, or finite Hilbert-space cutoff.
2. Implement the Hamiltonian as a sparse matrix or matrix-free operator.
3. Choose an initial vector with overlap with the desired state.
4. Run a Hermitian Krylov eigensolver for the lowest few Ritz pairs.
5. Check residuals, orthogonality, symmetry quantum numbers, and convergence under numerical refinement.

The initial vector matters. If it is exactly orthogonal to the ground state because of a symmetry, the algorithm will not find that ground state. This can be a feature when one intentionally works inside a symmetry sector, but it is a bug when the sector was chosen accidentally.

For near-degenerate low-energy states, computing one eigenpair at a time can be misleading. A block method or a request for several low-lying states is often safer, because the physically stable object may be the degenerate subspace rather than a particular numerical basis vector inside it.

Residual Diagnostics

For an approximate eigenpair $(\lambda,v)$, with $\lVert v\rVert=1$, the residual is

$$r = Hv-\lambda v.$$

The residual norm

$$\lVert r\rVert = \lVert Hv-\lambda v\rVert$$

is the most important diagnostic. A Ritz value that has stopped changing is not enough; the residual must be small at the scale required by the problem.

For a Hermitian $H$, if $\lambda=\langle v,Hv\rangle$, then

$$\lVert (H-\lambda I)v\rVert^2 = \langle v,H^2v\rangle - \langle v,Hv\rangle^2.$$

Thus the residual norm is the energy standard deviation of the trial state. This is a useful physics interpretation: an exact energy eigenstate has zero energy variance.

Convergence Checks

Report more than the final eigenvalue. Useful checks include:

• residual norms for each reported eigenpair;
• orthogonality of computed eigenvectors;
• stability of eigenvalues under tighter solver tolerances;
• stability under grid refinement, basis enlargement, or symmetry-sector checks;
• agreement with a small dense calculation when possible;
• convergence of expectation values, not only energies;
• absence of duplicate ghost eigenvalues caused by loss of orthogonality.

Solver tolerance is not the same as physical accuracy. A tiny Krylov residual only says that the finite matrix eigenproblem was solved accurately. It does not prove that the finite matrix accurately approximates the continuum or many-body problem.

Interior Eigenvalues and Shift-Invert

Extremal eigenvalues, such as the lowest energy, are usually easiest. Eigenvalues near an interior target $\sigma$ are harder. A standard transformation is shift-invert:

$$(H-\sigma I)^{-1}v = \mu v, \qquad \mu = \frac{1}{E-\sigma}.$$

Eigenvalues $E$ near $\sigma$ become large in magnitude after the transformation. A Krylov method can then target them as extremal eigenvalues of the transformed operator.

The price is that each iteration requires solving a linear system with $H-\sigma I$. That may require sparse factorization, preconditioning, or an inner iterative solve. Shift-invert can be powerful, but it changes the computational problem substantially.

Preconditioning and Davidson-Type Methods

Some sparse eigenvalue algorithms use approximate inverse information to accelerate convergence. Davidson and Jacobi-Davidson methods are important examples, especially when the Hamiltonian has a useful diagonal or block-diagonal approximation.

The basic idea is to correct a trial subspace using an approximate solution of a residual equation. These methods can be excellent in quantum chemistry and many-body calculations, but their reliability depends on the preconditioner and the spectral structure.

For a first pass, treat preconditioning as a controlled acceleration, not as a replacement for residual checks.

Degeneracy and Symmetry

Degeneracies require care. If several eigenvalues are equal or nearly equal, small numerical perturbations can rotate the computed eigenvectors inside the nearly degenerate subspace. The energies may be stable while individual eigenvectors are not.

Use symmetry quantum numbers, projectors, or block diagonalization when possible. If the physics depends on the subspace, check subspace convergence rather than the component-by-component agreement of individual eigenvectors.

This is the same mathematical warning that appears in dense diagonalization, but sparse iterative methods make it more visible because convergence can occur at different rates for different vectors in the cluster.

Common Mistakes

• Asking a sparse eigensolver for too many eigenpairs and expecting dense-diagonalization information cheaply.
• Reporting Ritz values without residual norms.
• Treating solver convergence as continuum convergence.
• Using an initial vector with no overlap with the desired symmetry sector.
• Ignoring loss of orthogonality in Lanczos iteration.
• Mistaking ghost eigenvalues for physical degeneracies.
• Computing only one state when a nearly degenerate multiplet is physically relevant.
• Using shift-invert without accounting for the cost and conditioning of the shifted linear solves.
• Comparing eigenvectors across runs without fixing phases, degeneracy conventions, and symmetry sectors.

Cross-Links

• Sparse Matrices
• Matrix Exponentials Numerically
• Matrix Diagonalization
• Conditioning and Stability
• Floating-Point Arithmetic
• Discretization
• Finite Difference Methods
• Eigenvalues and Eigenvectors
• Hermitian Operators
• Orthonormal Bases
• Spectral Decomposition
• Imaginary-Time Projection Notebook

References

• Y. Saad, Numerical Methods for Large Eigenvalue Problems, revised ed., SIAM, 2011.
• R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM, 1998.
• B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM, 1998.
• J. K. Cullum and R. A. Willoughby, Lanczos Algorithms for Large Symmetric Eigenvalue Computations, SIAM, 2002.
• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, 1997.

Exercises

1. Explain why every vector in $\mathcal K_m(H,q)$ can be written as $p(H)q$ for some polynomial $p$ of degree less than $m$.

Solution

By definition, a vector in the Krylov subspace is a linear combination

$$c_0q+c_1Hq+\cdots+c_{m-1}H^{m-1}q.$$

This equals $p(H)q$ with

$$p(z) = c_0+c_1z+\cdots+c_{m-1}z^{m-1}.$$

2. Let $v$ be normalized and set $\lambda=\langle v,Hv\rangle$ for Hermitian $H$. Show that the residual norm squared equals the energy variance.

Solution

Compute

$$\lVert (H-\lambda I)v\rVert^2 = \langle (H-\lambda I)v,(H-\lambda I)v\rangle.$$

Using Hermiticity and real $\lambda$,

$$\lVert (H-\lambda I)v\rVert^2 = \langle v,H^2v\rangle - 2\lambda\langle v,Hv\rangle + \lambda^2.$$

Since $\lambda=\langle v,Hv\rangle$ and $\lVert v\rVert=1$, this becomes

$$\langle v,H^2v\rangle - \langle v,Hv\rangle^2.$$

3. Why is Lanczos cheaper than Arnoldi for Hermitian matrices in exact arithmetic?

Solution

Hermiticity makes the projected Krylov matrix tridiagonal, so the new Lanczos vector only needs to be orthogonalized against the previous two basis directions in exact arithmetic. Arnoldi for a general matrix produces an upper-Hessenberg projected matrix and usually requires orthogonalization against all previous basis vectors.

4. In shift-invert iteration, why do eigenvalues near $\sigma$ become easier to target as extremal eigenvalues?

Solution

If $Hv=Ev$, then

$$(H-\sigma I)^{-1}v = \frac{1}{E-\sigma}v.$$

When $E$ is close to $\sigma$, the transformed eigenvalue $\mu=1/(E-\sigma)$ has large magnitude. Krylov methods that target large-magnitude extremal eigenvalues can therefore find eigenvalues of $H$ near $\sigma$, provided the shifted linear solves are accurate and well controlled.