[eigen] Numerical Stability Question

[Eric Chu <eytchu@xxxxxxxxx>]

Hi Eigen Folks,

First off, thanks for a (ridiculously) awesome library. I've managed to solve some robust least-squares problems in less than 2ms. Kudos to the developers.

Second, I've implemented a robust least-squares solver (solves Ax=b by minimizing a Huber loss function) by "translating" some pre-existing Matlab code, but am running into some numerical issues (it's not a big deal, since my application can tolerate small mistakes). I'd like to know if any of you have ideas about what's happening.



I'm running a 64-bit, quad core machine, and the following compiler flags (I'm using gcc 4.3.3):

-O3 -msse4 -m64 -Wall -DEIGEN_NO_DEBUG

In my test code, I've got the following loop

for (int i = 0; i < NUM_RUNS; i ++) {
  lh.linhuber(A, b, x, 1);

  cout << endl;
  cout << "Run: " << i+1 << endl;
  cout << x << endl;
}

Basically, it uses the same data to solve the least squares problem a couple of times. Running the same "type" of code in Matlab gives the same answer each time. However, when I run my program, I get the following output (assuming I run 5 times).


Run: 1
56.71
-20.04
-3606

Run: 2
56.82
-19.98
-3616

Run: 3
56.69
-20.05
-3604

Run: 4
56.67
-20.05
-3603

Run: 5
56.68
-20.05
-3604

The numbers aren't terribly "wrong," but I was wondering why each loop would give me a slightly different numerical answer. Is there any behavior in Eigen that *might* give rise to such a phenomenon?


Also, I checked the singular values of my matrix A (it's an 800x3 matrix). They came out to approx. 3000, 100, 0.5. If there's no other guess, my thought would be that the solver I've written is playing around in the "numerical" nullspace there.


Any feedback would be appreciated.

Once again, great work. Keep it up,
Eric.