Algorithm 6.4 Inverse vector iteration for computing (λj , uj )

1:The LU factorization of A − σI: LU = P (A − σI) and the eigenvectors u1, . . . , uj−1 are known.

2:Choose a starting vector x(0) such that uq x(0) = 0, q < j.

3:Set k := 0.

4:while kx(k) − y(k)/µ(k)k > tolky(k)k do

5:k := k + 1;

6:x(k) := y(k−1)/ky(k−1)k;

7:y(k) := U−1L−1P x(k);

6.8Rayleigh quotient iteration

We now assume that the matrix the eigenpairs of which we want to determine is Hermitian (or symmetric).

We have noticed that inverse iteration is an e ective way to compute eigenpairs, if a good approximation of the desired eigenvalue is known. This approximation is used as a shift. However, as we have seen earlier, if a good approximation of an eigenvector is available its Rayleigh quotient gives a very good approximation of its eigenvalue.

Indeed we have the following

Lemma 6.10 Let q be any nonzero vector. The number ρ that minimizes kAq − ρqk is the Rayleigh quotient

Proof. Let ρ R be the Rayleigh quotient (6.26) of q 6= 0 and let τ C be any number. Then we have

kAq − (ρ + τ)qk2 = q A2q − (2ρ + τ + τ¯) q Aq + |ρ + τ|2q q

= q A2q − 2ρ q Aq − 2Re(τ) q Aq + ρ2 q q + 2ρRe(τ) q q + |τ|2q q

= q A2q − (q Aq)2 + |τ|2 q q. q q

The last term is smallest if τ = 0.

The following algorithm 6.8 is a modification of inverse iteration. In each iteration step the shift is modified to be the Rayleigh quotient of the most recent eigenvector approximation. This is not a curse but a blessing [3] as we have seen in section 6.5.

In step 4 of this algorithm a close to singular system of equation is solved. This results in a very long solution whose norm is used a the convergence criterion.

The Rayleigh quotient iteration usually converges, however not always towards the desired solution. Therefore, to investigate the convergence rate we make the following

Assumption: y(k) −−−→ x with Ax = λx.

k→∞

This assumption garantees that there is at all convergence towards a certain eigenvector x. Let kxk = ky(k)k = 1 and let the angle between this eigenvector and its approximation be ϕ(k) = (x, y(k)). Then the assumption implies that {ϕ(k)}∞k=1 converges to zero. We can write

y(k) = x cos ϕ(k) + u(k) sin ϕ(k), kxk = ky(k)k = ku(k)k = 1.

Algorithm 6.5 Rayleigh quotient iteration (RQI)

1:Choose a starting vector y0 Fn with ky0k = 1 and a tolerance ε.

2:for k = 1, 2, . . . do

3:ρ(k) := y(k−1) Ay(k−1).

4: Solve (A − ρ(k)I)z(k) = y(k−1) for z(k).

5:σ(k) = kz(k)k.

6:y(k) := z(k)/σ(k).

7:if σ(k) > 10/ε then

8:return {y(k)}

9:end if

10:end for

Let

y(k) Ay(k)

ρ(k) = ρ(y(k)) = y(k) y(k) = y(k) Ay(k)

λ

=λ(1 − cos2 ϕk) − sin2 ϕkρ(uk)

=(λ − ρ(uk)) sin2 ϕk.

We now prove the

Theorem 6.11 (Local convergence of Rayleigh quotient iteration) With the above assumption we have

tan ϕk+1 = sin ϕk+1/cos ϕk+1

=sin ϕk k(A − ρkI)−1ukk (λ − ρk)/ cos ϕk

=(λ − ρk) k(A − ρkI)−1ukk tan ϕk

=(λ − ρ(uk)) k(A − ρkI)−1ukk sin2 ϕk tan ϕk.

kukk2=1

We define the gap between the eigenvalue λ and the rest of A’s spectrum by

γ := min|λi − ρk|.

λi=6 λ

The assumption implies that there must be a k0 N such that |λ − ρk| < γ2 for all k > k0, and, therefore,

Because tan ϕk ≈ sin ϕk ≈ ϕk if ϕk 1 we can deduce the cubic convergence rate.

We now look at the sequence {uk} more closely. We note from (6.28) that this sequence is obtained by “inverse iteration with variable shift ρk”. But since ρk −→ λ with a cubic rate of convergence we can for large k assume that ρk = λ and that uk x.

We now consider two cases, either {uk} converges, or it does not converge.

1. We assume that {uk} converges. Then the limit vector uˆ must be an eigenvector of

{ } ˆ

A in span x . (In general, uˆ is an eigenvector corresponding to the eigenvalue λ that is closest to λ.) Thus,

k→∞

2.Now we assume that {uk} does not converge. Then A has two eigenvalues of equal distance to λ and the cluster points of the sequence {uk} are two vectors in the plane that is spanned by two eigenvectors corresponding to these two eigenvalues λ ± δ, αxp + βxq, where α 6= 0, β 6= 0, and α2 + β2 = 1. Their Rayleigh quotients are

ρ(αxp + βxq) = α2λp + β2λq = α2(λ ± δ) + β2(λ δ) = λ ± δ(α2 − β2).

As k −→ ∞ the Rayleigh quotients of uk jump between these two values. Hence,