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

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

• intel MKL isn't a diagonalization method. A better comparison would be the Power method, which was used before Davidson became popular and is still used in a lot of applications (e.g. recommendation systems on the internet). Aug 13 at 16:17
• IIRC tne advantage of Davidson is that it's faster. 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. Aug 13 at 16:20
• @CodyAldaz Sounds like you have nearly enough to write an answer! If it's within the next 6 hours the user may get an HNQ (which may help attract some more opinions from people at the Computational Sciences SE, or the CS one, or Math one, or other places. Aug 13 at 17:56
• +1. This related question just happened to be at the top of the Computational Science SE page just now! scicomp.stackexchange.com/q/1681/27325 Aug 13 at 18:06
• Do you mean the Davidson method in the narrow sense, or do you mean Davidson-like methods in general? There are very recent sparse matrix diagonalization algorithms that on one hand cut down the memory cost significantly compared to Davidson, and on the other hand allows one to compute all eigenvalues within any energy interval, rather than only able to calculate the lowest ones. For example the iVI method: doi.org/10.1002/jcc.24907. Whether or not to include methods like iVI in your definition of the Davidson method will make the question quite different. Aug 13 at 19:35

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.

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.

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.