jeudi 26 juillet 2018

Eigen object value changes unexpectedly after expression

I'm struggling with a problem while using the Eigen3 library. This a very simplified version of my code:

#include <eigen3/Eigen/Dense>
#include <iostream>

using MyVec = Eigen::Matrix<double, 2, 1>;
using MyMat = Eigen::Matrix<double, 2, 2>;

class Foo {
public:
    Foo() = default;
    Foo(const MyVec &mean) : mean_(mean) { }
    Foo(MyVec&& mean) : mean_(std::move(mean)) { }

    Foo& operator+=(const MyVec& vec) 
    {
        mean_ += vec;
        return *this;
    }

    friend const Foo operator+(Foo lhs, const MyVec& rhs)
    {
        lhs += rhs;
        return lhs;
    }
private:
    MyVec mean_;
};

int main(){
    MyMat R(2*MyMat::Identity()),
          P(10*MyMat::Identity()),
          H(MyMat::Identity());
    MyVec residual, x_vec;

    x_vec << 2, 3;
    Foo x_pred(x_vec);

    residual << 0.1, 0.1;

    const auto S_inv = (H*P*H.transpose().eval() + R).inverse();

    const auto K = P*H*S_inv;
    std::cout << "K = " << K << std::endl;

    x_pred + K*residual;
    std::cout << "K = " << K << std::endl;

    return 0;
}

The output of the program is the following:

K = 0.833333        0
           0 0.833333
K = 0.00694444        0
    0.00694444 0.833333

Apparently, the value of K changes after the expression x_pred + K*residual, even if no one changed the variable, even if the variable is declared const. I did the following additional tests to try to find out the reason of this behavior:

  • using MyMat as type of the variable K;
  • replacing S_inv in the initialization of K with the expression (H*P*H.transpose().eval() + R).inverse();
  • changing the expression x_pred + K*residual with K*residual;
  • changing the expression x_pred + K*residual with x_vec + K*residual;

and none of them gave me the same unexpected behavior (i.e. the value of the variable K does not change).

Of course, I can adopt one of them as solution, but I'm curious to understand why this is happening. Does anyone have any idea?

I'm using g++-4.9.4 as compiler.

Thanks a lot in advance.

Aucun commentaire:

Enregistrer un commentaire