from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

8.2.Power iteration¶

Given that matrix-vector multiplication is fast for sparse matrices, let’s see what we might accomplish with only that at our disposal.

Example 8.2.1 (Power iteration)

Here we choose a random 5×5 matrix and a random 5-vector.

A = random.choice(range(10), (5, 5))
A = A / sum(A, 0)
x = random.randn(5)
print(x)

[ 0.29471445 -0.42990771 -1.36506003 -0.07878796  0.73791809]

Applying matrix-vector multiplication once doesn’t do anything recognizable.

y = A @ x
print(y)

[-0.02922966 -0.27181066 -0.00437814 -0.36560917 -0.17009553]

Repeating the multiplication still doesn’t do anything obvious.

z = A @ y
print(z)

[-0.22142726 -0.17437522 -0.17932751 -0.15087488 -0.11511829]

But if we keep repeating the matrix-vector multiplication, something remarkable happens: $\mathbf{A} \mathbf{x} \approx \mathbf{x}$ .

x = random.randn(5)
for j in range(6):
    x = A @ x
print(x)
print(A @ x)

[1.69425975 1.25621914 1.04154581 0.66726285 1.01234703]
[1.69375606 1.25615596 1.04177167 0.66773407 1.01221682]

This phenomenon is unlikely to be a coincidence!

There was a little cheating in Example 8.2.1 to make the story come out neatly (specifically, the normalization step after creating a random matrix). But it illustrates an important general fact that we investigate now.

8.2.1Dominant eigenvalue¶

Analysis of matrix powers is most straightforward in the diagonalizable case. Let $\mathbf{A}$ be any diagonalizable $n\times n$ matrix having eigenvalues $\lambda_1,\ldots,\lambda_n$ and corresponding linearly independent eigenvectors $\mathbf{v}_1,\ldots,\mathbf{v}_n$ . For our later convenience (and without losing any generality), we assume that the eigenvectors are normalized so that

\twonorm{\mathbf{v}_j} = 1 \quad \text{for } j=1,\ldots,n.

(8.2.1)

We also make an important assumption about the eigenvalue magnitudes that will hold for most but not all matrices.

In Example 8.2.1, for instance, $\lambda_1=1$ is the dominant eigenvalue.

Now let $\mathbf{x}$ be an $n$ -vector, let $k$ be a positive integer, and refer to (7.2.11):

\mathbf{A}^k \mathbf{x} = \mathbf{V}\mathbf{D}^k\mathbf{V}^{-1}\mathbf{x}.

(8.2.3)

Let $\mathbf{z}=\mathbf{V}^{-1}\mathbf{x}$ , and recall that $\mathbf{D}$ is a diagonal matrix of eigenvalues. Then

\begin{split} \mathbf{A}^k\mathbf{x} &= \mathbf{V}\mathbf{D}^k \mathbf{z} = \mathbf{V}\begin{bmatrix} \lambda_1^kz_1 \\[0.5ex] \lambda_2^kz_2 \\ \vdots \\ \lambda_n^kz_n \end{bmatrix} \\ &= \lambda_1^k \left[ z_1 \mathbf{v}_{1} + z_2 \left(\frac{\lambda_2}{\lambda_1}\right) ^k \mathbf{v}_{2} + \cdots + z_n \left(\frac{\lambda_n}{\lambda_1}\right)^k \mathbf{v}_{n} \right]. \end{split}

(8.2.4)

Since $\lambda_1$ is dominant, we conclude that if $z_1\neq 0$ ,

\begin{split} \twonorm{ \frac{ \mathbf{A}^k\mathbf{x}}{z_1 \lambda_1^k} - \mathbf{v}_1 } & \le \underbrace{\left|\frac{z_2}{z_1}\right|}_{c_2} \cdot \underbrace{ \left| \frac{\lambda_2}{\lambda_1} \right|^k}_{r_2^k} \twonorm{\mathbf{v}_{2}} + \cdots + \underbrace{\left|\frac{z_n}{z_1}\right|}_{c_n} \cdot \underbrace{\left|\frac{\lambda_n}{\lambda_1}\right|^k}_{r_n^k} \twonorm{\mathbf{v}_{n}} \\[1mm] & = \sum_{j=2}^n c_j r_j^k \cdot 1 \\ & \rightarrow 0 \text{ as $k\rightarrow \infty$}, \end{split}

(8.2.5)

since, by (8.2.2), $r_j = | \lambda_j / \lambda_1 | < 1$ for $j=2,\ldots,n$ . If we choose $\mathbf{x}$ randomly, then the probability that $z_1=0$ is zero, so we will not be concerned with that case.

8.2.2Algorithm¶

An important technicality separates us from an algorithm: unless $|\lambda_1|=1$ , the factor $\lambda_1^k$ tends to make $\|\mathbf{A}^k\mathbf{x}\|$ either very large or very small. In practice, we cannot easily normalize by $\lambda_1^k$ as we did in (8.2.5), since we don’t know $\lambda_1$ in advance.

This issue is resolved by alternating matrix–vector multiplications with renormalizations.

Observe that we can write

\mathbf{x}_{k} = (\alpha_1 \alpha_2 \cdots \alpha_k ) \mathbf{A}^k \mathbf{x}_{1},

(8.2.6)

where $\alpha_j = \twonorm{\mathbf{y}_j}^{-1}$ is the normalization factor at iteration $j$ . Thus, the reasoning of (8.2.5) implies that $\mathbf{x}_k$ converges to a scalar multiple of the dominant eigenvector $\mathbf{v}_1$ of $\mathbf{A}$ . Specifically, suppose that there is some number $\gamma$ with $|\gamma|=1$ such that

\mathbf{x}_k - \gamma \mathbf{v}_1 = \epsilon \mathbf{w},

(8.2.7)

where $\norm{\mathbf{w}} = 1$ and $\epsilon \ll 1$ . Then

\begin{split} \beta_k &= \mathbf{x}_k^* \mathbf{y}_k \\ &= \mathbf{x}_k^* \mathbf{A} \mathbf{x}_k \\ &= \bigl(\overline{\gamma} \mathbf{v}_1^* + \epsilon \mathbf{w}^*\bigr) \mathbf{A} \bigl(\gamma \mathbf{v}_1 + \epsilon \mathbf{w}\bigr) \\ & = |\gamma|^2 \mathbf{v}_1^* (\lambda_1 \mathbf{v}_1) + O(\epsilon) \\ & = \lambda_1 + O(\epsilon). \end{split}

(8.2.8)

That is, $\beta_k$ from the power iteration estimates $\lambda_1$ about as well as $\mathbf{x}_k$ estimates the dominant eigenvector.

Function 8.2.1 is our implementation of power iteration.

Algorithm 8.2.1 (poweriter)

poweriter.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def poweriter(A, numiter):
    """
    poweriter(A, numiter)

    Perform numiter power iterations with the matrix A, starting from a random vector, 
    and return a vector of eigenvalue estimates and the final eigenvector approximation.
    """
    n = A.shape[0]
    x = np.random.randn(n)
    x = x / np.linalg.norm(x)
    beta = np.zeros(numiter)
    for k in range(numiter):
        y = A @ x
        beta[k] = np.dot(x, y)
        x = y / np.linalg.norm(y)

    return beta, x

8.2.3Convergence rate¶

While we now feel confident that the sequence $\{\beta_k\}$ converges to the dominant eigenvalue $\lambda_1$ , we would like to know how fast this convergence is. Looking back at (8.2.4), we know that our normalizations make it so that

\mathbf{x}_{k} - \gamma \mathbf{v}_1 = b_2 \left(\frac{\lambda_2}{\lambda_1}\right) ^k \mathbf{v}_{2} + \cdots + b_n \left(\frac{\lambda_n}{\lambda_1}\right)^k \mathbf{v}_{n}

(8.2.10)

for some constants $b_2,\ldots,b_n$ . If we now make a stronger assumption that $\lambda_2$ dominates the rest of the eigenvalues, i.e.,

|\lambda_1| > |\lambda_2| > |\lambda_3| \ge \cdots \ge |\lambda_n|,

(8.2.11)

then expression on the right-hand side of (8.2.10) is dominated by its first term, because

\sum_{j=2}^n b_j \left(\frac{\lambda_j}{\lambda_1}\right) ^k \mathbf{v}_{j} = \left(\frac{\lambda_2}{\lambda_1}\right)^k \left[ b_2 \mathbf{v}_2 + \underbrace{\sum_{j=3}^n {b_j} \left(\frac{\lambda_j}{\lambda_2}\right) ^k \mathbf{v}_{j} }_{\to 0\, \text{ as } k\to\infty} \right] .

(8.2.12)

Therefore, (8.2.8) now implies that $|\beta_k - \lambda_1|$ is dominated by $|\lambda_2 / \lambda_1|^k$ , which is a case of linear convergence.

Example 8.2.2 (Convergence of power iteration)

We will experiment with the power iteration on a 5×5 matrix with prescribed eigenvalues and dominant eigenvalue at 1.

ev = [1, -0.75, 0.6, -0.4, 0]
A = triu(ones([5, 5]), 1) + diag(ev)    # triangular matrix, eigs on diagonal

We run the power iteration 60 times. The first output should be a sequence of estimates converging to the dominant eigenvalue—which, in this case, we set up to be 1.

beta, x = FNC.poweriter(A, 60)
print(beta)

[0.55788943 1.72091387 1.33729825 1.15014057 1.0998233  1.03993647
 1.03451344 1.01168549 1.0132165  1.00316611 1.00539449 1.00060579
 1.00232517 0.99992339 1.00105605 0.99980723 1.00050385 0.99983773
 1.00025103 0.99988934 1.00012955 0.99993077 1.00006868 0.99995855
 1.00003713 0.99997578 1.00002034 0.99998605 1.00001125 0.99999204
 1.00000626 0.99999548 1.00000349 0.99999744 1.00000196 0.99999856
 1.0000011  0.99999919 1.00000062 0.99999954 1.00000035 0.99999974
 1.00000019 0.99999985 1.00000011 0.99999992 1.00000006 0.99999995
 1.00000003 0.99999997 1.00000002 0.99999999 1.00000001 0.99999999
 1.00000001 1.         1.         1.         1.         1.        ]

We check for linear convergence using a log-linear plot of the error.

err = 1 - beta
semilogy(arange(60), abs(err), "-o")
ylim(1e-10, 1)
xlabel("$k$")
ylabel("$|\\lambda_1 - \\beta_k|$")
title("Convergence of power iteration");

The asymptotic trend seems to be a straight line, consistent with linear convergence. To estimate the convergence rate, we look at the ratio of two consecutive errors in the linear part of the convergence curve. The ratio of the first two eigenvalues should match the observed rate.

print(f"theory: {ev[1] / ev[0]:.5f}")
print(f"observed: {err[40] / err[39]:.5f}")

theory: -0.75000
observed: -0.75371

Note that the error is supposed to change sign on each iteration. The effect of these alternating signs is that estimates oscillate around the exact value.

print(beta[26:30])

[1.00002034 0.99998605 1.00001125 0.99999204]

In practical situations, we don’t know the exact eigenvalue that the algorithm is supposed to find. In that case we would base errors on the final $\beta$ that was found, as in the following plot.

err = beta[-1] - beta
semilogy(arange(60), abs(err), "-o")
ylim(1e-10, 1), xlabel("$k$")
ylabel("$|\\lambda_1 - \\beta_k|$")
title("Convergence of power iteration");

The results are very similar until the last few iterations, when the limited accuracy of the reference value begins to show. That is, while it is a good estimate of $\lambda_1$ , it is less good as an estimate of the error in nearby estimates.

The practical utility of (8.2.13) is limited: if we knew $\lambda_1$ and $\lambda_2$ , we wouldn’t be running the power iteration in the first place! Sometimes it’s possible to find estimates of or bounds on the ratio. If nothing else, though, it is useful to know that linear convergence is expected at a rate based solely on the dominant eigenvalues.

8.2.4Exercises¶

Exercise 8.2.3

⌨ In Exercise 2.3.5 we considered a mass-lumped model of a hanging string that led to a tridiagonal system of linear equations. Then, in Exercise 7.2.6, we found that eigenvectors of the same matrix correspond to vibrational modes of the string. The same setup can be applied to a membrane hanging from a square frame. Lumping the mass onto a Cartesian grid, each interacts with the four neighbors to the north, south, east, and west. If $n$ masses are used in each coordinate direction, we get an $n^2\times n^2$ sparse matrix $\mathbf{A}$ that can be constructed by FNC.poisson(n).

(a) Let $n=10$ and make a spy plot of $\mathbf{A}$ . What is the density of $\mathbf{A}$ ? Most rows all have the same number of nonzeros; find this number.

(b) Find the dominant $\lambda_1$ using eigs for $n=10,15,20,25$ .

(c) For each $n$ in part (b), apply 100 steps of Function 8.2.1. On one graph, plot the four convergence curves $|\beta_k-\lambda_1|$ using a semi-log scale. (They will not be smooth curves because this matrix has many repeated eigenvalues that complicate the convergence analysis.)