In [1]:
#import modules
import numpy as np
from matplotlib import pyplot as plt
In [2]:
# Help function
def is_zero_vector(v):
"""
Check whether vector v is a zero vector
Arguments:
- v : vector
Return:
- True if v is a nonzero vector. Otherwise, false
Exceptions:
- TypeError:
- v is not a vector
"""
if v.ndim != 1 and not np.isscalar(v):
raise TypeError('v is not a vector')
return not np.any(v)
Eigenvalue equation $x\lambda = Ax$, where $x$ is an approximate eigenvector and $\lambda$ is unknown.
The least squares solution of $x^{T}x\lambda = Ax$, or $$ \lambda = \frac{x^{T}Ax}{x^{T}x} $$ known as the Rayleigh quotient
In [3]:
def power_iteration(A, x, k=10):
"""
Compute the dominant eigenvalue and eigenvector of a matrix A
Arguments:
- A : A matrix
- x : initial and nonzero vector
- k : number of steps (default:10)
Returns:
- eigval : dominant eigenvalue
- eigvec : dominant eigenvector
Exceptions:
- ValueError:
- x is a zero vector
"""
if is_zero_vector(x):
raise ValueError('x is a zero vector')
for _ in range(k):
eigvec = x / np.linalg.norm(x)
x = np.matmul(A, eigvec)
eigval = np.matmul(np.matmul(eigvec.T, A), eigvec)
eigvec = x / np.linalg.norm(x)
return eigval, eigvec
In [4]:
A = np.array([
[1, 3],
[2, 2]
])
x = np.random.rand(2)
power_iteration(A, x, 25)
Out[4]:
Let A be an $m \times m$ matrix with real eigenvalues $\lambda_1,\cdots,\lambda_m$ satisfying $|\lambda_1| > |\lambda_2| \ge |\lambda_3| \ge \cdots \ge |\lambda_m|$. Assume that the eigenvectors of A span $R^m$. For almost every initial vector, Power Iteration converges linearly to an eigenvector associated to $\lambda_1$ with convergence rate constant $S = |\lambda_2/\lambda_1|$
Let the eigenvalues of the $m \times m$ matrix A be denoted by $\lambda_1,\lambda_2,\cdots,\lambda_m$.
(a) The eigenvalues of the inverse matrix $A^{-1}$ are $\lambda_1^{-1},\lambda_2^{-1},\cdots,\lambda_m^{-1}$, assuming that the inverse exists. The eigenvactors are the same as those of A.
(b) The eigenvalues of the shifted matrix $A - sI$ are $\lambda_1 - s,\lambda_2 - s,\cdots,\lambda_m - s$ and the eigenvectors are the same as those of A.
In [5]:
def inverse_power_iteration(A, x, s, k=10):
"""
Compute eigenvalue and eigenvector of a matrix A nearest to input s
Arguments:
- A : A matrix
- x : nonzero initial vector
- s : shift
- k : number of steps (default:10)
Returns:
- eigval : dominant eigenvalue
- eigvec : dominant eigenvector
"""
As = A - s * np.eye(A.shape[0])
for _ in range(k):
eigvec = x / np.linalg.norm(x)
x = np.linalg.solve(As, eigvec)
eigval = np.matmul(eigvec.T, x)
u = x / np.linalg.norm(x)
eigval = 1 / eigval + s
return eigval, eigvec
$A = \begin{bmatrix} 3 & 2 & 4 \\ 2 & 1 & 2 \\ 4 & 2 & 3 \\ \end{bmatrix} $
where
$
\left\{\begin{matrix}\begin{align*}
\lambda_{1} &= -1 \\
v_{1} &= (-x_1, 0, x_1)
\end{align*}\end{matrix}\right.
$
$
\left\{\begin{matrix}\begin{align*}
\lambda_{2} &= -\sqrt{17} + 4 \approx -0.1231056256176605 \\
v_{2} &= (x_2, \frac{-\sqrt{17}-3}{2} \times x_2, x_2)
\end{align*}\end{matrix}\right.
$
$
\left\{\begin{matrix}\begin{align*}
\lambda_{3} &= \sqrt{17} + 4 \approx 8.1231056256176605 \\
v_{3} &= (x_3, \frac{\sqrt{17}-3}{2} \times x_3, x_3)
\end{align*}\end{matrix}\right.
$
In [6]:
A = np.array([
[3, 2, 4],
[2, 1, 2],
[4, 2, 3]
])
x = np.random.rand(3)
print(inverse_power_iteration(A, x, -1.1))
print(inverse_power_iteration(A, x, 0))
print(inverse_power_iteration(A, x, 8))
In [7]:
def rayleigh_quotient_iteration(A, x, k=10):
"""
Compute the dominant eigenvalue and eigenvector of a matrix A
Arguments:
- A : A matrix
- x : nonzero initial vector
- k : number of steps (default:10)
Returns:
- eigval : dominant eigenvalue
- eigvec : dominant eigenvector
"""
for _ in range(k):
eigvec = x / np.linalg.norm(x)
eigval = np.matmul(np.matmul(eigvec.T, A), eigvec)
x = np.linalg.solve(A - eigval * np.eye(A.shape[0]), eigvec)
eigvec = x / np.linalg.norm(x)
eigval = np.matmul(np.matmul(eigvec.T, A), eigvec)
return eigval, eigvec
In [8]:
A = np.array([
[3, 2, 4],
[2, 1, 2],
[4, 2, 3]
])
x = np.random.rand(3)
print(rayleigh_quotient_iteration(A, x))
In [9]:
def normalized_simultaneous_iteration(A, k=10):
"""
Compute the eigenvalue and eigenvector of a symmetric matrix A
Arguments:
- A : A matrix
- k : number of steps (default:10)
Returns:
- eigval : dominant eigenvalue
- eigvec : dominant eigenvector
"""
m, n = A.shape
Q = np.eye(m)
for _ in range(k):
Q, R = np.linalg.qr(np.matmul(A, Q))
eigval = np.diag(np.matmul(np.matmul(Q.T, A), Q))
eigvec = Q
return eigval, eigvec
$A = \begin{bmatrix} 3 & 2 & 4 \\ 2 & 1 & 2 \\ 4 & 2 & 3 \\ \end{bmatrix} $
where
$
\left\{\begin{matrix}\begin{align*}
\lambda_{1} &= -1 \\
v_{1} &= (-x_1, 0, x_1)
\end{align*}\end{matrix}\right.
$
$
\left\{\begin{matrix}\begin{align*}
\lambda_{2} &= -\sqrt{17} + 4 \approx -0.1231056256176605 \\
v_{2} &= (x_2, \frac{-\sqrt{17}-3}{2} \times x_2, x_2)
\end{align*}\end{matrix}\right.
$
$
\left\{\begin{matrix}\begin{align*}
\lambda_{3} &= \sqrt{17} + 4 \approx 8.1231056256176605 \\
v_{3} &= (x_3, \frac{\sqrt{17}-3}{2} \times x_3, x_3)
\end{align*}\end{matrix}\right.
$
In [10]:
A = np.array([
[3, 2, 4],
[2, 1, 2],
[4, 2, 3]
])
eigval, eigvec = normalized_simultaneous_iteration(A)
print('eigenvalues : ')
print(eigval)
print()
print('eigenvectors : ')
print(eigvec)
print()
eigval, eigvec = np.linalg.eig(A)
print('eigenvalues (np.linalg.eig) : ')
print(eigval)
print()
print('eigenvectors (np.linalg.eig) : ')
print(eigvec)
In [11]:
def unshifted_qr(A, k=10):
"""
Compute the eigenvalue and eigenvector of a symmetric matrix A
Arguments:
- A : A matrix
- k : number of steps (default:10)
Returns:
- eigval : dominant eigenvalue
- eigvec : dominant eigenvector
"""
m, n = A.shape
Q = np.eye(m)
Qbar = Q
R = A
for _ in range(k):
Q, R = np.linalg.qr(np.matmul(R, Q))
Qbar = np.matmul(Qbar, Q)
eigval = np.diag(np.matmul(R, Q))
eigvec = Qbar
return eigval, eigvec
Assume that $A$ is a symmetric $m \times m$ matrix with eigenvalues $\lambda_i$ satisfying $|\lambda_1| > |\lambda_2| > \cdots > |\lambda_m| $. The unshifted QR algorithm converages linearly to the eigenvectors and eigenvalues of $A$. As $j \rightarrow \infty$, $A_j$ converages to a diagonal matrix containing the eigenvalues on the main diagonal and $\bar{Q_j} = Q_1 \cdots Q_j$ converges to an orthogonal matrix whose columns are the eigenvectors.
$A = \begin{bmatrix} 3 & 2 & 4 \\ 2 & 1 & 2 \\ 4 & 2 & 3 \\ \end{bmatrix} $
where
$
\left\{\begin{matrix}\begin{align*}
\lambda_{1} &= -1 \\
v_{1} &= (-x_1, 0, x_1)
\end{align*}\end{matrix}\right.
$
$
\left\{\begin{matrix}\begin{align*}
\lambda_{2} &= -\sqrt{17} + 4 \approx -0.1231056256176605 \\
v_{2} &= (x_2, \frac{-\sqrt{17}-3}{2} \times x_2, x_2)
\end{align*}\end{matrix}\right.
$
$
\left\{\begin{matrix}\begin{align*}
\lambda_{3} &= \sqrt{17} + 4 \approx 8.1231056256176605 \\
v_{3} &= (x_3, \frac{\sqrt{17}-3}{2} \times x_3, x_3)
\end{align*}\end{matrix}\right.
$
In [12]:
A = np.array([
[3, 2, 4],
[2, 1, 2],
[4, 2, 3]
])
eigval, eigvec = unshifted_qr(A)
print('eigenvalues : ')
print(eigval)
print()
print('eigenvectors : ')
print(eigvec)
print()
eigval, eigvec = np.linalg.eig(A)
print('eigenvalues (np.linalg.eig) : ')
print(eigval)
print()
print('eigenvectors (np.linalg.eig) : ')
print(eigvec)
In [13]:
def shifted_qr(A, tol=1e-14, max_count=1000):
"""
"""
m = A.shape[0] # row size
eigval = np.zeros(m)
n = m
while n > 1:
count = 0
while np.max(A[n-1, 0:n-1]) > tol and count < max_count:
count += 1
shift = A[n-1, n-1]
Q, R = np.linalg.qr(A - shift * np.eye(n))
A = np.matmul(R, Q) + shift * np.eye(n)
if count < max_count:
eigval[n-1] = A[n-1, n-1]
n -= 1
A = A[0:n, 0:n]
else:
disc = (A[n-2, n-2] - A[n-1, n-1])^2 + 4 * A[n-1, n-2] * A[n-2, n-1]
eigval[n-1] = (A[n-2, n-2] + A[n-1,n-1] + np.sqrt(disc)) / 2
eigval[n-2] = (A[n-2, n-2] + A[n-1,n-1] - np.sqrt(disc)) / 2
n -= 2
A = A[0:n, 0:n]
if n > 0:
eigval[0] = A[0, 0]
return eigval
In [14]:
A = np.array([
[3, 2, 4],
[2, 1, 2],
[4, 2, 3]
])
print(shifted_qr(A))