Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Eigen gives wrong result when not storing intermediate result

Tags:

c++

eigen

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.

like image 468
David Avatar asked Jan 03 '20 23:01

David


1 Answers

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 = []() {...}
like image 122
David Avatar answered Oct 18 '22 11:10

David