2
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • 6
    $\begingroup$ Is this the solution of a least squares problem min_x II Ax-b II ? I suspect that it is. In that case, there is a large number of solvers that can solve that straight-up, without multiplication with A transpose. This is called Linear Least-Squares, and wiki has methods for this, e.g. using QR decomposition. (Q,R) = qr(A), then x = R^(-1) * transpose(Q) * b. Of course, R^(-1) should be done with a triangular matrix solve, since R is upper triangular. $\endgroup$ – OscarB Apr 26 '18 at 10:05
  • 2
    $\begingroup$ hm, if A is 2000x5, then you end up with a tiny 5x5 matrix. $\endgroup$ – ggael Apr 26 '18 at 10:44
  • 3
    $\begingroup$ hmm, even if you have $A\in \mathbb R^{5\times 2000}$ - resulting in a $2000 \times 2000$ matrix after $A^T A$, in general, I don't see how it will be mostly zeroes. There is something special about $A$, I suspect. $\endgroup$ – Anton Menshov♦ Apr 26 '18 at 11:32
  • 2
    $\begingroup$ ...and the resulting matrix won't be invertible anyway. $\endgroup$ – ggael Apr 26 '18 at 13:06
  • 1
    $\begingroup$ @ggael Oops, that's completely correct! Therefore my divide and conquer idea makes zero sense... That's why my log files showed much less than 1 ms per iteration. Thinning the data out might be the only solution. (An upgrade of Eigen to a CUDA-compatible version is not easy due to lots of third-party dependencies.) $\endgroup$ – Ken Y-N Apr 28 '18 at 4:53
3
$\begingroup$

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.

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.