102

100

10−2

10−4

10−6

10−8

10−10

Figure 6.6: Simple vector iteration with τI40 − T40 and starting vector (1, 1, . . . , 1)T

Problem 6.6 When computing the smallest eigenvalue of Tn by the simple vector iteration we can find a better shift than τ above if the extremal points of the spectrum are known. Determine σ such that σIn − Tn exhibits the optimal convergence rate. Hint: On the one hand we would like the quotient (σ −λn−1)/(σ −λn) to be as small as possible. On the other hand |σ −λ1|/(σ −λn) must not become to big. Hint: Equate the two quantities.

6.5Inverse vector iteration

The previous examples have shown that the convergence of simple vector iteration is potentially very slow. The quotient of the second largest to the largest eigenvalue are very close to 1. We noticed this by means of a very simple and small eigenvalue problem. The situation gets much worse if the problems are big.

We have seen in (2.28) that a polynomial in A has the same eigenvectors as A. We therefore may try to find a polynomial that enhances the eigenvalue that we are looking for. This approach is however not successful in the most critical case when the wanted eigenvalue is very close to unwanted. In this situation, the shift-and-invert spectral transformation is most appropriate. Instead of a polynomial we transform the matrix by the rational function f(λ) = 1/(λ − σ) where σ is a so-called shift that is chosen close to the desired eigenvalue. Simple vector iteration with the matrix (A − σI)−1 is referred to as inverse vector iteration, see Algorithm 6.5.

The iteration converges towards the eigenvector with eigenvalue closest to σ. A linear system of equations has to be solved in each iteration step. Of course only one Cholesky or LU factorization has to be computed as the shift remains constants in all iterations. The stopping criterion is changed into

Table 6.1: Computing the lowest eigenvalue of the one-dimensinal Poisson equation by inverse iteration

known eigenvalue. The matrix that we consider is one of the so-called Wilkinson matrices

Wilkinson matrices are irreducible tridiagonal matrices that have very close eigenvalues. This matrix has the eigenvalues

The following MATLAB code constructs the SPARSE tridiagonal matrix T .

n = 40;

e = ones(n,1); f = abs([-n/2+1:n/2]’); T = spdiags([-e f -e], [-1:1], n, n); lam = sort(eig(T));

Computing the 20-th and 21-st eigenvectors could be done in the following way.

>>x = (T - lam(20)*eye(n))\e;

>>y = (T - lam(21)*eye(n))\e;

>>x = x/norm(x); y = y/norm(y);

>>x’*y

ans = 0.00140329005834

>> norm((T - lam(20)*eye(n))*x)

ans =

7.325760095786749e-15

>> norm((T - lam(21)*eye(n))*y) ans =

7.120036319503636e-15

The computed vectors x and y are good approximations in the sense that they give small residuals. However, the two vectors are nor mutually orthogonal at all. We try to improve orthogonality by apply a second step of inverse iteration

>>x = (T - lam(20)*eye(n))\x;

>>y = (T - lam(21)*eye(n))\y;

>>x = x/norm(x); y = y/norm(y);

>>x’*y

ans = -1.313592004487587e-05

Things have only slightly improved. Therefore, we orthogonalize y explicitely against x.

>>y = y - x*(x’*y);

>>x’*y

ans = -2.155571068436496e-17

>>norm((T - lam(21)*eye(n))*y) ans =

1.557058217172078e-15

>>norm((T - lam(20)*eye(n))*x) ans =

4.117116818055497e-16

This helped. The two eigenvectors are now perpendicular on each other, and the residuals are still fine.

Discussion of inverse iteration

We have seen that

•we can compute eigenvectors corresponding to any (simple and well separated) eigenvalue if we choose the shift properly, and that

•we have very good convergence rates, is the shift is close to an eigenvalue.

However, one may feel uncomfortable solving an ‘almost singular’ system of equation, after all σ ≈ λk means that the condition of A − σI is very big. From the analysis of linear systems of equations we know that this means large errors in the solution. Furtunately, the error that is su ered from when solving with A − σI points in the right direction. To see this, assume that the singular value decomposition of A − σI is given by

A − σI = UΣV , Σ = diag(σ1, . . . , σn), with σ1 ≥ · · · ≥ σn ≥ 0.

If A − σI is ‘almost’ singular then σn 1. If even σn = 0 then (A − σI)vn = 0, i.e., the last right singular vector is an eigenvector of A corresponding to the eigenvalue σ (the shift).