Eigen gives wrong result when not storing intermediate result
Asked Answered
F

1

6

Writing a function implementing a Jacobi Matrix for Newton's method I have noticed a really nasty error.

Calling the function

auto DF = [T](VectorXd y){
    return PhiAndW(y(0),y(1),T).second - MatrixXd::Identity(2,2);
 };

only returns the value of PhiAndW(y(0),y(1),T).second and omits the subtraction of MatrixXd::Identity(2,2). But if I change the code to

auto DF = [T](VectorXd y){
    MatrixXd mat = PhiAndW(y(0),y(1),T).second - MatrixXd::Identity(2,2);
    return mat;
 };

everything works seamlessly.

I have tried to reproduce it and this is not the exact same behavior but it also doesn't behave as wanted:

Consider the following code:

MatrixXd FF(MatrixXd y){
  return y;
}

int other(){

  auto DF = [](MatrixXd y){
    MatrixXd test = FF(y)  - MatrixXd::Identity(2,2);
    return test;
  };

  std::cout << DF(MatrixXd::Ones(2,2)) <<std::endl;
  std::cout << std::endl;
  std::cout << (MatrixXd::Ones(2,2) - MatrixXd::Identity(2,2))<< std::endl;

}

It will print

>  1 0
>  0 1 
>
>  1 0 
>  0 1 

to the console.

However if I change to function DF to

  auto DF = [](MatrixXd y){
    return FF(y)  - MatrixXd::Identity(2,2);
  };

The console will print

>  2.22045e-15 1.63042e-322
>  2.37152e-322    -0.999998

for the second matrix. (Which is just some uninitialized junk from memory).

Could someone explain what is happening with my code and my example problem. I have truly no idea what is going on here. I am especially interested why saving the result of the calculation in a temporary variable fixes the problem.

Fanciful answered 3/1, 2020 at 23:27 Comment(5)
Please read the section about using auto on this page: eigen.tuxfamily.org/dox-devel/TopicPitfalls.htmlOutrun
(in this case auto is implicitly the return type of the lambda)Commemorative
Did you try auto DF = []() -> MatrixXd {...} ? This should force the evaluation of return value.Sile
That's another nice example that auto might not be that cool like everyone else (except me) seems to believe. ;-) I assume an explicit cast would fix the first sample as well: auto DF = [](MatrixXd y){ return MatrixXd(FF(y) - MatrixXd::Identity(2,2)); };Butylene
@Scheff You're not the only one.Trenna
F
4

Since the comments have pretty much solved my issue (thank you very much) I thought I'd go ahead and answer my question so that other people see that this problem has been solved.

Why does the problem occur?

The problem is that, for example the result type of an Eigen multiplication of two matrices is not an Eigen matrix, but some internal object which represents the multiplication and references the two matrices we are trying to multiply.

Hence, if we use the auto keyword the compiler will very likely not give the variable we are setting the type MatrixXd but the type of some internal object.

For more information refer to the Eigen documentation which explicitly states:

In short: do not use the auto keywords with Eigen's expressions, unless you are 100% sure about what you are doing. In particular, do not use the auto keyword as a replacement for a Matrix<> type

How do I prevent it from happening?

  • do not use the auto keyword, use explicit types.
  • for a lambda function always specify the return type: Write auto DF = []() -> MatrixXd {...} instead of auto DF = []() {...}
Fanciful answered 4/1, 2020 at 11:23 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.