clear all
format short
set(0, 'defaultaxesfontsize', 12)
set(0, 'defaultlinelinewidth', 1.5)
set(0, 'defaultFunctionLinelinewidth', 1.5)
set(0, 'defaultscattermarkerfacecolor', 'flat')
gcf;
set(gcf, 'Position', [0 0 600 350], 'Theme', 'light')
addpath ../FNC_matlab

Starting MATLAB ...

Executing ...

8.3.Inverse iteration¶

Power iteration finds only the dominant eigenvalue. We next show that it can be adapted to find any eigenvalue, provided you start with a reasonably good estimate of it. Some simple linear algebra is all that is needed.

Proof

The equation $\mathbf{A}\mathbf{v}=\lambda \mathbf{v}$ implies that $(\mathbf{A}-s\mathbf{I})\mathbf{v} = \mathbf{A}\mathbf{v} - s\mathbf{I}\mathbf{v} = \lambda\mathbf{v} - s\mathbf{v} = (\lambda-s)\mathbf{v}$ . That proves the first part of the theorem, and half of part 3.

For the rest, we note that by assumption, $(\mathbf{A}-s\mathbf{I})$ is nonsingular, so $(\mathbf{A}-s\mathbf{I})\mathbf{v} = (\lambda-s) \mathbf{v}$ implies that $\mathbf{v} = (\lambda-s)^{-1} (\mathbf{A}-s\mathbf{I}) \mathbf{v}$ , or $(\lambda-s)^{-1} \mathbf{v} =(\mathbf{A}-s\mathbf{I})^{-1} \mathbf{v}$ .

Consider first part 2 of the theorem with $s=0$ , and suppose that $\mathbf{A}$ has a smallest eigenvalue $\lambda_1$ , such that

|\lambda_n| \ge |\lambda_{n-1}| \ge \cdots > |\lambda_1|.

(8.3.1)

Then clearly

|\lambda_1^{-1}| > |\lambda_{2}^{-1}| \ge \cdots \ge |\lambda_n^{-1}|,

(8.3.2)

and $\mathbf{A}^{-1}$ has dominant eigenvalue $\lambda_1^{-1}$ . Hence, power iteration on $\mathbf{A}^{-1}$ can be used to find the eigenvalue of $\mathbf{A}$ closest to zero. For nonzero values of $s$ , then we suppose there is an ordering

|\lambda_n-s| \ge \cdots \ge |\lambda_2-s| > |\lambda_1-s|.

(8.3.3)

Then it follows that

|\lambda_1-s|^{-1} > |\lambda_{2}-s|^{-1} \ge \cdots \ge |\lambda_n-s|^{-1},

(8.3.4)

and power iteration on the matrix $(\mathbf{A}-s\mathbf{I})^{-1}$ converges to $(\lambda_1-s)^{-1}$ , which is easily solved for $\lambda_1$ itself.

8.3.1Algorithm¶

A literal application of Definition 8.2.2 would include the step

\mathbf{y}_k = (\mathbf{A}-s\mathbf{I})^{-1} \mathbf{x}_k.

(8.3.5)

As always, however, we do not want to explicitly find the inverse of a matrix. Instead, we should implement this step as the solution of a linear system.

Each pass of inverse iteration requires the solution of a linear system of equations with the matrix $\mathbf{B}=\mathbf{A}-s\mathbf{I}$ . This solution might use methods we consider later in this chapter. Here, we use (sparse) PLU factorization and hope for the best. Since the matrix $\mathbf{B}$ is constant, the factorization needs to be done only once for all iterations. The details are in Function 8.3.1.

Algorithm 8.3.1 (inviter)

inviter.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
function [beta, x] = inviter(A, s, numiter)
% INVITER   Shifted inverse iteration for the closest eigenvalue.
% Input:
%   A         square matrix
%   s         value close to targeted eigenvalue (complex scalar)
%   numiter   number of iterations
% Output: 
%   beta      sequence of eigenvalue approximations (vector)
%   x         final eigenvector approximation

    n = length(A);
    x = randn(n, 1);
    x = x / norm(x, inf);
    B = A - s * eye(n);
    [L, U] = lu(B);
    beta = zeros(numiter, 1);
    for k = 1:numiter
        y = U \ (L \ x);
        beta(k) = (1 / dot(x, y)) + s;
        x = y / norm(y);
    end
end

8.3.2Convergence rate¶

The convergence is linear, at a rate found by reinterpreting Theorem 8.2.1 with $(\mathbf{A}-s\mathbf{I})^{-1}$ in place of $\mathbf{A}.$ With the eigenvalues ordered as in (8.3.3), in the general case we have

\frac{\abs{\beta_{k+1} - \lambda_1}}{\abs{\beta_{k} - \lambda_1}} \rightarrow \abs{\frac{ \lambda_1 - s } {\lambda_2 - s}}\quad \text{ as } \quad k\rightarrow \infty,

(8.3.7)

and in the hermitian case, we have

\frac{\abs{\beta_{k+1} - \lambda_1}}{\abs{\beta_{k} - \lambda_1}} \rightarrow \abs{\frac{ \lambda_1 - s } {\lambda_2 - s}}^2 \quad \text{ as } \quad k\rightarrow \infty.

(8.3.8)

Thus, the convergence is best when the shift $s$ is close to the target eigenvalue—specifically, when it is much closer to that eigenvalue than to any other.

Example 8.3.1 (Convergence of inverse iteration)

We set up a $5\times 5$ triangular matrix with prescribed eigenvalues on its diagonal.

ev = [1, -0.75, 0.6, -0.4, 0];
A = triu(ones(5, 5), 1) + diag(ev);

We run inverse iteration with the shift $s=0.7$ . The result should converge to the eigenvalue closest to 0.7, which we know to be 0.6 here.

s = 0.7;
[beta, x] = inviter(A, s, 30);
format short
beta(1:10)

The convergence is again linear.

err = 0.6 - beta;
semilogy(abs(err),'.-')
title('Convergence of inverse iteration')
xlabel('k'), ylabel(('|\lambda_j - \beta_k|'));

Let’s reorder the eigenvalues to enforce (8.3.3).

Tip

The second output of sort returns the index permutation needed to sort the given vector.

[~, idx] = sort(abs(ev - s));
ev = ev(idx)

Now it is easy to compare the theoretical and observed linear convergence rates.

theoretical_rate = (ev(1) - s) / (ev(2) - s)
observed_rate = err(26) / err(25)

8.3.3Rayleigh quotient iteration¶

There is a clear opportunity for positive feedback in Definition 8.3.1. The convergence rate of inverse iteration improves as the shift gets closer to the true eigenvalue—and the algorithm computes improving eigenvalue estimates! Updating the shift to $s=\beta_k$ after each iteration greatly accelerates the convergence. You are asked to implement this algorithm in Exercise 8.3.6.

If the eigenvalues are ordered by distance to $s$ , then (asymptotically) one step of inverse iteration reduces the error by the factor $|\lambda_1-s|/|\lambda_2-s|$ . As $s \to\lambda_1$ , the change in the denominator is negligible. So, if at one point the error $\abs{\lambda_1-s}$ is about $\epsilon$ , then the error in the next estimate is reduced by a factor $O(\epsilon)$ , making it $O(\epsilon^2)$ . That is, each step now squares the error, which is quadratic convergence.

Example 8.3.2 (Dynamic shift strategy)

ev = [1, -0.75, 0.6, -0.4, 0];
A = triu(ones(5, 5), 1) + diag(ev);

We begin with a shift $s=0.7$ , which is closest to the eigenvalue 0.6.

s = 0.7;
x = ones(5, 1);
y = (A - s * eye(5)) \ x;
beta = 1 / (x' * y) + s

Note that the result is not yet any closer to the targeted 0.6. But we proceed (without being too picky about normalization here).

s = beta;
x = y / norm(y);
y = (A - s * eye(5)) \ x;
beta = 1 / (x' * y) + s

Still not much apparent progress. However, in just a few more iterations the results are dramatically better.

format long
for k = 1:4
    s = beta;
    x = y / norm(y);
    y = (A - s * eye(5)) \ x;
    beta = 1 / (x' * y) + s
end

There is a price to pay for this improvement. The matrix of the linear system to be solved, $(\mathbf{A}-s\mathbf{I}),$ now changes with each iteration. That means that we can no longer do just one LU factorization for the entire iteration. The speedup in convergence usually makes this tradeoff worthwhile, however.

In practice power and inverse iteration are not as effective as the algorithms used by eigs and based on the mathematics described in the rest of this chapter. However, inverse iteration can be useful for turning an eigenvalue estimate into an eigenvector estimate.

8.3.4Exercises¶

Exercise 8.3.1

⌨ Use Function 8.3.1 to perform 10 iterations for the given matrix and shift. Compare the results quantitatively to the convergence given by (8.3.7).

(a) $\mathbf{A} = \begin{bmatrix} 1.1 & 1 \\ 0 & 2.1 \end{bmatrix}, \; s = 1 \qquad$ (b) $\mathbf{A} = \begin{bmatrix} 1.1 & 1 \\ 0 & 2.1 \end{bmatrix}, \; s = 2\qquad$

(c) $\mathbf{A} = \begin{bmatrix} 1.1 & 1 \\ 0 & 2.1 \end{bmatrix}, \; s = 1.6\qquad$ (d) $\mathbf{A} = \begin{bmatrix} 2 & 1 \\ 1 & 0 \end{bmatrix}, \; s = -0.33 \qquad$

(e) $\mathbf{A} = \begin{bmatrix} 6 & 5 & 4 \\ 5 & 4 & 3 \\ 4 & 3 & 2 \end{bmatrix}, \; s = 0.1$

Exercise 8.3.4

✍ When the shift $s$ is very close to an eigenvalue of $\mathbf{A}$ , the matrix $\mathbf{A}-s\mathbf{I}$ is close to a singular matrix. But then (8.3.6) is a linear system with a badly conditioned matrix, which should create a lot of error in the numerical solution for $\mathbf{y}_k$ . However, it happens that the error is mostly in the direction of the eigenvector we are looking for, as the following toy example illustrates.

Prove that $\displaystyle \begin{bmatrix} 1 & 1 \\ 0 & 0 \end{bmatrix}$ has an eigenvalue at zero with associated eigenvector $\mathbf{v}=[-1,1]^T$ . Suppose this matrix is perturbed slightly to $\displaystyle \mathbf{A} = \begin{bmatrix} 1 & 1 \\ 0 & \epsilon \end{bmatrix}$ , and that $\mathbf{x}_k=[1,1]$ in (8.3.6). Show that once $\mathbf{y}_k$ is normalized by its infinity norm, the result is within $\epsilon$ of a multiple of $\mathbf{v}$ .

Exercise 8.3.5

⌨ (Continuation of Exercise 8.2.3.) This exercise concerns the $n^2\times n^2$ sparse matrix defined by FNC.poisson(n) for integer $n$ . It represents a lumped model of a vibrating square membrane held fixed around the edges.

(a) The eigenvalues of $\mathbf{A}$ closest to zero are approximately squares of the frequencies of vibration for the membrane. Using eigs, find the eigenvalue $\lambda_m$ closest to zero for $n=10,15,20,25$ .

(b) For each $n$ in part (a), apply 50 steps of Function 8.3.1 with zero shift. On one graph, plot the four convergence curves $|\beta_k-\lambda_m|$ using a semi-log scale.

(c) Let v be the eigenvector (second output) found by Function 8.3.1 for $n=25$ . Make a surface plot of the vibration mode by reshaping v into an $n\times n$ matrix.

Exercise 8.3.6

⌨ This problem explores the use of Rayleigh quotient iteration.

(a) Modify Function 8.3.1 to change the value of the shift $s$ to be the most recently computed value in the vector $\beta$ . Note that the matrix B must also change with each iteration, so the LU factorization cannot be done just once.

(b) Define a $100\times 100$ matrix with values $k^2$ for $k=1,\ldots,100$ on the main diagonal and random values uniformly distributed between 0 and 1 on the first superdiagonal. (Since this matrix is triangular, the diagonal values are its eigenvalues.) Using an initial shift of $s=920$ , apply Rayleigh quotient iteration. Determine which eigenvalue was found and make a table of the log10 of the errors in the iteration as a function of iteration number. (These should approximately double, until machine precision is reached, due to quadratic convergence.)

(c) Repeat part (b) using a different initial shift of your choice.