×

### Dustin Lennon

##### Applied Scientistdlennon.org

(206) 291-8893

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/20210427_householder
April 2021

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

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