Ritz pairs turn out to be good approximate eigenpairs when one is interested in exterior eigenpairs. If we are interested in interior eigenpairs (somewhere inside the convex hull of the spectrum), we must ensure that we do not restart Jacobi-Davidson with a set of Ritz vectors that correspond to Ritz values converging to exterior eigenvalues. The simple criterion of Ritz values being close to our target, might give us a low-quality search space at restart. Maybe we could make a better judgement about which Ritz vectors to keep if we would inspect their residuals, but that strategy would be far too expensive.
Shift-and-invert is a well-known trick to solve the above problem: rather than solving the eigenvalue problem for close to , we reformulate the problem to solving with . The eigenvectors remain unchanged, but the eigenvalues close to our target are transformed to exterior eigenvalues. The message is that we do not change our methods, rather we change the problem itself.
However, changing the problem like this can be very costly: iterating with a matrix requires solving a system each iteration. In fact, the Arnoldi method requires accurate solves because otherwise the underlying Krylov subspace would be inexact.
Harmonic Ritz values
The main question is whether we can get higher quality approximations to eigenvalues (without doing much work), so that we can restart with a search space that does not contain partly converged exterior eigenvectors.
One way to go is to look at the Ritz values of with respect to a specially chosen subspace , so that the term drops out:
is equivalent to the requirement
The latter is a Petrov-Galerkin projection: the search space is and the test space is . We refer to a pair as a harmonic Ritz pair.
The important thing is that is identical in both expressions. The second shows a practical means to compute it, while the first shows its interpretation as being a Ritz value of .
Now the line of thought is as follows: if the Ritz values converge from within the convex hull of the spectrum of towards the exterior eigenvalues, then conversely its reciprocal becomes small only when it approximates an eigenvalue of well. This in turn means that can only be close to an eigenvalue of near if it is nearly converged. And that is what we wanted to ensure at restart.
All in all we work with the practical Petrov-Galerkin projection. Let’s write so that the last equation can be rewritten as the generalized eigenvalue problem
For stability reasons we like to work with orthonormal matrices. The Gram-Schmidt procedure allows us to perform the QR-decomposition of into with orthonormal and upper triangular. If we also define similarly to what we had before, then the generalized eigenvalue problem becomes
assuming that is invertible (which is the case a long as the columns of are independent).
So we can just obtain harmonic Ritz pairs by computing eigenpairs from the small generalized eigenvalue problem, and transform them back to approximate eigenpairs by setting and .
In the implementation we will explicitly build these two small matrices and .
Schur vectors vs eigenvectors
We have used a Schur decomposition in favour of an eigendecomposition of the low-dimensional problem to make it easier to shrink the search subspace. Let’s see how this works with harmonic Ritz values.
We compute the generalized Schur (QZ) decomposition and for unitary matrices and upper triangular matrices . The eigenvalues are now and can be reordered. In Julia it is the familiar function
ordschur!. We must order them from small to large in absolute magnitude. It is not hard to verify that the first column of is an eigenvector with as eigenvalue.
Now if the harmonic Ritz values we want to keep are moved to the front, we can easily restart by updating both our search space as and our test space as .
Rayleigh quotient vs harmonic Ritz value
Harmonic Ritz values are an excellent choice to determine which Ritz vectors to use at restart. They are however suboptimal for expansion. A few posts ago we wanted the residual to be perpendicular to the approximate eigenvector , so that when solving the correction equation with a Krylov method, would automatically be perpendicular to .
With harmonic Ritz values this might not be the case. Therefore we can choose to use harmonic Ritz values when it comes to shrinking the search space, yet use the Rayleigh quotients when it comes to expanding it.
We obtain it by setting ; in our case and . This gives after some manipulations that , where the bar denotes complex conjugation.