Dense decomposition of very non-square matrices

Question

I have inherited code that solves Eigen::Matrix problems using the code shown on this page:

(A.transpose() * A).ldlt().solve(A.transpose() * b)

This works, but often my matrix size can be $[200 \times 5]$. The $A^T.A$ expression gives me a huge square matrix with 400 times as much data, so I suspect that the solver wastes a lot of time handling that extra data.

Is there a "better" solver that avoids creating a square matrix, or can the problem be approached by a "divide and conquer strategy, where perhaps I could solve 40 50x5 problems, then solve the 40x5 results?

UPDATE: For further information, my problem is I have a point cloud from after doing a disparity check on a stereo image, and according to the function name, it is trying to detect a road surface. We solve the matrix, then use the calculated coefficients to weed out noise in each group of points.

It also strikes me that trying to accurately solve for all 2000 points seems like a lot of effort; perhaps picking each 10th point would give a roughly similar answer?

Furthermore, the 5 row values are something like (IIRC, I'm away from my computer for a week!) $X$, $Y$, $Z$, $X^2$ and $Y^2$

Finally, note that one operation in itself is not that inefficient, it's that I have over 300 groups of points for a grand total of over 100,000 points per frame.

Answer 1

As others have mentioned in comments, solving $A^TAx=A^Tb$ (the so-called least-norm solution) can be performed without explicitly forming $A^TA$ with the QR decomposition. It works like this:

$$A=QR$$ where $Q$ has the interesting properties that $Q^TQ=I$ and $R$ is square and triangular. Using this information, the least squares system $A^TAx=A^Tb$ is equivalent to: $$R^TQ^TQRx=R^TQ^Tb$$ Harnessing the properties of $Q$ and $R$, $$x=R^{-1}R^{-T}R^TQ^Tb=R^{-1}Q^Tb$$ Because $R$ is triangular, $R^{-1}$ is very easy to calculate. With Eigen, you don't have to worry about any of that.

A.colPivHouseholderQr().solve(b)

This is probably more accurate than forming the normal equations explicitly. If you need more speed, Intel MKL has a QR solver that will crush your problem. The API sucks compared to working with Eigen, however.