Behavior Change for FindInterpolatingPolynomial due to Change in Eigen FullPivLu

[weig...@light.co]

     Hi, we found that for some of our optimization problems Ceres (using LINE_SEARCH) have failed to converge if we combine it with Eigen of versions after 3.2.8.

     We looked into this problem and found that it is because the behavior of Eigen fullPivLu function has changed after 3.2.8 (You can check the change log of 3.2.8: http://eigen.tuxfamily.org/index.php?title=ChangeLog#Eigen_3.2.8 or more specifically: https://bitbucket.org/eigen/eigen/commits/84835f5f75c7a736e03a7e98a095f3ac7439defc?at=default ) especially for matrix that is numerically rank deficient, 

     This caused very different behavior for FindInterpolatingPolynomial (line 373 in internal/ceres/polynomial.cc) when the matrix lhs becomes numerically rank deficient, we found that for many optimization problems they converge better with the previous behavior of Eigen fullPivLu (version 3.2.7 and before). So we temporarily fix this problem by changing line 373 

 return lhs.fullPivLu().solve(rhs);

to the following:

 Eigen::FullPivLU<Matrix> lu_lhs( lhs );

  return lu_lhs.setThreshold( 0.0 ).solve( rhs );

The fix above so far works well for us, however, we believe you guys have a better fix for this since probably the solution given by fullPivLU is not the solution you are looking for when the matrix becomes rank deficient.

For your reference, some people encountered similar issues while using Eigen ( https://forum.kde.org/viewtopic.php?f=74&t=129439 ), we tried the suggestions Gael provided, however, they either converge much slower or failed to converge for our testing problem.

   We are looking forward to hearing back from you. Thanks a lot and happy holiday!

[Sameer Agarwal]

weiguang,

Thank you for the detailed bug report. I am surprised by these changes in Eigen's numerical behaviour.  I have filed 

https://github.com/ceres-solver/ceres-solver/issues/248

to track this and will take a look at this after the holidays.

Sameer

--
You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/fb95eb10-2360-4b87-9cca-32ae5a7e5b6e%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Sameer Agarwal]

One thing that would help here are some example problem or just matrices which are giving you different results.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/CAK0oyEr%3Dozt%2By%2Bv-VdGJaueEsvutXoE68VdJNg_4EaFeLhX8LA%40mail.gmail.com.

[Sameer Agarwal]

weiguang,

Could you try 

https://ceres-solver-review.googlesource.com/#/c/9280/

and let me know if this fixes the problem for you.

Sameer

[weig...@light.co]

    Thank you so much, Sameer, I'll try it this week and get back to you.

    Thanks again.

[weig...@light.co]

Hi, Sameer, thanks a lot for your attention and help. I tried this patch, it still failed the convergence after the first iteration.

1)   For our testing problem, this is printed information using the patch you provided (the optimization stops after the first iteration):

   0: f: 5.221393e+05 d: 0.00e+00 g: 5.24e+05 h: 0.00e+00 s: 0.00e+00 e:  0 it: 9.41e-03 tt: 9.49e-03

WARNING: Logging before InitGoogleLogging() is written to STDERR

W1229 11:02:50.069532 20645 line_search.cc:772] Line search failed: Wolfe zoom phase failed to find a point satisfying strong Wolfe conditions within specified max_num_iterations: 20, (num iterations taken for bracketing: 1).

Here is the initial cost and final cost:

summary.initial_cost = 522139.339699653 summary.final_cost = 522139.339699653

If we switched to the fix we came up with (lu_lhs.setThreshold( 0.0 ).solve( rhs );), here are the printed information

   0: f: 5.221393e+05 d: 0.00e+00 g: 5.24e+05 h: 0.00e+00 s: 0.00e+00 e:  0 it: 9.85e-03 tt: 9.93e-03

   1: f: 5.151351e+05 d: 7.00e+03 g: 3.86e+05 h: 6.91e-02 s: 1.21e-07 e:  2 it: 7.56e-03 tt: 1.76e-02

...

   760: f: 2.995621e-01 d: 6.87e-11 g: 6.53e-02 h: 8.57e-05 s: 1.00e+00 e:  1 it: 4.66e-03 tt: 3.39e+00

   761: f: 2.995621e-01 d: 3.56e-11 g: 4.51e-02 h: 2.51e-04 s: 1.00e+00 e:  1 it: 3.33e-03 tt: 3.39e+00

Here is the initial cost and final cost:

summary.initial_cost = 522139.339699653 summary.final_cost = 0.299562082

2) Here is one unit test I created to check the behavior of fullPivLu on a rank deficient matrix (this matrix is actually from the first iteration of our testing problem in FindInterpolatingPolynomial):

    ceres::Matrix lhs = ceres::Matrix::Zero( 4, 4 );

    ceres::Vector rhs = ceres::Vector::Zero( 4 );

    lhs << 0, 0, 0, 1,

        0, 0, 1, 0,

        6.94062e-18, 3.63858e-12, 1.90751e-06, 1,

        1.09157e-11, 3.81501e-06, 1, 0;

    rhs << 522139, -3.25765e+11, 1.81769e+06, -4.06082e+10;

    auto result1 = lhs.fullPivLu().solve( rhs );

    Eigen::FullPivLU<ceres::Matrix> lu_lhs( lhs );

    auto result2 = lu_lhs.setThreshold( 0.0 ).solve( rhs );

    std::cout << "result1: " << result1 << std::endl;

    std::cout << "result2: " << result2 << std::endl;

We use Eigen 3.3.1, the printed answer on my machine is the following:

result1:            

 0

  7.4746e+16

-3.25765e+11

      522139

result2: 

-4.74021e+23

 1.43104e+18

-3.25765e+11

      522139

result1 will fail the convergence, result2 can make the optimization continue.

3) I also have one question about the patch you provided, for a 4x4 matrix lhs, if lhs has rank 3, it will return a 3x1 vector, right? But it is supposed to return a 4x1 vector, will this cause any problem?

Thanks again.

On Wednesday, December 28, 2016 at 8:32:55 AM UTC-8, Sameer Agarwal wrote:

[Sameer Agarwal]

Thanks for the example Weiguang. This is a very interesting polynomial and  shows that my fix for the problem is broken. Because reducing the rank of the interpolating polynomial does not result in a better fit. In fact the QR factorization returns roughly the same result as the fullPivLU with its default threshold.

To investigate this a bit more carefully, would it be possible for you to share the arguments to MinimizeInterpolatingPolynomial that result in this matrix?

FunctionSample::ToDebugString() should do the trick.

finally you asked:

3) I also have one question about the patch you provided, for a 4x4 matrix lhs, if lhs has rank 3, it will return a 3x1 vector, right? But it is supposed to return a 4x1 vector, will this cause any problem?

The way the algorithm works is, we ask for an interpolating polynomial and then minimize it within the given bounds. The minimization algorithm does not care about the degree of the polynomial. 

Sameer

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/d5b3dac4-c75e-4a0e-9945-fa09bc9cea07%40googlegroups.com.

[weig...@light.co]

    Hi, Sameer, thanks for answering my question, here are the arguments to MinimizeInterpolatingPolynomial that result in the example matrix:

samples[0]

[x: 0.00000000e+00, value: 5.22139340e+05, gradient: -3.25764989e+11, value_is_valid: 1, gradient_is_valid: 1]

samples[1]

[x: 1.90750650e-06, value: 1.81769406e+06, gradient: -4.06082302e+10, value_is_valid: 1, gradient_is_valid: 1]

x_min: 0

x_max: 1.90751e-06

*optimal_x: 0

*optimal_value: 0

Please let me know if you need more information, thanks again.

[Sameer Agarwal]

Thanks Weiguang.

My initial experiments have left me more confused than before :/ so I need to look at this some more. Alex pointed out that most likely the same bug has brought down the performance of the LBFGS solver on the nist problems also. It was earlier able to solve 53 out of the 54 problems but is only able to solve 30 problems now.

So I have a local end to end reproduction of the failure too which I can use to debug this problem.

Sameer

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/436f42a0-e1a1-489c-873e-0bf92d3d836a%40googlegroups.com.

[Sameer Agarwal]

Weiguang,

I have committed the change suggested by you into the master branch while I figure out what is really going on.

Sameer

[weig...@light.co]

     Hi, Sameer, thanks a lot for taking care of this issue. Look forward to hearing more about this issue in future.

     Thanks again.