3x3 Householder QR Demo

This demo constructs a $3\times 3$ QR factorization using Householder reflectors.

import numpy as np
import numpy.linalg as la

n = 3

e1 = np.array([1,0,0])
e2 = np.array([0,1,0])
e3 = np.array([0,0,1])

A = np.random.randn(n, n)
A

array([[ 1.06427586,  1.03785412,  1.19643376],
       [-1.06405176,  0.48693156, -0.83335032],
       [-0.20227028, -1.10857722,  0.27986689]])

Householder reflector: $$I-2\frac{vv^T}{v^Tv}$$

Choose $v=a-\|a\|e_1$.

a = A[:, 0]
v = a-la.norm(a)*e1

H1 = np.eye(3) - 2*np.outer(v, v)/(v@v)

A1 = H1 @ A
A1

array([[  1.51848691e+00,   5.33870203e-01,   1.38523070e+00],
       [  3.40401588e-16,  -6.93719959e-01,  -3.91067578e-01],
       [  9.99443798e-17,  -1.33301245e+00,   3.63942360e-01]])

NB: Never build full Householder matrices in actual code! (Why? How?)

a = A1[:, 1].copy()
a[0] = 0
v = a-la.norm(a)*e2

H2 = np.eye(3) - 2*np.outer(v, v)/(v@v)

R = H2 @ A1
R

array([[  1.51848691e+00,   5.33870203e-01,   1.38523070e+00],
       [ -2.45801147e-16,   1.50272073e+00,  -1.42307423e-01],
       [ -2.55820086e-16,   1.61523465e-16,   5.14914060e-01]])

Q = np.dot(H2, H1).T
la.norm(np.dot(Q, R) - A)

8.5458331069102901e-16