How to solve a small least-squares problem

Question

This question is not very deep. Suppose I have a small rectangular matrix $A$, with number of rows and columns between $50$-$100$, respectively.

Given a right-hand side $b$, I want to solve the least-squares problem $Ax = b$, i.e., I want to find $x$ with norm-minimal residual $r = b - Ax$ and $x$ having minimum norm among these minimizers.

Since the system is very small, and $A$ has a good generalized condition number, I consider solving the normal equations $A^t A x = A^t b$ instead, say, by Cholesky decomposition.

At this stage, efficiency is less important than getting accurate solutions. Would you argue, from your experience, that this is a good choice, rather than, say, computing the SVD?

Answer 1

If need be, you can use the "seminormal equations"; the key here is that the $\mathbf R$ factor in the QR decomposition of $\mathbf A$ is the same (up to changes in signs of entries) as the Cholesky triangle of $\mathbf A^\top\mathbf A$, since

$\mathbf A^\top\mathbf A=(\mathbf Q\mathbf R)^\top(\mathbf Q\mathbf R)=\mathbf R^\top\mathbf Q^\top\mathbf Q\mathbf R=\mathbf R^\top\mathbf R$.

This is especially convenient if you are unable or unwilling to store the orthogonal matrix, either explicitly or as its decomposition in Householder vectors.

In short, you can solve the equations $\mathbf R^\top\mathbf R\mathbf x=\mathbf A^\top\mathbf b$ to obtain your least-squares solution.

Still, nothing beats the diagnostic ability afforded by the singular-value decomposition; if you can afford to use it, you should, unless you know well in advance that your least-squares problems will always be well-conditioned.

Answer 2

I'd guess a QR decomposition is better than solving the normal equations and faster than SVD.

There are some class notes that compare the three approaches.

Also: The QR decomposition for least-squares solution of overdetermined linear systems has been previously discussed here.

Answer 3

Sure, go with inverting $A^TA$. That matrix is symmetric and positive definite, so you can do that with CG as well. Alternatively, you can solve the original problem by computing a decomposition of the original, rectangular matrix, using an algorithm like QR.