Shrinking the search subspace in Jacobi-Davidson

Once the residual norm of a Ritz pair is small enough, we say a Ritz pair has converged. Until now I have not mentioned how we can find the next eigenpair, but it is clear that the converged eigenvector should not be part of the search space anymore, so that we can find eigenvector corresponding to eigenvalues farther away from our target.

Another reason to shrink the search subspace is when it simply becomes too big. A large search space not only requires much memory, but is computationally expensive as well. We want to keep the computational costs of orthogonalization and extraction fixed, and therefore we simply impose a maximum dimension $m_{max}$.

Hermetian problems and shrinking the search space

Hermetian matrices have orthogonal eigenvectors for distinct eigenvalues. This is useful, because we work with an orthonormal basis for our search space as well. In fact this orthogonality property is most pronounced when we shrink the search space.

Suppose $A$ is Hermetian and the columns of the orthonormal matrix $V$ span the search subspace. Then the Galerkin projection of $A$ onto the search subspace gives rise to a Hermetian matrix $M = V^*AV$. So if its eigenvalues are distinct, its has an eigendecomposition $MS = S\Lambda$ where $S$ is unitary and $\Lambda = diag(\theta_1, \cdots, \theta_m)$.

Suppose that the Ritz pair $(\theta_1, Vs_1)$ is converged, then we can remove it from our search space by choosing a new basis $V$ in terms of all Ritz vectors excluding $Vs_1$. This is equivalent as updating $V$ as

$V \leftarrow VS[:, 2 : m].$

Similarly for restarts. If the search subspace becomes to large so that $m = m_{max}$, we decide to keep only the most promising $m_{min}$ Ritz vectors. That would for instance mean updating $V$ as

$V \leftarrow VS[:, 1 : m_{min}].$

Each time we update our basis of the search space $V$, we must also update our Galerkin approximation of $A$. Since our new basis consists of Ritz vectors, our Galerkin approximation will be diagonal:

$(VS)^*A(VS) = S^*MS = S^*S\Lambda = \Lambda.$

Note that we could also restart with our current best approximate eigenvector (that is $m_{min} = 1$), but usually it’s a shame to throw away so many good directions. IRAM is another eigensolver that restarts with a small basis of Ritz vectors, but this restart is more subtle, as one must also have that the shrunken search space is a Krylov subspace.

Change of basis without temporaries

The following is probably not a real issue these days, but if a restart is necessary because of memory limitations, you really would not want to store temporarily both the old basis $V$ together with the new basis $VS$. It is however possible to compute this matrix-matrix product in-place.

To do so (assuming $S$ is square), compute the LU decomposition $S = LU$, and compute the product

$V \leftarrow (V * L) * U$

column-wise from left to right with $L$ and from right to left with $U$. If we must compute the product with only a few columns of $S$, then we could perform only a partial LU decomposition. Lastly, for a stable LU decomposition we need pivoting, but I will not go into that detail.

Non-Hermetian problems

Non-Hermatian matrices do not necessarily have orthonormal eigenvectors, and since we really want to work with orthonormal vectors, an alternative is to use the Schur decomposition. More about this in a next post.