Power and inverse iterations#

import numpy as np

Power iteration#

We slightly re-arrange the Algorithm 27.1 so as to perform only one matrix-vector product per iteration.

# Algorithm 27.1
def power_iter(A, v, maxiter=1000, tol=1.0e-12):
    v /= np.linalg.norm(v)
    w = A @ v
    lam = np.dot(v, w)
    print(0, lam)
    for k in range(maxiter):
        v = w /np.linalg.norm(w)
        w = A @ v
        lam = np.dot(v, w)
        print(k+1, lam)
        if np.linalg.norm(A@v - lam*v) < tol:
            return v, lam
    print('Did not converge !!!')
    return v, lam

Example 27.1#

\[\begin{split} A = \begin{bmatrix} 2 & 1 & 1 \\ 1 & 3 & 1 \\ 1 & 1 & 4 \end{bmatrix} \end{split}\]

The eigenvalues as given by Matlab are

1.324869129433354
2.460811127189110
5.214319743377535

We start with

\[\begin{split} v = \frac{1}{\sqrt{3}} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \end{split}\]

The eigenvector closest to \(v\) has eigenvalue \(\lambda = 5.214319743377...\)

A = np.array([[2, 1, 1],
              [1, 3, 1],
              [1, 1, 4]])
v = np.array([1, 1, 1]) / np.sqrt(3.0)
v, lam = power_iter(A, v)
0 5.000000000000002
1 5.181818181818181
2 5.208192771084338
3 5.213028887981392
4 5.214037052110615
5 5.21425709431699
6 5.214305810348445
7 5.214316641508852
8 5.214319052610873
9 5.214319589534747
10 5.214319709113864
11 5.214319735746316
12 5.214319741677905
13 5.214319742998992
14 5.214319743293227
15 5.2143197433587565
16 5.214319743373353
17 5.214319743376604
18 5.2143197433773265
19 5.214319743377489
20 5.214319743377524
21 5.214319743377533
22 5.214319743377534
23 5.214319743377534
24 5.214319743377534
25 5.214319743377535
26 5.214319743377535
27 5.214319743377535
28 5.214319743377535
29 5.214319743377534
30 5.214319743377534
31 5.214319743377534
32 5.214319743377536
33 5.214319743377536
34 5.214319743377535
35 5.214319743377535
36 5.214319743377535
37 5.214319743377534

The eigenvalues converge quickly due to quadratic convergence, but eigenvectors converge slowly due to linear convergence.

Inverse iteration#

# Algorithm 27.2
def inv_iter(A, mu, maxiter=1000, tol=1.0e-12):
    v = 2 * np.random.rand(A.shape[0]) - 1
    v /= np.linalg.norm(v)
    lam = np.dot(v, A @ v)
    print(0, lam)
    for k in range(maxiter):
        B = A - mu * np.eye(A.shape[0])
        w = np.linalg.solve(B, v)
        v = w /np.linalg.norm(w)
        lam = np.dot(v, A @ v)
        print(k+1, lam)
        if np.linalg.norm(A@v - lam*v) < tol:
            return v, lam
    print('Did not converge !!!')
    return v, lam

We apply this on Example 27.1 starting with a random initial \(v\) and \(\mu = 5\).

A = np.array([[2, 1, 1],
              [1, 3, 1],
              [1, 1, 4]])
mu = 5.0
v, lam = inv_iter(A, mu)
0 2.3514908932274903
1 4.965915149671642
2 5.212582365992642
3 5.214308024559074
4 5.214319662155532
5 5.214319742806599
6 5.214319743373493
7 5.214319743377507
8 5.214319743377534
9 5.214319743377535
10 5.214319743377535
11 5.214319743377535
12 5.214319743377536
13 5.214319743377536

Rayleigh quotient iteration#

# Algorithm 27.3
def ray_iter(A, v=None, lam=None, maxiter=1000, tol=1.0e-12):
    m = A.shape[0]
    if v is None:
        v = 2 * np.random.random(m) - 1
    v /= np.linalg.norm(v)
    if lam is None:
        lam = np.dot(v, A @ v)
    print(0, lam)
    for k in range(maxiter):
        B = A - lam * np.eye(A.shape[0])
        w = np.linalg.solve(B, v)
        v = w /np.linalg.norm(w)
        lam = np.dot(v, A @ v)
        print(k+1, lam)
        if np.linalg.norm(A@v - lam*v) < tol:
            return v, lam
    print('Did not converge !!!')
    return v, lam

We apply this on Example 27.1.

A = np.array([[2, 1, 1],
              [1, 3, 1],
              [1, 1, 4]])
v = np.array([1, 1, 1]) / np.sqrt(3.0)
v, lam = ray_iter(A, v=v)
0 5.000000000000002
1 5.213114754098361
2 5.214319743184031
3 5.214319743377534

Let us try the same problem with a random initial \(v\).

v = 2 * np.random.rand(3) - 1
v, lam = ray_iter(A, v=v)
0 1.5046591617491027
1 1.326330044406299
2 1.3248691312831105
3 1.3248691294333537

Depending on which eigenvector is closest to the starting \(v\), we may converge to any of the three eigenvalues.

Specify the eigenvalue estimate

v, lam = ray_iter(A, lam=1.0)
0 1.0
1 1.344496753950572
2 1.3248741863822042
3 1.324869129433354
4 1.324869129433354
v, lam = ray_iter(A, lam=3.0)
0 3.0
1 2.446611000745533
2 2.4608086355417433
3 2.4608111271891113
4 2.460811127189111
v, lam = ray_iter(A, lam=5.0)
0 5.0
1 3.06835792546707
2 2.5740834268863555
3 2.458812422425412
4 2.4608111184188464
5 2.460811127189111

Choosing the shift close to an eigenvalue may not always converge to that eigenvalue, it also depends on the initial vector.