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& 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 eigen 3.2.0-8 and g++-4.9.4 as compiler.

Thanks a lot in advance.

EDIT

I tried to compile the same code with g++-4.7.3, clang++-3.6.0 and clang++-3.9.1; for all of them the code works as expected (K does not change). However, if I move to a newer version of g++ (I tried g++-4.9.4, g++-5.4.1 and g++-7.2.0), this weird behavior still remains. Could it be some compiler optimization then?

Solution

Short answer: Avoid using auto with Eigen expressions, unless you really know what you are doing (https://eigen.tuxfamily.org/dox/TopicPitfalls.html). Instead just write:

const MyMat S_inv = (H*P*H.transpose() + R).inverse(); // .eval() really did not help here
const MyMat K = P*H*S_inv;

Using .inverse() instead of solving a linear system is actually fine here, since there is an explicit formula used for 2x2 matrices.

A few more details: The following line

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

is equivalent to:

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

and auto will be some Eigen expression template (different, depending on the version of Eigen). Now H.transpose().eval() will create a temporary, which will get destructed at the end of the line, but is still referred to from the expression S_inv (which again is referred to from K. This means undefined behavior, i.e., you may or may not get a correct result, but your program could also just crash.

If you remove the .eval() your expression should actually work with Eigen 3.3 or higher (I think). It is still not a good idea to use auto, since expressions will be re-evaluated every time they are used.