Deflation

We now consider the following implementation which incorporates a deflation process. So far we have described algorithms that compute only one eigenpair. In case several eigenpairs are sought, there are two possible options.

The first is to take $v_1$ to be a linear combination of the approximate eigenvectors when we restart. For example, if we need to compute the $p$ rightmost eigenvectors, we may take

$$\hat v_1 = \sum_{i=1}^p \rho_i \tilde u_i,$$

where the eigenvalues are numbered in decreasing order of their real parts. The vector $v_1$ is then obtained from normalizing $\hat v_1$. The simplest choice for the coefficients $\rho_i$ is to take $\rho_i = 1, i=1,\ldots,p$. There are several drawbacks to this approach, the most important of which being that there is no easy way of choosing the coefficients $\rho_i$ in a systematic manner. The result is that for hard problems, convergence is difficult to achieve.

A more reliable alternative is to compute one eigenpair at a time and use deflation. The matrix $A$ can be deflated explicitly by constructing progressively the first $k$ Schur vectors. If a previous orthogonal basis $U_{k-1} = [u_1,\ldots,u_{k-1}]$ of the invariant subspace has already been computed, then to compute the eigenvalue $\lambda_{k}$, we can work with the matrix

$$\tilde A = A - U_{k-1} \Sigma U_{k-1}^{\ast},$$

in which $\Sigma = \diag(\sigma_i)$ is a diagonal matrix of shifts. The eigenvalues of $\tilde A$ consist of two groups. Those eigenvalues associated with the Schur vectors $u_1, \ldots, u_{k-1}$ will be shifted to $\tilde \lambda_i = \lambda_i - \sigma_i$ and the others remain unchanged. If the eigenvalues with largest real parts are sought, then the shifts are selected so that $\lambda_k$ becomes the next eigenvalue with largest real part of $\tilde A$. It is also possible to deflate by simply projecting out the components associated with the invariant subspace spanned by $U_{k-1}$; this would lead to operating with the matrix

$$\tilde A = A (I - U_{k-1} U_{k-1}^{\ast}).$$

Note that if $A U_{k-1} = U_{k-1} R_{k-1}$ is the partial Schur decomposition associated with the first $k-1$ Ritz values, then $\tilde A = A - U_{k-1} R_{k-1} U_{k-1}^{\ast}$. Those eigenvalues associated with the Schur vectors $u_1, \ldots, u_{k-1}$ will now all be moved to zero.

A better implementation of deflation, which fits in well with the Arnoldi procedure, is to work with a single basis $v_1, v_2, \ldots , v_m$ whose first vectors are the Schur vectors that have already converged. Suppose that $k-1$ such vectors have converged and call them $v_1, v_2 ,\ldots,v_{k-1}$. Then we start by choosing a vector $v_k$ which is orthogonal to $v_1,\ldots,v_{k-1}$ and of norm 1. Next we perform $m-k$ steps of an Arnoldi procedure in which orthogonality of the vector $v_j$ against all previous $v_i$'s, including $v_1,\ldots,v_{k-1}$ is enforced. This generates an orthogonal basis of the subspace

$${\rm span}\{ v_1,\ldots, v_{k-1} , v_k, A v_k, \ldots, A^{m-k} v_k \} \ .$$ (125)

Thus, the dimension of this modified Krylov subspace is constant and equal to $m$ in general. A sketch of this implicit deflation procedure combined with the Arnoldi method appears in the following.

Note that in the Loop, the Schur vectors associated with the eigenvalues $\la_1,\ldots,\la_{k-1}$ will not be touched in subsequent steps. They are sometimes referred to as ``locked vectors.'' Similarly, the corresponding upper triangular matrix corresponding to these vectors is also locked.

$$\underbrace{\left[v_1, v_2, \ldots, v_{k-1}\right.}_{Locked}, \underbrace{\left. v_k, v_{k+1}, \ldots v_m \right] }_{Active}$$

When a new Schur vector converges, step (10) computes the $k$th column of $R$ associated with this new basis vector. In the subsequent steps, the approximate eigenvalues are the eigenvalues of the $m \times m$ Hessenberg matrix $H_m$ defined in the algorithm and whose $k \times k$ principal submatrix is upper triangular. For example, when $m=6$ and after the second Schur vector, $k=2$, has converged, the matrix $H_m$ will have the form

$$H_m ~ = ~ \left[ \begin{array}{cccccc} * & * & * & * & * &... ... & & & * & * & * \\ & & & & * & * \\ \end{array} \right].$$ (126)

In the subsequent steps, only the eigenvalues not associated with the $2 \times 2$ upper triangular matrix need to be considered.

It can be shown that, in exact arithmetic, the $(n-k) \times (n-k)$ Hessenberg matrix in the lower $(2 \times 2)$ block is the same matrix that would be obtained from an Arnoldi run applied to the matrix

$$\tilde A = (I-U_{k-1} U_{k-1}^{\ast} ) A.$$

Thus, we are implicitly projecting out the invariant subspace already computed from the range of $A$.