Least squares and QR

We can us this to find the least square solution to \(\mathbf{A}\mathbf{x} = \mathbf{b}\) by substituting the QR factorization of \(\mathbf{A}\) into the normal equations:

\[ \begin{split} \mathbf{A}^T\mathbf{A} \mathbf{x} &= \mathbf{A}^T \mathbf{b}, \\ \hat{\mathbf{R}}^T \hat{\mathbf{Q}}^T \hat{\mathbf{Q}} \hat{\mathbf{R}} \mathbf{x} &= \hat{\mathbf{R}}^T \hat{\mathbf{Q}}^T \mathbf{b}, \\ \hat{\mathbf{R}}^T \hat{\mathbf{R}} \mathbf{x}& = \hat{\mathbf{R}}^T \hat{\mathbf{Q}}^T \mathbf{b}. \end{split} \]

As long as \(\mathbf{A}\) is not rank deficient, we then have \(\hat{\mathbf{R}} \mathbf{x}=\hat{\mathbf{Q}}^T \mathbf{b}\). Since \(R\) is upper triangular we can solve this using back subsitution!

function lsqrfact(A,b)
    Q,R = qr(A)
    c = Q'*b
    x = FNC.backsub(R,c)
    return x
end

Does this improve our previous solution to demo in 3.2?

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 = 2.1513528812733333e-9
# error_bound = κ * eps() = 4.053030227715619e-9