1

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?

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
Mastiff
  • 2,101
  • 2
  • 23
  • 38
  • Possibly relevant: [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken). – Jesper Juhl Aug 22 '18 at 16:49
  • 1
    I presume that MATLAB uses a different algorithm than `ColPivHouseholderQr()` to solve this? – Cris Luengo Aug 22 '18 at 17:46
  • 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 Aug 22 '18 at 18:08
  • 2
    Are you sure that you use doubles (not floats) in C++/Eigen ? – Mikhail Genkin Aug 22 '18 at 19:22
  • Everything is MatrixXd and VectorXd. I'm assuming that carries through to all the internal computations. – Mastiff Aug 22 '18 at 20:03
  • 1
    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 svd(P); cout << P.singularValues()(0)/P.singularValues()(P.cols()-1) << "\n";`) – ggael Aug 22 '18 at 20:46
  • In this particular case, P is symmetric anyway (a covariance matrix), so the transpose shouldn't matter, unless it's somehow lossy, which I sure hope it isn't. I'll ask Eigen what it thinks the condition number is... – Mastiff Aug 22 '18 at 21:02

1 Answers1

3

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.

chtz
  • 17,329
  • 4
  • 26
  • 56
  • Thanks for the reply. I think just a misunderstanding regarding the wrongness. My implementation of x/P is as described (I overloaded forward slash). For finding x'/P, replace x with x.transpose() in the messy expression. I'll try your suggestion for reduced fragility. – Mastiff Aug 23 '18 at 00:24