This small notebook shows a naive implementation of an algorithm to compute the eigen values and eigen vectors of a matrix, written in the Julia programming language.
For a real or complex valued square matrix $A \in \mathbb{M}_{n,n}(\mathbb{K})$ (on a field $\mathbb{K}=\mathbb{R}$ or $\mathbb{C}$), an eigenvalue $\lambda$ (of multiplicity $k\in\mathbb{N}^*$) and its associated eigenvector $v$ satisfy $$ (A - \lambda I)^k v = 0.$$
We will need at least two examples of matrix, to test our algorithms.
A = [ 1. 2 3; 4 5 6; 7 8 9 ]
B = [12 -51 4; 6 167 -68; -4 24 -41]
Based on the definition, it is not difficult to test if a pair $\lambda, v$ is an eigenvalue, eigenvector pair.
We can write a small function that will check if the eigen values and eigen vectors are correct.
function check_DV(A, D, V)
I = eye(size(A)[1])
for i = 1:size(A)[1]
check = (A - D[i] * I) * V[:,i]
println("For the ", i, "th eigen value and vector, (A - lambda I)^k * v = ", check, "\n\twhich has a L2 norm = ", norm(check))
assert(norm(check) <= 1e-14)
end
end
D_A, V_A = eig(A)
check_DV(A, D_A, V_A)
A second example, showing that $B$ has two eigen values that are $0$ (numerically found very close to $0$).
D_B, V_B = eig(B)
check_DV(B, D_B, V_B)
Let's compute the eigenpair with largest eigenvalue, through the power iteration method.
"Power iteration method on A, to find eigenpair with largest eigen value."
function poweriteration(A, maxiter=20, userandomstart=true)
n = size(A)[1]
if userandomstart
X = randn((n, 1))
else
X = ones((n, 1))
end
for i = 1:maxiter
next_X = A * X
X = next_X / norm(next_X)
end
# lambda = (X^* ⋅ A ⋅ X) / (X^* ⋅ X)
lambda = (transpose(conj(X)) * (A * X)) / (transpose(conj(X)) * X)
lambda, X
end
For instance, we obtain correctly the largest eigenpair for our example $A$. (upto a $-1$ sign in $v$)
poweriteration(A)
D_A, V_A = eig(A)
(D_A[1], -1 * V_A[:,1])
What is the typical dependency on the number of iterations?
poweriteration
println(poweriteration(A, 1))
println(poweriteration(A, 2))
println(poweriteration(A, 5))
println(poweriteration(A, 10))
The QR method can be used to compute eigenvalues, and then X method to compute an approximation of a correspond eigenvector for each eigenvalue.
First, we need to implement the QR decomposition, which can be based on the Gram-Schmidt process.
It requires a generic inner-product function:
"Inner product of two vectors v, w. Can be vertical (n,1) or horizontal (1,n) vectors."
function inner(v, w)
assert(size(v) == size(w))
nm = size(v)
if length(nm) == 1
return transpose(conj(v)) * w
elseif nm[1] == 1
return conj(v) * transpose(w)
end
end
It works for both $(1,n)$ and $(n,1)$ vectors:
inner([1 1 1], [2 3 4])
inner([1; 1; 1], [2; 3; 4])
Then a projection operator:
"projection(a, e): projection operator of vector a onto base vector e. Uses inner(e, a) and inner(e, e)."
function projection(a, e)
return (inner(e, a) / inner(e, e)) * e
end
Then the Gram-Schmidt process:
"gramschmidt(A): Gram-Schmidt orthogonalization operator, returns us, es (unormalized matrix, base matrix)."
function gramschmidt(A, verbose=false)
assert(size(A)[1] == size(A)[2])
n = size(A)[1]
if verbose
println("n = ", n)
end
us = zeros((n, n))
es = zeros((n, n))
for i = 1:n
if verbose
println("i = ", i)
end
us[:,i] = A[:,i]
for j = 1:i-1
if verbose
println("\tj = ", j)
println("\tus[:,i] = ", us[:,i])
end
us[:,i] -= projection(A[:,i], us[:,j])
end
if verbose
println("us[:,i] = ", us[:,i])
end
es[:,i] = us[:,i] / norm(us[:,i])
if verbose
println("es[:,i] = ", es[:,i])
end
end
return us, es
end
Example:
A
us, es = gramschmidt(A)
us
es
We can check that $a_1 = <e_1,a_1> e_1$, $a_2 = <e_1,a_2> e_1 + <e_2,a_2> e_2$ and so on:
inner(es[:,1], A[:,1]) * es[:,1]
inner(es[:,1], A[:,2]) * es[:,1] + inner(es[:,2], A[:,2]) * es[:,2]
It's not true for the last one, as $A$ was not full rank.
rank(A)
inner(es[:,1], A[:,3]) * es[:,1] + inner(es[:,2], A[:,3]) * es[:,2] + inner(es[:,3], A[:,3]) * es[:,3]
$E$ is an orthogonal matrix indeed:
[ norm(es[:,i]) for i = 1:size(es)[1] ]
(es' * es)[1:2,1:2]
Second example:
us, es = gramschmidt(B)
[ norm(es[:,i]) for i = 1:size(es)[1] ]
es * es'
Last example:
V = [3 2; 1 2]
us, es = gramschmidt(V)
es' * es
It is very easy from the Gram-Schmidt decomposition:
"QR decomposition, returns (Q, R): Q is orthogonal and R is upper-triangular, such that A = Q R."
function QR(A)
assert(size(A)[1] == size(A)[2])
n = size(A)[1]
us, es = gramschmidt(A)
Q = copy(es)
R = zeros((n, n))
for i = 1:n
for j = i:n
R[i,j] = inner(es[:,i], A[:,j])
end
end
return Q, R
end
First example
Q, R = QR(A)
A
Q
R
Q * R
Second example
Q2, R2 = QR(B)
B
Q2
R2
Q2 * R2
Due to numerical errors, $A \neq QR$ but almost:
B == Q2 * R2
assert(norm(B - (Q2 * R2)) <= 1e-13)
"Apply the QR method for maxiter steps. Should return a triangular matrix similar to A (same eigenvalues)."
function QR_method(A, maxiter=50)
Ak = A
for k = 1:maxiter
Qk, Rk = QR(Ak)
Ak = Rk * Qk
end
return Ak
end
It should produce matrix which are almost triangular!
Ak = QR_method(A)
"Truncate to zero values under the diagonal (have to be smaller than tolerance)"
function truncate_to_zero_below_diag(A, tolerance=1e-12)
assert(size(A)[1] == size(A)[2])
n = size(A)[1]
for j = 1:n
for i = j+1:n
assert(norm(A[i,j]) <= tolerance)
A[i,j] = 0
end
end
return A
end
truncate_to_zero_below_diag(Ak)
Ak = QR_method(B)
truncate_to_zero_below_diag(Ak)
Finally, we can extract the eigen values from the diagonal:
"Extract the diagonal from the QR method."
function QR_eigenvalues(A, maxiter=50)
return diag(QR_method(A, maxiter))
end
It does not work for matrices which are not full rank:
QR_eigenvalues(A)
eigvals(A)
But it works fine for full rank matrix:
QR_eigenvalues(B)
eigvals(B)
Now, we have a numerical algorithm that gives an approximation of the eigen values. For each eigen value $\lambda$, we can use the inverse iteration to find the corresponding eigen vector.
"Inverse iteration method to find the eigen vector corresponding to a given eigen value."
function inverseIteration(A, val, maxiter=20, userandomstart=true)
mu_I = val * eye(A)
inv_A_minus_muI = inv(A - mu_I)
n = size(A)[1]
if userandomstart
X = randn((n, 1))
else
X = ones((n, 1))
end
for i = 1:maxiter
next_X = inv_A_minus_muI * X
X = next_X / norm(next_X)
end
X
end
Now, putting it all together:
"Approximation of D, V: D contains the eigen values (vector) and V the eigen vectors (column based)."
function QR_eigen(A, maxiter=20, userandomstart=true)
n = size(A)[1]
D = QR_eigenvalues(A, maxiter)
V = zeros((n,n))
for i = 1:n
V[:,i] = inverseIteration(A, D[i], maxiter, userandomstart)
end
return D, V
end
D2, V2 = QR_eigen(B)
D2s, V2s = eig(B)
Up to a permutation, and up to multiplication by $-1$ of some eigen vectors, we got the good values!
It is well known that finding eigenpairs and finding roots of a polynomial are equivalent problems (thanks to companion matrices). The Abel-Ruffini theorem states that there is no exact algorithm to find roots of polynomial of degrees larger than $4$, and therefore there is no hope of finding an exact algorithm for the eigenpair problem.
Hence, we focussed only on one numerical (i.e., iterative) method.
It's sad to see that there is not more details about general algorithms for the eigenpair problem.
That's it for today, folks! See here for other notebooks I wrote.