Re: [eigen] Numerical Stability Question

[Eric Chu <eytchu@xxxxxxxxx>]

Hi Benoit,

The linhuber() function is basically an iterative solver (uses Newton's method), so it works as follows:

- pick initial starting point (x0)
- compute the gradient and Hessian at x0
- line search to find step length, t

- x := x + t*x0

And iterate.

I used the SVD to compute the initial starting point. The gradient and Hessian are computed using QR factorizations.

I'll dig in my code some more. Mark's suggestion (about uninitialized memory) may be a good place to start looking.

Thanks!
Eric

On Mon, Jun 29, 2009 at 3:21 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Hi,

2009/6/30 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.

Thanks :)

> In my test code, I've got the following loop
>
> for (int i = 0; i < NUM_RUNS; i ++) {
>   lh.linhuber(A, b, x, 1);

It's hard to say anything without knowing what's inside linhuber(). At
least, what decomposition(s) is it using ?

Some of our decompositions aren't top quality in 2.0. That's
especially the case of the SVD. We're working on this for 2.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

I'm not aware of anything involving random numbers in any of our
decompositions algorithms, so I'm as surprised as you are. Again it's
hard to say anything more without looking at linhuber().

Cheers,
Benoit