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