QR factorization

Suppose \(m \gg n\) and visualize the QR factorization, \[\mathbf{A} = \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 & \cdots & \mathbf{q}_m \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & \cdots & r_{1m} \\ 0 & r_{22} & \cdots & r_2m \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & r_{nn} \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix}.\]

\(\mathbf{Q}\) is an orthogonal matrix and \(\hat{\mathbf{Q}}\) is an ONC matrix.

opnorm(Q'*Q - I)

7.784676762789553e-16

opnorm(Q̂'*Q̂ - I)

4.980709636773836e-16

Least squares and QR

Suppose we have a thin QR factorization of \(\mathbf{A} = \hat{\mathbf{Q}} \hat{\mathbf{R}}\) and we are solving least-squares via the normal equations:

\[\begin{align*} \mathbf{A}^\top \mathbf{A} \mathbf{x} &= \mathbf{A}^\top \mathbf{b} \\ \hat{\mathbf{R}}^\top \hat{\mathbf{Q}}^\top \hat{\mathbf{Q}} \hat{\mathbf{R}} \mathbf{x} &= \hat{\mathbf{R}}^\top \hat{\mathbf{Q}}^\top \mathbf{b} \\ \hat{\mathbf{R}}^\top \hat{\mathbf{R}} \mathbf{x} &= \hat{\mathbf{R}}^\top \hat{\mathbf{Q}}^\top \mathbf{b} \end{align*}\]

This algorithm is essentially what is implemented in the Julia \ operator.

"""
    lsqrfact(A,b)
Solve a linear least-squares problem by QR factorization.  Returns the minimizer of ||b - Ax||.
"""
function lsqrfact(A,b)
    Q,R = qr(A)
    z = Q'*b
    x = FNC.backsub(R,z)
    return x
end

lsqrfact

Recall previously we saw that the normal equations can be unstable – this method doesn’t change that fact, but allows us to solve them up to the loss of accuracy predicted due to the instability.

t = range(0,3,length=400)
f = [ x->sin(x)^2, x->cos((1+1e-7)*x)^2, x->1. ]
A = [ f(t) for t in t, f in f ]
x = [1., 2, 1]
b = A*x;

observed_error = norm(lsqrfact(A,b) - x)/norm(x);
@show observed_error;
κ = cond(A)
@show error_bound = κ*eps();

observed_error = 6.4042143577854036e-9
error_bound = κ * eps() = 4.053030227715619e-9

Householder reflections

Note that \(\mathbf{P}\mathbf{x} = \mathbf{x} - 2 \mathbf{v}(\mathbf{v}^\top\mathbf{x})\). Visualizing this equation explains why these are called reflectors.

Now, we may choose \(\mathbf{v}\) so that this reflection \(\mathbf{P}\mathbf{z}\) is very sparse. In fact, we choose \(\mathbf{v}\) so that \[\mathbf{P}\mathbf{z} = \begin{bmatrix} \pm \|\mathbf{z}\| \\ 0 \\ \vdots \\ 0 \end{bmatrix} = \pm \|\mathbf{z}\| \mathbf{e}_1.\]

Factorization algorithm demo

A = rand(float(1:9),6,4)
m,n = size(A)

(6, 4)

z = A[:,1];
v = normalize(z - norm(z)*[1;zeros(m-1)])
P₁ = I - 2v*v'

6×6 Matrix{Float64}:
 0.263117   0.613941   0.175412    0.526235   0.350823    0.350823
 0.613941   0.48849   -0.146146   -0.438437  -0.292292   -0.292292
 0.175412  -0.146146   0.958244   -0.125268  -0.0835119  -0.0835119
 0.526235  -0.438437  -0.125268    0.624196  -0.250536   -0.250536
 0.350823  -0.292292  -0.0835119  -0.250536   0.832976   -0.167024
 0.350823  -0.292292  -0.0835119  -0.250536  -0.167024    0.832976

P₁*z

6-element Vector{Float64}:
 11.401754250991381
  0.0
  0.0
  4.440892098500626e-16
 -4.440892098500626e-16
 -4.440892098500626e-16

Now we recurse! We redefine \(\mathbf{A}\) and work on the submatrix that does not include its first row or column.

A = P₁*A

6×4 Matrix{Float64}:
 11.4018       12.1911     9.99846  11.8403
  0.0          -1.8245    -1.4982   -1.86589
 -1.11022e-16   6.05014    2.28623  -0.390253
  0.0           0.150431   5.85869   4.82924
  0.0           1.10029    2.57246   0.219493
  0.0           1.10029    2.57246   4.21949

z = A[2:m,2]
v = normalize(z - norm(z)*[1;zeros(m-2)])
P₂ = I - 2v*v'

5×5 Matrix{Float64}:
 -0.280271    0.929396    0.0231086    0.169021     0.169021
  0.929396    0.325317   -0.0167754   -0.122699    -0.122699
  0.0231086  -0.0167754   0.999583    -0.00305079  -0.00305079
  0.169021   -0.122699   -0.00305079   0.977686    -0.0223142
  0.169021   -0.122699   -0.00305079  -0.0223142    0.977686

A[2:m,:] = P₂*A[2:m,:]
A

6×4 Matrix{Float64}:
 11.4018       12.1911        9.99846  11.8403
 -1.03184e-16   6.50976       3.5497    1.02213
 -3.61174e-17   1.5525e-15   -1.37823  -2.48678
  1.86244e-18   3.75082e-17   5.76758   4.77711
  1.36223e-17   1.25457e-16   1.90604  -0.161783
  1.36223e-17   2.7015e-16    1.90604   3.83822

for j in 3:n
    z = A[j:m,j]
    v = normalize(z - norm(z)*[1;zeros(m-j)])
    P = I - 2v*v'
    A[j:m,:] = P*A[j:m,:]
end

R = triu(A)

6×4 Matrix{Float64}:
 11.4018  12.1911   9.99846  11.8403
  0.0      6.50976  3.5497    1.02213
  0.0      0.0      6.51386   5.83174
  0.0      0.0      0.0       3.12311
  0.0      0.0      0.0       0.0
  0.0      0.0      0.0       0.0

Now, note that we didn’t explicitly build \(\hat{\mathbf{Q}}\) here. \(\hat{\mathbf{Q}}\) is derived from the products of the Householder reflectors, but these matrices are getting smaller as they only need to act on a submatrix of \(\mathbf{A}\). To explicitly build \(\hat{\mathbf{Q}}\), we have to pad out the Householder reflectors \(\mathbf{P}_k\) to make them appropriately sized.

Define \(\mathbf{Q}_k = \begin{bmatrix} \mathbf{I}_{k-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{P}_k \end{bmatrix}\). First, note that \(\mathbf{Q}_k\) is orthogonal and that \[\mathbf{Q}_n \mathbf{Q}_{n-1} \cdots \mathbf{Q}_1 \mathbf{A} = \mathbf{R}.\] Thus, we have that \(\mathbf{Q}^\top = \mathbf{Q}_n \mathbf{Q}_{n-1} \cdots \mathbf{Q}_1\). Now, rather than explicitly building these \(\mathbf{Q}_k\) matrices explicitly and constructing their product, we can simply begin from \(\mathbf{I}\) and do the actions of the \(\mathbf{Q}_k\) matrices (the \(\mathbf{P}_k\) matrices) in the same way we did on \(\mathbf{A}\).

Q-less QR

Now, note that in many applications we don’t explicitly need the matrix \(\mathbf{Q}\).