When we spoke about linear systems we discovered how reducing matricies to alternate forms (LU decomposition) simplified computaiton. The same is true for the Eigenvalue problem.
Householder’s algorithm¶
The Householder algorithm is an iterative method that reduces a matrix to a simpler form which preserves the eigenvalues. It has properties that also prevent accumulation of roundoff error.
These matricies are then used in other methods.
The Householder algorithm reflects vectors in order to zero out elements of a matrix. Define the Householder reflection:
where
Apply to zero out elements below the diagonal. It can be used in a few ways:
Form a Hessenberg matrix (one subdiagonal) useful for input to a QR decomposition
Used in QR decomposition to find the factors.
For symmetric matricies, form a tridiagonal matrix for iterative Eigenvalue solvers
# Example of the Householder algorithm on a 4x4 matrix
import numpy as np
def householder_reflection(A):
"""Perform Householder reflections to reduce A to upper triangular form."""
m, n = A.shape
R = A.copy()
Q = np.eye(m)
for i in range(min(m, n)):
x = R[i:, i]
e = np.zeros_like(x)
e[0] = np.linalg.norm(x)
u = x - e
v = u / np.linalg.norm(u)
Q_i = np.eye(m)
Q_i[i:, i:] -= 2.0 * np.outer(v, v)
R = Q_i @ R
Q = Q @ Q_i
return Q, R
# Example usage
A = np.random.rand(4, 4)
Q, R = householder_reflection(A)
print("Matrix A:")
print(A)
print("Orthogonal matrix Q:")
print(Q)
print("Upper triangular matrix R:")
print(R)Hessenberg Matrix (R):
[[ 0.75627235 0.98090482 0.28772365 0.9381365 ]
[-1.09215513 -0.3857163 -0.80066735 -1.05307442]
[ 0. 0.34742699 -0.52132642 0.81655068]
[-0. 0. 0.02013476 0.31739987]]
Householder Matrix 1:
[[ 1. 0. 0. 0. ]
[ 0. -0.47738597 -0.82734878 -0.29596729]
[ 0. -0.82734878 0.53667761 -0.16574421]
[ 0. -0.29596729 -0.16574421 0.94070836]]
Householder Matrix 2:
[[ 1. 0. 0. 0. ]
[ 0. 1. 0. 0. ]
[ 0. 0. -0.93826854 0.34590772]
[ 0. 0. 0.34590772 0.93826854]]
QR decomposition¶
QR decomposition is another matrix decomposition technique that finds
where
is a real square matrix
is orthogonal;
is upper triangular.
If is invertible, the decomposition is unique.
# prompt: Build the QR decomposition for the Householder reflections above
import numpy as np
def qr_decomposition(A):
"""
Performs QR decomposition of a matrix A using Householder reflections.
Args:
A: A NumPy array representing the input matrix.
Returns:
A tuple (Q, R) where Q is an orthogonal matrix and R is an upper triangular matrix.
"""
m, n = A.shape
Q = np.eye(m)
R = A.copy()
for j in range(n):
x = R[j:, j]
e1 = np.zeros_like(x)
e1[0] = 1
norm_x = np.linalg.norm(x)
u = x + np.sign(x[0]) * norm_x * e1
u = u / np.linalg.norm(u)
H = np.eye(m)
H[j:, j:] -= 2 * np.outer(u, u)
R = H @ R
Q = Q @ H
return Q, R
# Example usage (using the previously defined 'A'):
A = np.random.rand(4, 4) # Example matrix
Q, R = qr_decomposition(A)
print("Q:\n", Q)
print("\nR:\n", np.round(R,14))
# Verify the decomposition
print("\nQ @ R:\n", np.dot(Q, R)) # Should be close to the original matrix A
print("\nA:\n", A)
print("\n From numpy: \n")
print("Q:\n",np.linalg.qr(A)[0])
print("R:\n",np.linalg.qr(A)[1])Q:
[[-0.33081541 0.14663562 0.9322079 0.00689863]
[-0.6638009 -0.31647359 -0.1809512 -0.65304632]
[-0.65049107 0.0125749 -0.23815611 0.72109979]
[-0.16367925 0.9371149 -0.20378117 -0.23129635]]
R:
[[-1.04217786 -1.5190878 -0.85009162 -1.25888391]
[ 0. 0.76285964 0.29873252 0.72804249]
[ 0. -0. 0.53931419 0.43086392]
[-0. -0. 0. -0.15868068]]
Q @ R:
[[0.3447685 0.61440006 0.82778119 0.92377524]
[0.6917986 0.76694692 0.37216107 0.63090254]
[0.67792739 0.99774594 0.42829257 0.61101032]
[0.17058289 0.96353029 0.30918698 0.83721295]]
A:
[[0.3447685 0.61440006 0.82778119 0.92377524]
[0.6917986 0.76694692 0.37216107 0.63090254]
[0.67792739 0.99774594 0.42829257 0.61101032]
[0.17058289 0.96353029 0.30918698 0.83721295]]
From numpy:
Q:
[[-0.33081541 0.14663562 0.9322079 -0.00689863]
[-0.6638009 -0.31647359 -0.1809512 0.65304632]
[-0.65049107 0.0125749 -0.23815611 -0.72109979]
[-0.16367925 0.9371149 -0.20378117 0.23129635]]
R:
[[-1.04217786 -1.5190878 -0.85009162 -1.25888391]
[ 0. 0.76285964 0.29873252 0.72804249]
[ 0. 0. 0.53931419 0.43086392]
[ 0. 0. 0. 0.15868068]]
The QR eigenvalue algorithm¶
The QR eigenvalue algorithm uses QR decomposition to find
and then reverses the matrix product to create a new iteration:
which reduces off diagonal components until a threshold is reached at which point the diagonal reveals the eigenvalues.
import numpy as np
def qr_eigenvalues(A, tol=1e-10, max_iter=1000):
n = A.shape[0]
Ak = A.copy()
for i in range(max_iter):
# Perform QR decomposition
Q, R = np.linalg.qr(Ak)
# Update Ak
Ak = R @ Q
print('\n Iteration ',i,'\n', np.round(Ak,10))
# Check for convergence (off-diagonal elements close to zero)
off_diagonal = np.sqrt(np.sum(np.tril(Ak, -1)**2))
if off_diagonal < tol:
break
# Eigenvalues are the diagonal elements
eigenvalues = np.diag(Ak)
return eigenvalues
# Example Usage
A = np.random.rand(3, 3) # Example matrix
eigenvalues = qr_eigenvalues(A)
print("Eigenvalues:", eigenvalues)
Iteration 0
[[ 1.5298972 -0.07394911 0.2732753 ]
[ 0.62182393 0.39698186 0.27647066]
[-0.01050214 -0.02906924 0.00747495]]
Iteration 1
[[ 1.55818229e+00 -5.71665984e-01 -3.26819374e-01]
[ 1.48377706e-01 3.55473939e-01 -1.79128054e-01]
[ 1.31967700e-04 1.49389580e-03 2.06977858e-02]]
Iteration 2
[[ 1.50740020e+00 -6.82647358e-01 3.39881389e-01]
[ 3.86707898e-02 4.05671771e-01 1.48795305e-01]
[-1.79430000e-06 -8.02663000e-05 2.12820459e-02]]
Iteration 3
[[ 1.49016550e+00 -7.10534758e-01 -3.43451134e-01]
[ 1.08490501e-02 4.22879224e-01 -1.40106898e-01]
[ 2.54000000e-08 4.08340000e-06 2.13092871e-02]]
Iteration 4
[[ 1.48501518e+00 -7.18270863e-01 3.44455197e-01]
[ 3.11623800e-03 4.28028220e-01 1.37606669e-01]
[-4.00000000e-10 -2.03900000e-07 2.13106095e-02]]
Iteration 5
[[ 1.48350981e+00 -7.20485908e-01 -3.44742857e-01]
[ 9.01357000e-04 4.29533523e-01 -1.36883738e-01]
[ 0.00000000e+00 1.01000000e-08 2.13106747e-02]]
Iteration 6
[[ 1.48307222e+00 -7.21126030e-01 3.44825945e-01]
[ 2.61243600e-04 4.29971117e-01 1.36674262e-01]
[-0.00000000e+00 -5.00000000e-10 2.13106779e-02]]
Iteration 7
[[ 1.48294520e+00 -7.21311512e-01 -3.44850014e-01]
[ 7.57619000e-05 4.30098130e-01 -1.36613519e-01]
[ 0.00000000e+00 0.00000000e+00 2.13106780e-02]]
Iteration 8
[[ 1.48290835e+00 -7.21365299e-01 3.44856993e-01]
[ 2.19751000e-05 4.30134980e-01 1.36595901e-01]
[-0.00000000e+00 -0.00000000e+00 2.13106781e-02]]
Iteration 9
[[ 1.48289766e+00 -7.21380899e-01 -3.44859017e-01]
[ 6.37430000e-06 4.30145669e-01 -1.36590791e-01]
[ 0.00000000e+00 0.00000000e+00 2.13106781e-02]]
Iteration 10
[[ 1.48289456e+00 -7.21385425e-01 3.44859604e-01]
[ 1.84900000e-06 4.30148770e-01 1.36589309e-01]
[-0.00000000e+00 -0.00000000e+00 2.13106781e-02]]
Iteration 11
[[ 1.48289366e+00 -7.21386737e-01 -3.44859774e-01]
[ 5.36400000e-07 4.30149670e-01 -1.36588879e-01]
[ 0.00000000e+00 0.00000000e+00 2.13106781e-02]]
Iteration 12
[[ 1.48289340e+00 -7.21387118e-01 3.44859824e-01]
[ 1.55600000e-07 4.30149931e-01 1.36588754e-01]
[-0.00000000e+00 -0.00000000e+00 2.13106781e-02]]
Iteration 13
[[ 1.48289333e+00 -7.21387229e-01 -3.44859838e-01]
[ 4.51000000e-08 4.30150006e-01 -1.36588718e-01]
[ 0.00000000e+00 0.00000000e+00 2.13106781e-02]]
Iteration 14
[[ 1.48289331e+00 -7.21387261e-01 3.44859842e-01]
[ 1.31000000e-08 4.30150028e-01 1.36588707e-01]
[-0.00000000e+00 -0.00000000e+00 2.13106781e-02]]
Iteration 15
[[ 1.48289330e+00 -7.21387270e-01 -3.44859844e-01]
[ 3.80000000e-09 4.30150035e-01 -1.36588704e-01]
[ 0.00000000e+00 0.00000000e+00 2.13106781e-02]]
Iteration 16
[[ 1.48289330e+00 -7.21387273e-01 3.44859844e-01]
[ 1.10000000e-09 4.30150036e-01 1.36588703e-01]
[-0.00000000e+00 -0.00000000e+00 2.13106781e-02]]
Iteration 17
[[ 1.48289330e+00 -7.21387273e-01 -3.44859844e-01]
[ 3.00000000e-10 4.30150037e-01 -1.36588703e-01]
[ 0.00000000e+00 0.00000000e+00 2.13106781e-02]]
Iteration 18
[[ 1.48289330e+00 -7.21387274e-01 3.44859844e-01]
[ 1.00000000e-10 4.30150037e-01 1.36588703e-01]
[-0.00000000e+00 -0.00000000e+00 2.13106781e-02]]
Eigenvalues: [1.4828933 0.43015004 0.02131068]