Matrix decomposition methods#
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:
\[ H = \rm{I} - 2 \frac{u u^T}{u^T u} \]
where \( u = x- \pm \|x\| e_1\)
Apply \(A \leftarrow H A\) 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
# prompt: Can you give me an example of the householder algorithm on a 4x4 matrix?
import numpy as np
def householder_reflection(x):
"""Computes the Householder reflection matrix for a given vector x."""
x = np.asarray(x, dtype=float)
norm_x = np.linalg.norm(x)
if norm_x == 0:
return np.identity(len(x))
e1 = np.zeros_like(x)
e1[0] = 1
u = x + np.sign(x[0]) * norm_x * e1 # Choose the sign to avoid cancellation
u = u / np.linalg.norm(u)
H = np.identity(len(x)) - 2 * np.outer(u, u)
return H
def householder_algorithm(A):
"""Applies the Householder algorithm to a 4x4 matrix A."""
n = A.shape[0]
H = []
R = A.copy()
for k in range(n - 2): # Iterate until a Hessenberg form is achieved.
x = R[k+1:, k]
H_k = householder_reflection(x)
# Embed H_k in a larger matrix.
H_k_full = np.identity(n)
H_k_full[k+1:, k+1:] = H_k
# Apply the Householder reflection
R = H_k_full @ R
H.append(H_k_full)
return R, H # Return Hessenberg matrix R and accumulated reflections H.
# Example 4x4 matrix
A = np.random.rand(4,4)
# Apply the Householder algorithm
R, H = householder_algorithm(A)
# Print the Hessenberg matrix
print("Hessenberg Matrix (R):\n", np.round(R,14))
#Print the householder matrices
for i, h in enumerate(H):
print(f"\nHouseholder Matrix {i+1}: \n", h)
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 $\(A = QR\)$
where
\(A\) is a real square matrix
\(Q\) is orthogonal; \(Q^T = Q^{-1}\)
\(R\) is upper triangular.
If \(A\) 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 $\( A = QR\)\( and then reverses the matrix product to create a new iteration: \)\( A^{k+1} = R^k Q^k\)$
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]