Dustin Lennon
Applied Scientist
dlennon.org
numerical linear algebra
householder
Numerical Linear Algebra: Householder a'la Golub and Van Loan
Exploring properties of Householder transformations
Dustin Lennon
April 2021
https://dlennon.org/archive/20210427_householder
April 2021
Abstract
Abstract
Abstract
This is a document for exploring properties of Householder transformations.
import numpy as np
from pydlennon.contexts.seeded_rng import SeededRng
np.set_printoptions( precision=3, linewidth = 132, suppress = True )
m = 8
n = 6
I_m = np.identity(m)
I_n = np.identity(n)
with SeededRng(111) as s:
data = s.negative_binomial(10, 0.25, size=m*n)
A = np.array(data, dtype=np.float64).reshape((m,n))
data = s.negative_binomial(10, 0.25, size=n*n)
B = np.array(data, dtype=np.float64).reshape((n,n))
data = s.negative_binomial(10, 0.25, size=n*n)
Z,_ = np.linalg.qr(
np.array(data, dtype=np.float64).reshape((n,n))
)
def house(x):
"""
Golub and Van Loan, 3rd Ed.
Algorithm 5.1.1
"""
n = len(x)
sigma = x[1:] @ x[1:]
v = np.concatenate([[1], x[1:]])
if sigma == 0:
beta = 0
else:
mu = np.sqrt(x[0]**2 + sigma)
if x[0] <= 0:
v[0] = x[0] - mu
else:
v[0] = -sigma / (x[0] + mu)
beta = 2 * v[0]**2 / (sigma + v[0]**2)
v = v/v[0]
return(v, beta)
x = A[:,0]
v,beta = house( x )
P = (I_m - beta * np.outer(v,v))
P @ x
array([102.113, 0. , 0. , 0. , -0. , -0. , -0. , 0. ])
R = A.copy()
V = np.identity(m)
Beta = np.zeros(n)
for j in range(n):
x = R[j:,j]
v,beta = house( x )
s = np.s_[j:,j:]
# Atmp = P @ Atmp
w = beta * R[s].T @ v
R[s] = R[s] - np.outer(v,w)
# Factored form representation
Beta[j] = beta
V[j,j] = 1.0
V[j+1:,j] = v[1:]
# Reconstruct Q from factored form; backward accumulation
Q = I_m.copy()
for j in range(n)[::-1]:
v = V[j:,j].copy()
beta = Beta[j]
s = np.s_[j:,j:]
P = np.identity(m-j) - beta * np.outer(v,v)
Q[s] = P @ Q[s] # Q_1 ... Q_k ... Q_n
(Q.T @ A, Q @ R)
(array([[102.113, 75.671, 65.506, 74.163, 47.115, 94.631],
[ 0. , 65.642, 41.91 , 42.107, 28.85 , 30.943],
[ -0. , 0. , 54.419, 2.924, 12.635, 8.219],
[ 0. , 0. , -0. , 28.023, 10.352, 8.562],
[ 0. , 0. , -0. , 0. , 14.215, -18.268],
[ 0. , -0. , -0. , -0. , -0. , 18.734],
[ -0. , -0. , 0. , -0. , 0. , -0. ],
[ 0. , -0. , 0. , 0. , 0. , 0. ]]),
array([[26., 26., 26., 22., 17., 14.],
[50., 21., 15., 23., 15., 46.],
[11., 40., 32., 39., 21., 30.],
[27., 60., 29., 26., 18., 40.],
[42., 22., 63., 20., 25., 34.],
[29., 41., 20., 44., 34., 25.],
[60., 23., 5., 35., 10., 47.],
[14., 32., 43., 36., 19., 44.]]))
Hessenberg Reduction
Hessenberg Reduction
Hessenberg Reduction
B0 = Z.T @ B @ Z
# B0 = B
Btmp = B0.copy()
U = I_n.copy()
for k in range(n-2):
v,beta = house(Btmp[k+1:n, k])
vp = np.pad(v, (k+1,0))
P = I_n - beta * np.outer(vp,vp)
Btmp = P.T @ Btmp @ P
U = U @ P
H = U.T @ B0 @ U
(B0, H)
(array([[164.533, 37.299, 36.604, -17.722, 30.584, -27.13 ],
[ 39.477, 16.303, 3.115, -10.634, 3.408, -10.418],
[ 24.111, 29.105, 14.486, 9.866, 9.773, 12.609],
[ -2.569, 4.244, 1.124, 8.255, -5.707, -10.352],
[ 12.692, 12.145, -3.231, -6.987, 2.092, 20.604],
[-13.23 , -18.403, 11.126, 11.049, -16.502, -0.67 ]]),
array([[164.533, 63.174, 3.406, -20.94 , 13.669, 9.417],
[ 49.825, 32.936, -18.525, 6.785, 4.039, 1.218],
[ 0. , 17.621, -0.586, 9.554, 2.336, -14.628],
[ 0. , 0. , 25.482, -1.78 , -9.527, -1.698],
[ 0. , 0. , 0. , 24.062, 10.328, -2.967],
[ 0. , 0. , 0. , -0. , 9.138, -0.432]]))
e1 = np.zeros(n)
e1[0] = 1
(Z @ e1, Z @ U @ e1, U @ e1)
(array([-0.548, -0.43 , -0.252, -0.548, -0.296, -0.252]),
array([-0.548, -0.43 , -0.252, -0.548, -0.296, -0.252]),
array([1., 0., 0., 0., 0., 0.]))