What are the advantages of the Davidson diagonalization method over other sparse matrix diagonalization methods?

Question

I am interested to understand the advantages of the Davidson diagonalization method over other sparse matrix diagonalization routines. For instance, Intel MKL.

Answer 1

The Davidson-Method is the best algorithm when you want a comparatively small number of eigenvectors/eigenvalues from a sparse, diagonally dominant matrix.

MKL doesn’t have a sparse matrix diagonalization method, but it’s a good idea to use it for the Linear Algebra operations used when implementing Davidson. If you want ALL the eigenvectors then you definitely want a dense eigensolver like what MKL offers.

Previous to the Davidson Algorithm the methods used to calculate a few eigenvectors were the Power method, and the Raleigh-Ritz Method. The power method is still widely used today in applications but IIRC the advantage of Davidson is that it's faster than these methods because If you have many roots you need to solve for, you don't need to update them all if some have converged, which saves a lot of time.

Answer 2

Following Nike Dattani's suggestion, I now add some comparisons of the Davidson method with methods that are more "modern" than it, to supplement the existing answers which compared the Davidson method with methods that predates it. However my answer will be more or less off-topic because I will mainly talk about disadvantages of the Davidson method compared to these methods. (It would sound strange if I introduce some new methods just to reveal their disadvantages!)

The Davidson method was developed in the 1970s (E. R. Davidson, J. Comput. Phys. 1975, 17, 87), and basically consists of:

1. Projecting the matrix to be diagonalized, $H$, onto a basis $\{X\}$ whose dimension is much smaller than $H$ itself.
2. Perform exact diagonalization on the projected matrix, giving eigenvalues $\{E\}$ and eigenvectors $\{V\}$. If the desired eigenvectors of $H$ can be exactly spanned by the basis $\{X\}$, then $\{E\}$ are exact eigenvalues of $H$, and exact eigenvectors of $H$ can be obtained from $\{V\}$: \begin{align} C_k = \sum_i X_i V_{ik} \tag{1} \end{align} However as $\{X\}$ has a much lower dimension than $H$ does, this condition is normally not exactly fulfilled, and we must expand $\{X\}$ until the desired eigenvectors are entirely contained in that space.
3. Compute approximate eigenvectors $C$ according to Eq. (1), and check if they have converged. The latter is done by verifying if $HC - EC$ is (approximately) the zero vector.
4. If the eigenvectors have not converged yet, expand the basis $\{X\}$ in the hope that we can expand the exact eigenvectors better with the enlarged basis. To expand the basis in the most efficient way, we resort to the residual vector $HC - EC$, which roughly tells us in which direction $C$ is wrong. We then selectively include some basis vectors that point to roughly the same direction as $HC - EC$ does (but not quite, since even better results are given by first scaling the residual vector by a matrix called the preconditioner - this is the essence of the Davidson method).
5. Go to step 1, using the enlarged basis instead.

(Note that the very fact that the Davidson method is faster than dense matrix diagonalization routines, such as DSYEV in LAPACK, is because we only want a small subset of the matrix's eigenvalues and eigenvectors. Otherwise, the Davidson method, at least in its original form, cannot be faster than DSYEV even if the matrix is sparse. Conversely, if we want only a few eigenvalues and eigenvectors, then Davidson is faster than DSYEV even if the matrix is dense.)

Despite the great success of the Davidson method, it has two drawbacks:

1. The memory cost increases during the iterations. For large scale calculations, the basis $\{X\}$ thus must be periodically reduced by discarding some of the basis vectors, in order to avoid memory bottlenecks, but this frequently throws important basis vectors away, leading to poor convergence.
2. It's hard, or even impossible, to converge some eigenvalues in the middle of the spectrum of $H$, without calculating either all eigenvalues that are smaller than them, or all eigenvalues that are larger than them. For example, suppose you want to calculate the carbon K-edge X-ray absorption spectrum of an organic molecule with TDDFT. With the Davidson method, you have to either (1) calculate all the eigenvalues from the visible region to UV to X-ray, until you reach the carbon K-edge; or (2) calculate all the eigenvalues from the highest energy X-ray absorption downward, until you reach the carbon K-edge. Neither is affordable.

One remedy of the first problem is to throw away certain linear combinations of the basis vectors in $\{X\}$ that are known to be relatively unimportant. If this is done at each iteration, rather than only when the memory bottleneck is reached, then it is found that the second problem can also be solved. To see this, suppose that the user wants to calculate all eigenvalues within an interval $[ E_1, E_2 ]$. At one point the expansion space $\{X\}$ may contain a basis vector that is a linear combination of two eigenvectors, the eigenvalue of one of them being smaller than $E_1$, the other larger than $E_2$. When diagonalizing the projected matrix, this basis vector will give a spurious eigenvalue that may fall within the window $[ E_1, E_2 ]$. If $\{X\}$ approaches a fixed size, then as the basis $\{X\}$ becomes better and better, such spurious eigenvectors will eventually disappear, as the basis will eventually be good enough for the algorithm to discover that the eigenvalue is spurious. However, if $\{X\}$ keeps enlarging, then more spurious eigenvalues keep appearing in the energy window while we remove the existing spurious eigenvalues, which easily makes the program resemble Sisyphus.

Based on this idea (remove certain linear combination of basis vectors at every iteration), several successors of the Davidson method emerged in recent decades, e.g. GPLHR, LOBPCG and iVI. Of these, iVI is the most general, and includes LOBPCG and a variant of GPLHR as special cases. The details of the algorithm (http://doi.org/10.1002/jcc.24907) are rather technical so I'll not reproduce them here. Instead I only mention some practical benefits of iVI:

1. iVI can be used to calculate core excitation spectra (X-ray absorption spectra), which the Davidson method struggles to converge;
2. When doing a SCF calculation of a material that has many energy bands, one can parallelize the calculation with respect to the energy bands, so that each computer node solves a few bands, and during the calculation there is no internode communication; this is in contrast to most programs where the calculation can only be parallelized on the level of atoms, basis functions and/or grid points, thereby involving a lot of communication overhead.
3. Likewise, a TDDFT calculation can be parallelized based on wavelength ranges. Suppose you have a large molecule that has hundreds of excited states in the visible region. You can compute the $> 600 \textrm{ nm}$ part of the spectrum on one node, the $500 \sim 600 \textrm{ nm}$ part on a second node, the $420 \sim 500 \textrm{ nm}$ part on a third one, and the $380 \sim 420 \textrm{ nm}$ part on a fourth one. Finally you combine the four spectra to yield the whole visible absorption spectrum.
4. You can compute the absorption spectrum within a wavelength range by directly specifying the wavelength range itself, and you do not need to guess the number of roots that you need to calculate - the program will automatically find this out.
5. iVI requires less memory than the Davidson method, typically by 2x or more.

Finally, back to topic. What are the advantages of the Davidson method compared to iVI etc.? One can name a few:

1. The Davidson method is easier to implement than iVI, and open-source implementations of the former are easier to find;
2. For TDDFT calculations, iVI sometimes suffer from the accumulation of numerical errors that affect the orthogonality of the basis vectors, and require additional matrix-vector product calculations to remove, although this problem is specific to full TDDFT, and is absent when the matrix to be diagonalized is Hermitian, for example in SCF or the Tamm-Dancoff approximation of TDDFT. Note that the Davidson method also suffer from similar orthogonality errors, but they are very cheap to remove.

Answer 3

To supplement Cody's excellent answer, Davidson diagonalization originates from quantum chemistry, where solving configuration interaction problems requires diagonalizing the Hamiltonian in the space of electron configurations. The Davidson algorithm and refinements thereof enabled full configuration interaction calculations with matrix sizes of $10^9 \times 10^9$ already in the early 1990s. The key in the algorithm is that the matrix is extremely sparse, since matrix elements are nonzero only in the case where two configurations differ by a double excitation, at most; moreover, the diagonal elements are energies of individual configurations which increase monotonically, and are larger than the off-diagonal values.

An alternative to Davidson, which appears to be popular among physicists, is the Lanczos algorithm. Davidson and Lanczos are somewhat similar algorithms, in that both operate on a vector subspace which is used to estimate the eigenvalues and eigenvectors of the full matrix.

However, according to an older colleague who worked on large configuration interaction calculations in the 1990s, the Davidson algorithm has much better and more reliable convergence properties, because the vectors added at each step lead to an optimal decrease of the energy; this is especially true for the reformulations due to e.g. Jeppe Olsen that employ better preconditioning and eliminate potential numerical instabilities that appear for gigantic matrices. At variance, the Lanczos method may stagnate for a long time (for tens to hundreds of iterations), and then change suddenly as the algorithm lands on a new, more important basis vector. This makes it hard to trust the results of Lanczos calculations, as seemingly converged values may turn out to be not converged. In contrast, you can always trust Davidson if it has been converged to a suitably small threshold.