Implementation

This leads us to a way to solve our linear system, we just use our previous methods to solve the normal equations as a \(n\times n\) system. Since \(\mathbf{A}^T\mathbf{A}\) is symmetric and positive definite, we can use the Cholesky (\(\mathbf{N} = \mathbf{R}^T\mathbf{R}\)) factorization:

function lsnormal(A,b)
    N = A'*A;  z = A'*b;
    R = cholesky(N).U
    w = FNC.forwardsub(R',z)                   # solve R'z=c
    x = FNC.backsub(R,w)                       # solve Rx=z
    return x
end

This takes \(\sim (mn^2 + \frac{1}{3}n^3)\) flops