0

I'm using Eigen v3.2.7.

I have a medium-sized rectangular matrix X (170x17) and row vector Y (170x1) and I'm trying to solve them using Eigen. Octave solves this problem fine using X\Y, but Eigen is returning incorrect values for these matrices (but not smaller ones) - however I suspect that it's how I'm using Eigen, rather than Eigen itself.

auto X = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>{170, 17};
auto Y = Eigen::Matrix<T, Eigen::Dynamic, 1>{170};

// Assign their values...

const auto theta = X.colPivHouseholderQr().solve(Y).eval(); // Wrong!

According to the Eigen documentation, the ColPivHouseholderQR solver is for general matrices and pretty robust, but to make sure I've also tried the FullPivHouseholderQR. The results were identical.

Is there some special magic that Octave's mldivide does that I need to implement manually for Eigen?

Update

This spreadsheet has the two input matrices, plus Octave's and my result matrices.

Replacing auto doesn't make a difference, nor would I expect it to because construction cannot be a lazy operation, and I have to call .eval() on the solve result because the next thing I do with the result matrix is get at the raw data (using .data()) on tail and head operations. The expression template versions of the result of those block operations do not have a .data() member, so I have to force evaluation beforehand - in other words theta is the concrete type already, not an expression template.

The result for (X*theta-Y).norm()/Y.norm() is:

2.5365e-007

And the result for (X.transpose()*X*theta-X.transpose()*Y).norm() / (X.transpose()*Y).norm() is:

2.80096e-007

As I'm currently using single precision float for my basic numerical type, that's pretty much zero for both.

cmannett85
  • 21,725
  • 8
  • 76
  • 119
  • it would be better if you provide data to reproduce the error. –  Feb 03 '16 at 23:02
  • 1
    if the question is not how to fix it , but why this is so, then look on this answer: http://stackoverflow.com/a/18553768/4651282 –  Feb 03 '16 at 23:07
  • You want me to paste a 170x17 floating point matrix!? – cmannett85 Feb 03 '16 at 23:07
  • you can't upload them somewhere and leave a link? –  Feb 03 '16 at 23:09
  • 1
    Just a guess (I'm still waking up): Don't use `auto` with Eigen. See e.g. [here](http://stackoverflow.com/questions/31099246/wrong-results-using-auto-with-eigen) and [here](http://stackoverflow.com/a/26925302/2899559). – Avi Ginsburg Feb 04 '16 at 04:40
  • 1
    How did you check that the result is incorrect? First, as Avi said, replace auto by `Eigen::VectorXd`, then print `(X*theta-Y).norm()/Y.norm()` as well as `(X.transpose()*X*theta-X.transpose()*Y).norm() / (X.transpose()*Y).norm()`. – ggael Feb 04 '16 at 08:20
  • @Atomic_alarm - I've updated my question with your response. Thanks for your time! – cmannett85 Feb 04 '16 at 09:54
  • @ggael - I've updated my question with your response. Thanks for your time! – cmannett85 Feb 04 '16 at 09:55

1 Answers1

1

According to your verifications, the solution you get is perfectly fine. If you want more accuracy, then use double floating point numbers. Note that MatLab/Octave use double precision by default.

Moreover, it might also likely be that your problem is not full rank, in which case your problem admit an infinite number of solution. ColPivHouseholderQR picks one, somehow arbitrarily. On the other hand, mldivide will pick the minimal norm one that you can also obtain with Eigen::BDCSVD (Eigen 3.3), or the slower Eigen::JacobiSVD.

ggael
  • 28,425
  • 2
  • 65
  • 71
  • So are you saying that although the result I get from Eigen is not the same as Octave - it is equivalent? I tried using `JacobiSVD`, but it gave me a different result to Octave (regardless of the 'thin'/'full' U/V state). – cmannett85 Feb 04 '16 at 10:23
  • Yes or no. The values you get confirm that Eigen find one solution of the problem you asked. But maybe you are not asking for the exact same problem. If you want to dig deeper, then 1) make sure you use double for the comparison, 2) print the norm of the solution for both Eigen and Octave, 3) print the singular values `(X.jacobiSvd().singularValues().transpose())`. – ggael Feb 04 '16 at 10:29
  • Thanks for your help, after plotting the frequency and phase response of what I'm getting from Eigen against Octave (my result is a polynomial ratio representing an IIR), they're pretty close. Close enough that it is probably a precision problem, I'll investigate further. – cmannett85 Feb 04 '16 at 10:33
  • 2
    Even using single precision, when I tried to reproduce your issues I got the "exact" same result as Octave (both with and without auto). See [link](https://ideone.com/jYEjpB), just edit it to copy paste to your favorite editor. What compiler etc. are you using? – Avi Ginsburg Feb 04 '16 at 10:39
  • I also get the same as Avi. Maybe you are miss-parsing the `-9.11E-05` and similar entries if using a text file to exchange data? – ggael Feb 04 '16 at 10:50
  • I'm converting a MATLAB script to C++, and the `X` and `Y` matrices are generated from other input data. I've compared the C++ matrices to Octave's and they are exactly the same, so I know the code up until point I solve them is correct. Interestingly, if I change to double precision then my frequency response graphs then match Octave's accurately - implying that the solve results are equivalent, but not the same. What version of Eigen are you using? – cmannett85 Feb 04 '16 at 11:18
  • @cmannett85 Same as you, 3.2.7. I compiled using VS2013. I can compile with g++ (Windows / Ubuntu) if that matches your toolchain. Did you copy and compile the MCVE? – Avi Ginsburg Feb 04 '16 at 13:14
  • @AviGinsburg I get the same results with v3.2.5 on MinGW64 and 3.2.7 on Linux (both using g++ 5.2.1). I haven't tried using MSVC. – cmannett85 Feb 04 '16 at 13:38
  • Interesting. Using double precision I get the same results between VS2013 and MinGW g++ 4.8.1. Using single precision, I get different results with the two compilers. – Avi Ginsburg Feb 04 '16 at 14:01