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.
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.
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
auto
keyword, use explicit types.auto DF = []() -> MatrixXd {...}
instead of auto DF = []() {...}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With