Example 11.2.2.

Consider the model problem from finite difference discretization of Poisson's equation on the unit square. The main diagonal is constant (all 4's), so diagonal preconditioning is ineffective. It merely shifts and scales the spectrum, without changing the relative separation of the eigenvalues. However, incomplete factorization techniques [326] have been shown effective for preconditioning the corresponding linear equations, and some theoretical results have been proved. With no preconditioning, the conjugate gradient method converges in $O(n^{1/2})$ iterations. With incomplete factorization IC(0), convergence is faster, but still approximately $O(n^{1/2})$. However with MIC(0), modified incomplete factorization [212], convergence is $O(n^{1/4})$. It is shown in [330] that the same holds true for computing the smallest eigenvalue. The modified incomplete factorization is of $A$, not $A-\theta I$. The theoretical result can be backed up by computations. Table 11.1 gives numbers of iterations as the problem size increases. The Lanczos algorithm (corresponding to both no preconditioning and diagonal preconditioning) and IC(0) and MIC(0) preconditionings are used. The starting vector has random values from the interval (-1,1), and the convergence test is residual norm less than $10^{-8}$. No restarting is used. For large problems, MIC(0) preconditioning is much better. In a log-log plot, the slopes are roughly 0.56 for Lanczos and 0.43 for IC(0). So both are near to the expected 0.5, although IC(0) is significantly better. With MIC(0) preconditioning, the slope is about 0.23, which is near the value of 0.25 predicted by the $O(n^{1/4})$ convergence result.

For computing several eigenvalues, MIC is not as good. We look at computing the smallest five eigenvalues with $n=1600$. Again the preconditioners are from incomplete factorizations of $A$, with no attempt to adjust for the particular eigenvalue being computed. This time we restart the Davidson methods, with maximum subspaces of dimension 20, but do not restart Lanczos (if eigenvectors are desired, Lanczos might also need to be restarted). Relatively speaking, Lanczos does better at computing several eigenvalues, because the Davidson method targets one eigenvalue at a time. IC(0) uses 159 iterations, MIC(0) requires 158, and Lanczos 172. However, it should be noted that Lanczos misses the fact that one of the small eigenvalues is double, while the Davidson approaches do compute two eigenvalues at that point.

Table 11.1: Preconditioning the model problem

$n$	Lanczos	IC(0)	MIC(0)
25	13	11	12
49	25	14	15
100	36	17	17
196	51	24	20
400	69	27	22
784	99	40	26
1600	137	49	30
3136	182	64	34
6400	287	101	40

For a sparse matrix such as in Example 11.2.2, more than just the number of iterations must be considered. The Lanczos algorithm can be much cheaper per iteration than the Davidson approach. This motivates the development of the preconditioned Lanczos method which takes advantage of preconditioning while retaining much of the efficiency of the Lanczos approach. This will be discussed in §11.2.6. Efficiency issues also lead to connections with the Jacobi-Davidson method. One way to reduce the number of iterations required by the Davidson method and thus reduce the overhead for the Rayleigh-Ritz procedure is to develop an extremely good preconditioner. An iterative method can be used to approximately solve $(A-\theta B)w_j = r$ for the new vector $w_j$ in step (3) of Algorithm 11.2. And this solution can be made as accurate as is desired. In the symmetric case, the conjugate gradient method can be used for the iterative method, and it can be preconditioned itself. As with the preconditioned Lanczos method, we then have an inner loop that is efficient in terms of computations. But it is also efficient in terms of storage requirements. There is, however, danger in developing a preconditioner that is too accurate. Exact solution of $(A-\theta B)w_j = r$ gives $w_j=y$, the approximate eigenvector which is already in the subspace. This could be dealt with by instead solving $(A-\mu B)w = r$, with $\mu \ne \theta$, as in the inexact Cayley transform. However, in the Jacobi-Davidson method a different approach is taken; the approximate eigenvector is deflated from the solution of $(A-\theta B)w_j = r$. Discussion of connections with Jacobi-Davidson is continued in the next subsection.

Davidson's method can have difficulty when the starting vector or starting subspace is not very accurate. In [335] and [172], it is suggested that initially $M$ be set to $P - \mu I$, where $P$ is an approximation to $A$ and $\mu$ is near the desired eigenvalues. As mentioned in §11.1, another possibility is to use a Krylov subspace method in an initial phase to generate starting vectors for the Davidson method.