Computing Eigenvalues

Computing Eigenvalues#

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

Power Iteration#

Suppose a matrix \(A\) has a dominant eigenvalue \(\lambda_1\). This means that \(\lambda_1\) has multiplicity 1 and \(| \lambda_1 | > | \lambda |\) for all other eignvalues \(\lambda\). Choose a random vector \(\mathbf{x}_0\) and compute the (normalized) power iteration algorithm:

\[ \mathbf{x}_{k+1} = \frac{A \mathbf{x}_k}{\| A \mathbf{x}_k \|_{\infty}} \]

The result (usually) converges to an eigenvector for \(\lambda_1\). Let us implement the algorithm with pauses at each iteration so that we can see the vector \(\mathbf{x}_k\) at each step.

import time

A = np.array([[1,1,0],[1,1,1],[0,1,1]])
x0 = np.array([[1],[0],[0]])
iterations = 5
print('\nIteration k=0\n\nxk =\n\n{}'.format(x0))

xk = x0
for k in range(1,iterations+1):
    xk = A@xk
    xk = xk/np.max(np.abs(xk))
    print('\nIteration k={0}\n\nxk =\n\n{1}'.format(k,xk))
    time.sleep(1)

Iteration k=0

xk =

[[1]
 [0]
 [0]]

Iteration k=1

xk =

[[1.]
 [1.]
 [0.]]

Iteration k=2

xk =

[[1. ]
 [1. ]
 [0.5]]

Iteration k=3

xk =

[[0.8]
 [1. ]
 [0.6]]

Iteration k=4

xk =

[[0.75      ]
 [1.        ]
 [0.66666667]]

Iteration k=5

xk =

[[0.72413793]
 [1.        ]
 [0.68965517]]

Now let us implement the method without pauses and increase the number of iterations.

A = np.array([[1,1,0],[1,1,1],[0,1,1]])
x0 = np.array([[1],[0],[0]])
iterations = 20
xk = x0

for k in range(0,iterations):
    xk = A@xk
    xk = xk/np.max(np.abs(xk))

print(xk)

[[0.70710681]
 [1.        ]
 [0.70710675]]

Compute the Rayleigh quotient to approximate the corresponding eigenvalue:

xk.T @ A @ xk / (xk.T @ xk)

array([[2.41421356]])

Compare to the function scipy.linalg.eig.

evals,evecs = la.eig(A)

print(evals)

[-0.41421356+0.j  1.        +0.j  2.41421356+0.j]

print(evecs)

[[ 5.00000000e-01  7.07106781e-01  5.00000000e-01]
 [-7.07106781e-01  4.91354351e-16  7.07106781e-01]
 [ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]]

Indeed, power iteration produces an approximation of the eigenvector for eigenvalue \(\lambda_1 = 1 + \sqrt{2}\).

Inverse Iteration#

Suppose matrix \(A\) is invertible. Then the smallest eigenvalue \(\lambda_n\) (in absolute value) of \(A\) corresponds to the dominant eigenvalue \(1/\lambda_n\) of \(A^{-1}\). Choose a random vector \(\mathbf{x}_0\) and compute the (normalized) inverse iteration algorithm:

\[\begin{split} A \tilde{\mathbf{x}}_{k+1} = \mathbf{x}_k \\ \mathbf{x}_{k+1} = \tilde{\mathbf{x}}_k / \| \tilde{\mathbf{x}}_{k+1} \|_{\infty} \end{split}\]

The result (usually) converges to an eigenvector for \(A\) for \(\lambda_n\) (which is also an eigenvector for \(A^{-1}\) for \(1/\lambda_n\)). Let us implement the algorithm.

A = np.array([[1,1,0],[1,1,1],[0,1,1]])
x0 = np.array([[1],[0],[0]])
iterations = 20
xk = x0
LU,P = la.lu_factor(A)

for k in range(0,iterations):
    xk = la.lu_solve((LU,P),xk)
    xk = xk/np.max(np.abs(xk))

print(xk/la.norm(xk))

[[ 0.50000002]
 [-0.70710678]
 [ 0.49999998]]

This is the eigenvector for \(A\) for \(\lambda = 1 - \sqrt{2}\).

Random matrices#

Let us implement the power iteration method for \(B^TB\) for a random matrix \(B\) of size \(N\). Since \(B^TB\) is symmetric we know the eigenvalues are real. Compute the Rayleigh quotient of the result and compare with the function scipy.linalg.eig.

N = 100
B = np.random.randn(N,N)
A = B.T@B
x0 = np.random.randn(N,1)
iterations = 150
xk = x0

for k in range(0,iterations):
    xk = A@xk
    xk = xk/np.max(np.abs(xk))

lambda_1 = xk.T@A@xk/(xk.T@xk)
print("Power method: ", lambda_1[0,0])

evals,evecs = la.eig(A)
print("SciPy: ", np.max(np.abs(evals.real)))

Power method:  375.75262486612826
SciPy:  376.1225114616083

The result is pretty good for 150 iterations!

QR Iteration#

Start with \(A_0 = A\), for each step \(k\), compute the QR decomposition \(A_k = Q_k R_k\), define \(A_{k+1} = R_k Q_k\) and repeat. The matrices \(A_k\) are similar therefore they have the same eigenvalues. The result is an upper (block) triangular matrix with eigenvalues on the diagonal.

A = np.array([[1,1,0],[1,1,1],[0,1,1]])
iterations = 5
print('\nIteration k=0\n\nA =\n\n{}'.format(A))

Ak = A
for k in range(1,iterations+1):
    Q,R = la.qr(Ak)
    Ak = R@Q
    print('\nIteration k={0}\n\nAk =\n\n{1}'.format(k,Ak))
    time.sleep(1)

Iteration k=0

A =

[[1 1 0]
 [1 1 1]
 [0 1 1]]

Iteration k=1

Ak =

[[ 2.00000000e+00 -7.07106781e-01  2.02930727e-17]
 [-7.07106781e-01  1.00000000e+00  7.07106781e-01]
 [ 0.00000000e+00  7.07106781e-01  0.00000000e+00]]

Iteration k=2

Ak =

[[ 2.33333333e+00 -3.33333333e-01 -3.20655176e-16]
 [-3.33333333e-01  1.00000000e+00  3.33333333e-01]
 [ 0.00000000e+00  3.33333333e-01 -3.33333333e-01]]

Iteration k=3

Ak =

[[ 2.40000000e+00 -1.41421356e-01  2.07336995e-16]
 [-1.41421356e-01  1.00000000e+00  1.41421356e-01]
 [ 0.00000000e+00  1.41421356e-01 -4.00000000e-01]]

Iteration k=4

Ak =

[[ 2.41176471e+00 -5.88235294e-02 -2.50505098e-16]
 [-5.88235294e-02  1.00000000e+00  5.88235294e-02]
 [ 0.00000000e+00  5.88235294e-02 -4.11764706e-01]]

Iteration k=5

Ak =

[[ 2.41379310e+00 -2.43829925e-02  2.24955710e-16]
 [-2.43829925e-02  1.00000000e+00  2.43829925e-02]
 [ 0.00000000e+00  2.43829925e-02 -4.13793103e-01]]

Let us implement the method without pauses and increase the number of iterations.

A = np.array([[1,1,0],[1,1,1],[0,1,1]])
iterations = 50
Ak = A

for k in range(0,iterations):
    Q,R = la.qr(Ak)
    Ak = R@Q

print(Ak)

[[ 2.41421356e+00 -3.77213497e-16 -2.29396765e-16]
 [-1.45293347e-19  1.00000000e+00 -1.23531192e-16]
 [ 0.00000000e+00  1.45293347e-19 -4.14213562e-01]]

We can see the eigenvalues on the diagonal of the matrix \(A_k\).