---
# Section 5.3: The Power Method and Some Simple Extensions
---

Let $A \in \mathbb{C}^{n \times n}$ be a matrix with _linearly independent_ eigenvectors

$$
v_1, \ldots, v_n
$$

and corresponding eigenvalues

$$
\lambda_1, \ldots, \lambda_n
$$

(i.e., $A v_i = \lambda_i v_i$, for $i=1,\ldots,n$) ordered such that

$$
|\lambda_1| \ge |\lambda_2| \ge \cdots \ge |\lambda_n|.
$$

We say that $A$ has a **dominant eigenvalue** if
 
$$
|\lambda_1| > |\lambda_2|.
$$

---

## The Power Method

The basic idea of the **power method** is to pick a vector $q \in \mathbb{C}^n$ and compute the sequence

$$
q,\ A q,\ A^2 q,\ A^3 q,\ \ldots.
$$

Since the eigenvectors $v_1,\ldots,v_n$ form a basis for $\mathbb{C}^n$, we have that

$$
q = c_1 v_1 + \cdots + c_n v_n.
$$

For a random $q$, we expect $c_1 \ne 0$.

Then

$$
\begin{align}
A q
&= c_1 A v_1 + \cdots + c_n A v_n \\
&= c_1 \lambda_1 v_1 + \cdots + c_n \lambda_n v_n
\end{align}
$$

and

$$
\begin{align}
A^2 q
&= c_1 \lambda_1 A v_1 + \cdots + c_n \lambda_n A v_n \\
&= c_1 \lambda_1^2 v_1 + \cdots + c_n \lambda_n^2 v_n.
\end{align}
$$

In general, we have

$$
A^j q = c_1 \lambda_1^j v_1 + \cdots + c_n \lambda_n^j v_n
$$

and

$$
\frac{A^j q}{\lambda_1^j} = c_1 v_1 + c_2 \left(\frac{\lambda_2}{\lambda_1}\right)^j v_2 + \cdots + c_n \left(\frac{\lambda_n}{\lambda_1}\right)^j v_n.
$$

Letting

$$
q_j = \frac{A^j q}{\lambda_1^j},
$$

we have

$$
\begin{align}
\| q_j - c_1 v_1 \|
&= \left\| c_2 \left(\frac{\lambda_2}{\lambda_1}\right)^j v_2 + \cdots + c_n \left(\frac{\lambda_n}{\lambda_1}\right)^j v_n \right\| \\
&\le |c_2| \left|\frac{\lambda_2}{\lambda_1}\right|^j \|v_2\| + \cdots + |c_n| \left|\frac{\lambda_n}{\lambda_1}\right|^j \|v_n\| \\
&\le \left|\frac{\lambda_2}{\lambda_1}\right|^j \big(|c_2| \|v_2\| + \cdots + |c_n| \|v_n\|\big).
\end{align}
$$

Now suppose $|\lambda_1| > |\lambda_2|$. Then

$$
\left|\frac{\lambda_2}{\lambda_1}\right| < 1.
$$

Therefore,

$$
\left|\frac{\lambda_2}{\lambda_1}\right|^j \to 0 \quad \text{as} \ j \to \infty.
$$

Thus, $\| q_j - c_1 v_1 \| \to 0$ as $j \to \infty$, so we conclude that

$$
q_j \to c_1 v_1 \quad \text{as $j \to \infty$.}
$$

The rate of the convergence of the power method is generally linear ($\|q_{j+1} - c_1 v_1\| \approx r \|q_j - c_1 v_1\|$ for all $j$ sufficiently large) with convergence ratio

$$
r = \left|\frac{\lambda_2}{\lambda_1}\right|.
$$

Thus, the larger the gap between $|\lambda_1|$ and $|\lambda_2|$, the smaller the convergence ratio and the faster the convergence.

---

## Scaling

Since we usually do not know $\lambda_1$ while running the power method, we will not be able to compute $q_j = A^j q/\lambda_1^j$. However, it is important that we scale $A^j q$ since $\|A^j q\| \to \infty$ if $|\lambda_1| > 1$ and $\|A^j q\| \to 0$ if $|\lambda_1| < 1$.

A simple choice is to scale $A^j q$ so that its largest entry is equal to one. Thus, we let

$$
q_{j+1} = \frac{A q_j}{s_{j+1}},
$$

where $s_{j+1}$ is the component of $A q_j$ which has the largest absolute value.

---

## Algorithm

Given $q_0 \in \mathbb{C}^n$, we iterate

1. $\hat{q} = A q_j$

2. $s_{j+1} =$ entry of $\hat{q}$ with largest absolute value

3. $q_{j+1} \gets \hat{q}/s_{j+1}$

for $j = 0, 1, 2, \ldots$.

Then $q_j$ approaches a multiple of $v_1$ and $s_j$ approaches the eigenvalue $\lambda_1$.

If $A$ is a dense $n \times n$ matrix, then each iteration of this algorithm will require $2n^2 + O(n)$ flops. However, if $A$ is sparse and has at most $k$ nonzeros on each row, then each iteration will require approximately $2 k n$ flops. Therefore, the power method is very well suited for computing the dominant eigenvalue and associated eigenvector of large sparse matrices.

---

## `power_method`


```julia
using LinearAlgebra, SparseArrays
```


```julia
function scale!(q)
    maxval, idx = maximum((abs(q[i]),i) for i=1:length(q))
    s = q[idx]
    q ./= s
    return s
end
```


```julia
function power_method(A; tol=sqrt(eps())/2, maxiter=100_000)
    m, n = size(A)
    n == m || error("Matrix must be square.")
    
    q = randn(n)
    s = scale!(q)
    lam = s

    qold = similar(q)

    k = 0
    done = false
    while !done && k < maxiter
        k += 1
        copy!(qold, q)        # qold = q
        mul!(q, A, qold)      # q = A*qold
        s = scale!(q)         # q = q/s
        lam = s
        done = norm(A*q - lam*q)/norm(q) <= tol
    end

    if done
        println("Converged after $k iterations.")
    else
        println("Failed to converge.")
    end

    return lam, q
end
```


```julia
n = 1_000
k = 10
density = (k - 1)/n    # density = (k*n - n)/n^2

A = triu(sprand(n, n, density), 1)
A = A + A' + I

# Expect nnz(A) ≈ k*n
@show nnz(A)

if n <= 1000
    λ = eigvals(Matrix(A))
    abseig = abs.(λ) |> sort
    r = abseig[end-1]/abseig[end]
    @show r
end

println()
@time lam, q = power_method(A)
@show lam
@show norm(A*q - lam*q)/norm(q);
```

---

## Google PageRank Algorithm

Google uses its [PageRank](https://en.wikipedia.org/wiki/PageRank) algorithm to determine its ranking of webpages in search results.

The [Google matrix](https://en.wikipedia.org/wiki/Google_matrix) represents how webpages on the Internet link to one another.

PageRank uses the power method to compute the dominant eigenvector of the Google matrix, and this dominant eigenvector is then used to rank the importance of webpages.

By design, the convergence ratio of the Google matrix is

$$
\left|\frac{\lambda_2}{\lambda_1}\right| = 0.85,
$$

so the number of power method iterations is reasonable.

---

## The Inverse Power Method

Let $A \in \mathbb{C}^{n \times n}$ be nonsingular. Since $A$ is nonsingular, all of its eigenvalues are nonzero. 

Since

$$
A v = \lambda v \quad \implies \quad A^{-1} v = \lambda^{-1} v,
$$

the eigenvalues of $A^{-1}$ are $\lambda_1^{-1},\ldots,\lambda_n^{-1}$ and the corresponding eigenvectors are $v_1,\ldots,v_n$.

Since

$$
|\lambda_1| \ge |\lambda_2| \ge \cdots \ge |\lambda_n|,
$$

we have that

$$
\left|\lambda_1^{-1}\right| \le \left|\lambda_2^{-1}\right| \le \cdots \le \left|\lambda_n^{-1}\right|.
$$

If $|\lambda_{n-1}| > |\lambda_n|$, then $\left|\lambda_n^{-1}\right| >  \left|\lambda_{n-1}^{-1}\right|$, so the **inverse power method**,

$$
q,\ A^{-1} q,\ A^{-2} q,\ A^{-3} q,\ \ldots,
$$

will generate a sequence $q_j$ that converges to a multiple of $v_n$ (i.e., the eigenvector corresponding to the _smallest_ eigenvalue of $A$).

---

## `inverse_power_method`


```julia
function inverse_power_method(A; tol=sqrt(eps())/2, maxiter=100_000)
    m, n = size(A)
    n == m || error("Matrix must be square.")
    
    F = lu(A)
    
    q = randn(n)
    s = scale!(q)
    lam = 1/s

    qold = similar(q)

    k = 0
    done = false
    while !done && k < maxiter
        k += 1
        copy!(qold, q)        # qold = q
        ldiv!(q, F, qold)     # q = F\qold
        s = scale!(q)         # q = q/s
        lam = 1/s
        done = norm(A*q - lam*q)/norm(q) <= tol
    end

    if done
        println("Converged after $k iterations.")
    else
        println("Failed to converge.")
    end

    return lam, q
end
```


```julia
n = 1000
k = 5
density = (k - 1)/n    # density = (k*n - n)/n^2

A = triu(sprand(n, n, density), 1)
A = A + A' + I

# Expect nnz(A) ≈ k*n
@show nnz(A)

if n <= 1000
    λ = eigvals(Matrix(A))
    abseig = abs.(λ) |> sort
    r = abseig[1]/abseig[2]
    @show r
end

println()
@time lam, q = inverse_power_method(A)
@show lam
@show norm(A*q - lam*q)/norm(q);
```

---

## The Shift-and-Invert Method

If $A v = \lambda v$, then

$$
\big( A - \rho I \big) v = \big( \lambda - \rho \big) v.
$$

Therefore, using the inverse power method on $A - \rho I$, we can compute an eigenvector with eigenvalue closest to the shift $\rho$.

That is, if

$$
|\lambda_i - \rho| \ll |\lambda_j - \rho|, \quad \forall j \ne i,
$$

then the **shift-and-invert method**,

$$
q,\ (A - \rho I)^{-1} q,\ (A - \rho I)^{-2} q,\ (A - \rho I)^{-3} q,\ \ldots,
$$

will generate a sequence $q_j$ that converges to a multiple of $v_i$.

The rate of convergence is

$$
\left| \frac{\lambda_i - \rho}{\lambda_k - \rho} \right|,
$$

where $\lambda_k - \rho$ is the second smallest eigenvalue of $A - \rho I$ in absolute value.

Once we have an $LU$ decomposition of $A - \rho I$ (which requires $\frac{2}{3}n^3 + O(n^2)$ flops), we can compute

$$
q \gets (A - \rho I)^{-1} q
$$

each iteration in $2 n^2$ flops.

---

## `inverse_power_method` with shift $\rho$


```julia
function inverse_power_method(A; ρ=0.0, tol=sqrt(eps())/2, maxiter=100_000)
    m, n = size(A)
    n == m || error("Matrix must be square.")
    
    F = lu(A - ρ*I)
    
    q = randn(n)
    s = scale!(q)
    lam = 1/s + ρ

    qold = similar(q)

    k = 0
    done = false
    while !done && k < maxiter
        k += 1
        copy!(qold, q)        # qold = q
        ldiv!(q, F, qold)     # q = F\qold
        s = scale!(q)         # q = q/s
        lam = 1/s + ρ
        done = norm(A*q - lam*q)/norm(q) <= tol
    end

    if done
        println("Converged after $k iterations.")
    else
        println("Failed to converge.")
    end

    return lam, q
end
```


```julia
n = 1000
k = 5
density = (k - 1)/n    # density = (k*n - n)/n^2

A = triu(sprand(n, n, density), 1)
A = A + A' + I

ρ = 2.0

# Expect nnz(A) ≈ k*n
@show nnz(A)

if n <= 1000
    λ = eigvals(Matrix(A))
    abseig = abs.(λ .- ρ) |> sort
    r = abseig[1]/abseig[2]
    @show r
end

println()
@time lam, q = inverse_power_method(A, ρ=ρ)
@show ρ, lam
@show norm(A*q - lam*q)/norm(q);
```

---

## Rayleigh Quotient Iteration

Suppose $q \in \mathbb{C}^n$ approximates an eigenvector of $A$. If $A q = \rho q$, then $\rho$ is an eigenvalue of $A$. Otherwise, we want to find the value of $\rho$ that minimizes

$$
\| A q - \rho q \|_2.
$$

The _normal equations_ for this least squares problem is

$$
(q^* q) \rho = q^* A q,
$$

where $q^*$ is the **conjugate transpose** of $q$.

For example, if

$$
q = \begin{bmatrix} 1 + 3 i \\ 4 - 2 i \end{bmatrix},
$$

then

$$
q^* = \begin{bmatrix} 1 - 3 i & 4 + 2 i \end{bmatrix}.
$$

Note that $q^* q = \|q\|_2^2$.

The solution of the normal equations is

$$
\rho = \frac{q^* A q}{q^* q}
$$

and is called the **Rayleigh quotient**.

The **Rayleigh quotient iteration** uses

$$
\rho_j = \frac{q_j^* A q_j}{q_j^* q_j}
$$

as the _shift_ in each iteration of the inverse power method.

Since the shift changes each iteration, we need to compute an $LU$ decomposition _each iteration_. This can be very expensive since each iteration will now cost $\frac{2}{3}n^3 + O(n^2)$ flops.

To make the Raleigh quotient iteration practical, we can first compute a "simple" matrix $H$ that is _similar_ to $A$, such as an **upper Hessenberg** matrix

$$
H =
\begin{bmatrix}
* & * & * & * & * \\
* & * & * & * & * \\
  & * & * & * & * \\
  &   & * & * & * \\
  &   &   & * & * \\
\end{bmatrix}
$$

or a **tridiagonal** matrix

$$
H =
\begin{bmatrix}
* & * &   &   &   \\
* & * & * &   &   \\
  & * & * & * &   \\
  &   & * & * & * \\
  &   &   & * & * \\
\end{bmatrix}.
$$

Computing $LU$ decomposition of an upper Hessenberg matrix only needs $O(n^2)$ flops, and the same for a tridiagonal matrix only needs $O(n)$ flops. We will return to this topic in Section 5.5.

---

## `rayleigh`


```julia
using SuiteSparse

function rayleigh(A; ρ0=0.0, tol=sqrt(eps())/2, maxiter=100)
    m, n = size(A)
    n == m || error("Matrix must be square.")
    
    q = randn(n)
    normalize!(q)
    ρ = ρ0
    lam = dot(q, A*q)
    
    qold = similar(q)
    F = SuiteSparse.UMFPACK.UmfpackLU(
        Ptr{Nothing}(), Ptr{Nothing}(), 0, 0, 
        Int[], Int[], Float64[], 0)

    k = 0
    done = false
    while !done && k < maxiter
        k += 1
        copy!(qold, q)       # qold = q
        if k == 1
            F = lu(A - ρ*I)  # Creates symbolic factorization F
        else
            lu!(F, A - ρ*I)  # Overwrites F with new factorization
        end
        ldiv!(q, F, qold)    # q = (A - ρ*I)\qold
        normalize!(q)
        lam = dot(q, A*q)
        if k > 1
            ρ = lam
        end
        done = norm(A*q - lam*q) <= tol
    end

    if done
        println("Converged after $k iterations.")
    else
        println("Failed to converge.")
    end

    return lam, q
end
```


```julia
n = 4000
k = 5
density = (k - 1)/n    # density = (k*n - n)/n^2

A = triu(sprand(n, n, density), 1)
A = A + A' + I

ρ = 2.0

# Expect nnz(A) ≈ k*n
@show nnz(A)

if n <= 1000
    λ = eigvals(Matrix(A))
    abseig = abs.(λ .- ρ) |> sort
    r = abseig[1]/abseig[2]
    @show r
end

println()
println("Inverse power method:")
@time lam, q = inverse_power_method(A, ρ=ρ)
@show ρ, lam
@show norm(A*q - lam*q)/norm(q);

println()
println("Rayleigh quotient method:")
@time lam, q = rayleigh(A, ρ0=ρ)
@show ρ, lam
@show norm(A*q - lam*q);
```

---

## Quadratic convergence of the Raleigh Quotient Iteration

> ### Theorem: (Raleigh Quotient Approximates Eigenvalue)
>
> Let $A \in \mathbb{C}^{n \times n}$. Let $v$ be an eigenvector of $A$ with eigenvalue $\lambda$, and $\|v\|_2 = 1$.
>
> Let $q \in \mathbb{C}^n$ with $\|q\|_2 = 1$ and
>
> $$ \rho = q^* A q $$
>
> be the Raleigh quotient of $q$. Then
>
> $$ |\lambda - \rho| \le 2 \|A\|_2 \|v - q\|_2. $$

Therefore, if $\|v - q\|_2 = O(\varepsilon)$, then $|\lambda - \rho| = O(\varepsilon)$.

Let $q_0 \in \mathbb{C}^n$ such that $\|q_0\|_2 = 1$, and let $q_j$, for $j=1,2,\ldots$, be defined by

$$
\rho_j = q_j^* A q_j,
\qquad
(A - \rho_j I) \hat{q}_{j+1} = q_j,
\qquad
q_{j+1} = \hat{q}_{j+1}/\|\hat{q}_{j+1}\|_2.
$$

Then $\|q_j\|_2 = 1$, for all $j$.

1. Suppose that $q_j \to v_i$ as $j \to \infty$. Then $\|v_i\|_2 = 1$ and $\rho_j \to \lambda_i$ as $j \to \infty$.

2. Let $\lambda_k$ be the closest eigenvalue to $\lambda_i$.

3. Suppose that $\rho_j \approx \lambda_i$.

Then

$$
\begin{align}
\|v_i - q_{j+1}\|_2
&\approx \left| \frac{(\lambda_k - \rho_j)^{-1}}{(\lambda_i - \rho_j)^{-1}} \right| \|v_i - q_j\|_2 \\
&= \left| \frac{\lambda_i - \rho_j}{\lambda_k - \rho_j} \right| \|v_i - q_j\|_2 \\
&\le \frac{2 \|A\|_2 \|v_i - q_j\|_2}{|\lambda_k - \rho_j|} \|v_i - q_j\|_2 \\
&\approx \frac{2 \|A\|_2}{|\lambda_k - \lambda_i|} \|v_i - q_j\|_2^2. \\
\end{align}
$$

Thus, we obtain the estimate

$$ \|v_i - q_{j+1}\|_2 \approx C \|v_i - q_j\|_2^2, $$

where $C = 2 \|A\|_2 / |\lambda_k - \lambda_i|$. This indicates that the Rayleigh quotient iteration typically converges _quadratically_ when it does converge.

Moreover, if $A$ is a symmetric matrix, then $\|v - q\|_2 = O(\varepsilon)$ implies that $|\lambda - \rho| = O(\varepsilon^2)$, which indicates _cubic_ convergence:

$$ \|v_i - q_{j+1}\|_2 \approx C \|v_i - q_j\|_2^3. $$

---

a2e7dec56a1e9d28831c2f9702f53adf0292aa00
ipynb
Jupyter Notebook
Section 5.3 - The Power Method and Some Simple Extensions.ipynb
math434/fall2021math434
6317ce76de1eb7dbfdc3ea37a21dc5e1e3228316
MIT
2021-08-31T21:01:22.000Z
2021-08-31T21:01:22.000Z
Section 5.3 - The Power Method and Some Simple Extensions.ipynb
math434/fall2021math434
6317ce76de1eb7dbfdc3ea37a21dc5e1e3228316
MIT
Section 5.3 - The Power Method and Some Simple Extensions.ipynb
math434/fall2021math434
6317ce76de1eb7dbfdc3ea37a21dc5e1e3228316
MIT
2021-11-16T19:28:56.000Z
2021-11-16T19:28:56.000Z
Qwen/Qwen-72B
1. YES
2. YES