Solving eigenvalue problems numerically

Cours ENPC — Pratique du calcul scientifique

Objective

Goal of this chapter

Study numerical methods to find the eigenpairs $(\mathbf{\boldsymbol{v}}_i, \mathbf{\boldsymbol{\lambda}}_i)$ satisfying \[ \mathsf A \mathbf{\boldsymbol{v}} = \lambda \mathbf{\boldsymbol{v}}, \] where $\mathsf A \in \mathbf C^{n \times n}$ and $b \in \mathbf C^n.$

$~$

Simplifying assumption throughout this chapter: the matrix $\mathsf A$ is diagonalizable.

Notation: eigenvalues $|\lambda_1| \geqslant\dotsc \geqslant|\lambda_n|$, associated normalized eigenvectors $\mathbf{\boldsymbol{v}}_1, \dotsc, \mathbf{\boldsymbol{v}}_n$. \[ \mathsf A \mathsf V = \mathsf V \mathsf D, \qquad \mathsf V = \begin{pmatrix} \mathbf{\boldsymbol{v}}_1 & \dots & \mathbf{\boldsymbol{v}}_n \end{pmatrix}, \qquad \mathsf D = {\rm diag}(\lambda_1, \dotsc, \lambda_n). \]

Finding all the eigenvalues of a matrix is expensive:
```
using LinearAlgebra
A = rand(2000, 2000)
@time eigen(A)
```
```
6.694395 seconds
```
${\color{green} \rightarrow}$ We first focus on methods for finding just one or a few eigenpairs.

Program of today

Simple vector iterations, for calculating just one eigenpair;
Subspace iterations: generalization to calculate multiple eigenpairs at the same time;
Projection methods: approximate large eigenvalue problem by a smaller one.

Remark: can we calculate the eigenvalues as roots of the characteristic polynomial?

Cost of calculating the characteristic polynomial scales as ${\color{red}n!}$

${\color{green} \rightarrow}$ not a viable approach…

Before we begin:

We present a motivating example;
We discuss sparse matrix formats

Motivating example: Random walk on a graph

Consider a random walk on the graph:

Assumption: jumps to adjacent nodes equi-probable.

The associated adjacency matrix $\mathsf A \in \{0, 1\}^{n \times n}$ is given by \[ a_{ij} = \begin{cases} 1 & \text{if there is an edge between $i$ and $j$,} \\ 0 & \text{otherwise.} \end{cases} \]

$~$

Probability of transition $i \rightarrow j$ given by \[ p_{ij} = \frac{a_{ij}}{n_e(i)}, \] where $n_e(i)$ is the number of edges connected to node $i$.

${\color{green} \rightarrow}$ note that $\mathsf A$ is symmetric but $\mathsf P := (p_{ij})$ is not.

$~$

Steady state probability distribution $\mathbf{\boldsymbol{v}}_1 \in \mathbf R_{\geqslant 0}^n$ satisfies \[ \mathsf P^T\mathbf{\boldsymbol{v}}_1 = \mathbf{\boldsymbol{v}}_1, \qquad \mathbf{\boldsymbol{1}}^T\mathbf{\boldsymbol{v}}_1 = 1. \] (it is possible to show that $1$ is the dominant eigenvalue)

${\color{green} \rightarrow}$ eigenvalue problem for matrix $\mathsf P^T$.

Constructing the adjacency matrix

function adj(m)
  x = vcat([collect(0:i) for i in 0:m-1]...)
  y = vcat([fill(m-i, i) for i in 1:m]...)
  row_lengths = 1:m
  row_indices = [0; cumsum(row_lengths)]
  A = zeros(Int, length(x), length(x))
  for i in 1:m-1
      for j in 1:row_lengths[i]
          r1 = row_indices[i]
          r2 = row_indices[i+1]
          A[r2 + j, r2 + j + 1] = 1  # Horizontal
          A[r1 + j, r2 + j]     = 1  # Vertical
          A[r1 + j, r2 + j + 1] = 1  # Diagonal
      end
  end
  x, y, A + A'
end

Here ... is the splat operator (?... for more info)

sum3(a, b, c) = a + b + c
numbers = [1, 2, 3]
sum3(numbers...)

Also works with keyword arguments (use ; and NamedTuple)

plotargs = (color=:blue, label="example")
plot(sin; plotargs...)

Exercise: write function transition(m) that builds $\mathsf P$.

Notice that $\mathsf A$ is a sparse matrix:

# Adjancency matrix
x, y, A = adj(10);
Plots.spy(A, xlabel=L"i", ylabel=L"j")
Plots.plot!(size=(600, 600))

This is the case in many applications:

Discretization of partial differential equations;
Graph applications.

Calculating the steady state

Steady state probability distribution satisfies \[ \mathsf P^T\mathbf{\boldsymbol{v}}_1 = \mathbf{\boldsymbol{v}}_1 \]

Illustration with units in ‰, rounded to integers:

using LinearAlgebra
v₁ = eigvecs(P', sortby=real)[:, end] |> real
v₁ = v₁ / sum(v₁)  # Normalization
plotgraph(v₁)

Strategy for finding $\mathbf{\boldsymbol{v}}_1$: evolve probability vector until equilibrium: \[ \mathbf{\boldsymbol{x}}^{(k+1)} = \mathsf P^T\mathbf{\boldsymbol{x}}^{(k)} \qquad \Leftrightarrow \qquad x^{(k+1)}_i = \sum_{j=1}^{n_n} p_{ji} x^{(k)}_j \]

$x^{(k)}_i$ is the probability of being at node $i$ at iteration $k$;
This is an example of the power iteration;
Below we take $\mathbf{\boldsymbol{x}}^{(0)} = \mathbf{\boldsymbol{e}}_1$.