C++ Eigen linear system solve, numerical issues?

Question

I'm getting unexpected and invalid results from a linear solve computation in Eigen. I have a vector x and matrix P.

In MATLAB notation, I'm doing this:

x'/P*x

(similar to x'*inv(P)*x, but no direct inversion)

which is quadratic form and should yield a positive result. The matrix P is positive definite and not ill conditioned. In MATLAB the result is positive though large.

In C++ Eigen, my slash operator is like this:

x/P is implemented as:

(P.transpose()).ColPivHouseholderQr().solve(x.transpose).transpose()

This gives the correct answer in general, but appears to be more fragile than MATLAB. In the case I'm looking at now, it gives a very large negative number for x'/P*x as though there was overflow and wraparound or somesuch.

Is there a way to defragilize this?

EDIT: Doing some experimentation I see that Eigen also fails to invert this matrix, where MATLAB has no problem. The matrix condition number is 3e7, which is not bad. Eigen gives the same (wrong) answer with the linear solve and a simple x.transpose() * P.inverse() * x. What's going on?

I presume that MATLAB uses a different algorithm than ColPivHouseholderQr() to solve this? — Cris Luengo
I'm not certain what MATLAB has under the hood. Quite possibly they look at the matrix and choose the best algorithm to use. The disturbing thing here is the inability to invert a matrix that I wouldn't consider that badly conditioned. Incidentally, if I change from ColPivHouseholderQr() to llt(), I can get it to work here, but still can't invert it using that method. — Mastiff
Are you sure that you use doubles (not floats) in C++/Eigen ? — Mikhail Genkin
Have you tried the simpler: x.transpose() * P.colPivHouseholderQr().solve(x)?? Since it's doing column-pivoting, working on P or P' might be quite different. Nonetheless, with 3e7 as condition number even an explicit computation of P.ibverse() should produce a correct result. Have you check the condition number on the C++ side? ( e.g., BDVSVD<MatrixXd> svd(P); cout << P.singularValues()(0)/P.singularValues()(P.cols()-1) << "\n";) — ggael

chtz chtz · Accepted Answer · 2018-08-22T22:24:40

The following is wrong (besides the () you missed in your question):

(P.transpose()).ColPivHouseholderQr().solve(x.transpose()).transpose();

because

x'*inv(P) = ((x'    *inv(P))')'
          = (inv(P)'* (x')'  )'
          = (inv(P) *  x     )'  % Note: P=P'

Your expression should actually raise an assertion in Eigen when built without -DNDEBUG, since x.transpose() has only one row but P has (usually) more.

To compute x'*inv(P)*x for symmetric positive definite P, I suggest using the Cholesky decomposition L*L'=P which gives you x'*inv(P)*x = (L\x)'*(L\x) which in Eigen is:

typedef Eigen::LLT<Eigen::MatrixXd> Chol;
Chol chol(P); // Can be reused if x changes but P stays the same
// Handle non positive definite covariance somehow:
if(chol.info()!=Eigen::Success) throw "decomposition failed!";
const Chol::Traits::MatrixL& L = chol.matrixL();
double quadform = (L.solve(x)).squaredNorm();

For the matrix Eigen failed to invert (which Matlab does invert), it would be interesting to see it.

C++ Eigen linear system solve, numerical issues?

1 Answers