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