★ Eigenvalues And Singular Values ★


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)

12.1 power Iteration methods

Definition

Let A be an $m \times m$ matrix. A dominant eigenvalue of A is an eigenvalue $\lambda$ whose magnitude is greater than all other eigenvalues of A. If it exists, an eigenvector associated to $\lambda$ is called a dominant eigenvector

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

Example

Find the dominant eigenvalue and eigenvector for $A = \begin{bmatrix} 1 & 3 \\ 2 & 2 \end{bmatrix} $


In [4]:
A = np.array([
    [1, 3],
    [2, 2]
])

x = np.random.rand(2)

power_iteration(A, x, 25)


Out[4]:
(4.0, array([0.70710678, 0.70710678]))

Theorem

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|$

Inverse Power Iteration

Lemma

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

Example

$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))


(-1.0000000000000004, array([-7.07106782e-01,  4.63492990e-09,  7.07106780e-01]))
(-0.12310562561766053, array([ 0.26095647, -0.92941026,  0.26095647]))
(8.123105625617661, array([0.6571923 , 0.36904818, 0.6571923 ]))

Rayleigh Quotient Iteration


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

Example


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))


(8.12310562561766, array([-0.6571923 , -0.36904818, -0.6571923 ]))

12.2 QR Algorithm

Normalized Simultaneous Iteration


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

Example

$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)


eigenvalues : 
[ 8.12310563 -1.         -0.12310563]

eigenvectors : 
[[-6.57192300e-01  7.07106781e-01  2.60956472e-01]
 [-3.69048184e-01 -2.33071126e-09 -9.29410263e-01]
 [-6.57192299e-01 -7.07106781e-01  2.60956475e-01]]

eigenvalues (np.linalg.eig) : 
[ 8.12310563 -1.         -0.12310563]

eigenvectors (np.linalg.eig) : 
[[ 6.57192300e-01  7.07106781e-01  2.60956474e-01]
 [ 3.69048184e-01 -8.76089486e-17 -9.29410263e-01]
 [ 6.57192300e-01 -7.07106781e-01  2.60956474e-01]]

Unshifted QR Algorithm


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

Theorem

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.

Example

$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)


eigenvalues : 
[ 8.12310563 -1.         -0.12310563]

eigenvectors : 
[[ 6.57192300e-01 -7.07106781e-01  2.60956472e-01]
 [ 3.69048184e-01  2.33071110e-09 -9.29410263e-01]
 [ 6.57192299e-01  7.07106781e-01  2.60956475e-01]]

eigenvalues (np.linalg.eig) : 
[ 8.12310563 -1.         -0.12310563]

eigenvectors (np.linalg.eig) : 
[[ 6.57192300e-01  7.07106781e-01  2.60956474e-01]
 [ 3.69048184e-01 -8.76089486e-17 -9.29410263e-01]
 [ 6.57192300e-01 -7.07106781e-01  2.60956474e-01]]

Real Schur form and the QR algorithm

Definition

A matrix $T$ has real Schur form if it is upper triangular, except possibly for $2 \times 2$ blocks on the main diagonal

Theorem

Let $A$ be a square matric with real entries. Then there exists an orthogonal matrix $Q$ and a matrix $T$ in real Schur form such that $A = Q^{T}TQ$


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

Example


In [14]:
A = np.array([
    [3,   2,  4],
    [2,   1,  2],
    [4,   2,  3]
])

print(shifted_qr(A))


[ 8.12310361 -0.99999799 -0.12310563]