c++ - eigen auto type deduction in general product

Question

I have the following piece of code (I apologize for the slightly larger code snippet, this is the minimal example I was able to reduce my problem to):

#include <Eigen/Dense>
#include <complex>
#include <iostream>
#include <typeinfo>

// Dynamic Matrix over Scalar field
template <typename Scalar> 
using DynMat = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;

// Dynamic column vector over Scalar field
template <typename Scalar>
using DynVect = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;

// Returns the D x D Identity matrix over the field Derived::Scalar
// deduced from the expression Eigen::MatrixBase<Derived>& A
template<typename Derived>
DynMat<typename Derived::Scalar> Id(const Eigen::MatrixBase<Derived>& A, std::size_t D)
{   
    DynMat<typename Derived::Scalar> result =
            DynMat<typename Derived::Scalar>::Identity(D, D);

    return result;
}

int main()
{
    //using ScalarField = std::complex<double>; // same issue even if I use complex numbers
    using ScalarField = double; // we use doubles in this example

    // A double dynamic matrix (i.e. MatrixXd)
    DynMat<ScalarField> Foo; // used to deduce the type in Id<>()

    // A double dynamic column vector (i.e. VectorXd)
    DynVect<ScalarField> v(4);
    v << 1., 0. , 0. ,0.; // plug in some values into it

    // Make sure that Id(Foo, 4) correctly deduces the template parameters
    std::cout << "Id(Foo, 4) is indeed the 4 x 4 identiy matrix over the ScalarField of "
              << "typeid().name(): " << typeid(ScalarField).name() << std::endl;
    std::cout << Id(Foo, 4) << std::endl; // Indeed the 4 x 4 complex Identity matrix

    // Use auto type deduction for GenMatProduct, junk is displayed. Why?!
    std::cout << std::endl << "Use auto type deduction for GenMatProduct,
                 sometimes junk is displayed. Why?!" << std::endl;
    auto autoresult = Id(Foo, 4) * v; // evaluated result must be identically equal to v
    for(int i=0; i<10; i++)
    {
            std::cout << autoresult.transpose(); // thought 1 0 0 0 is the result, but NO, junk
            std::cout << " has norm: " << autoresult.norm() << std::endl; // junk
    }

    // Use implicit cast to Dynamic Matrix, works fine
    std::cout << std::endl << "Use implicit cast to Dynamic Matrix, works fine" << std::endl;
    DynMat<ScalarField> castresult = Id(Foo, 4) * v; // evaluated result must be identically equal to v
    for(int i=0; i<10; i++)
    {
            std::cout << castresult.transpose(); // 1 0 0 0, works ok
            std::cout << " has norm: " << castresult.norm() << std::endl; // ok
    }
}

The main idea is that the template function Id<>() takes an Eigen expression A as a parameter, together with a size D, and produces the identity matrix over the scalar field of the expression A. This function by itself works fine. However, when I use it in an Eigen product with auto deduced type, such as in the line auto autoresult = Id(Foo, 4) * v, I would expect to multiply the vector v by the identity matrix, so the net result should be an expression which, when evaluated, should be identically equal to v. But this is not the case, see the first for loop, whenever I display the result and computes its norm, I get most of the time junk. If, on the other hand, I implicitly cast the Eigen product Id(Foo, 4) * v to a dynamic matrix, everything works fine, the result is properly evaluated.

I use Eigen 3.2.2 on OS X Yosemite, and get the same weird behaviour both with g++4.9.1 and Apple LLVM version 6.0 (clang-600.0.54) (based on LLVM 3.5svn)

QUESTION:

• I do not understand what is happening in the first for loop, why isn't the product evaluated when I use std::cout, or even when I use the norm method? Am I missing something? There is no aliasing involved here, and I am really puzzled by what is going on. I know that Eigen uses lazy evaluation, and evaluates the expression when needed, but this doesn't seem to be the case here. This problem is extremely important for me, as I have lots of functions of the same flavour as Id<>(), which when used in auto deduced expressions may fail.

The problem occurs quite often, but not always. However, if you run the program 3-4 times, you will definitely see it.

The command I use to compile and run it is:

clang++ (g++) -std=c++11 -isystem ./eigen_3.2.2/ testeigen.cpp -otesteigen; ./testeigen

A typical output I got in a real run is:

Id(Foo, 4) is indeed the 4 x 4 identiy matrix over the ScalarField of typeid().name(): d
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

Use GenMatProduct, sometimes junk is displayed. Why?!
1 0 0 0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf

Use implicit cast to Dynamic Matrix, works fine
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1

Even if I use eval() in

  std::cout << autoresult.eval().transpose(); // thought 1 0 0 0 is the result, but NO, junk
  std::cout << " has norm: " << autoresult.eval().norm() << std::endl; // junk

I am getting the same weird behaviour.

Answer 1

The problem is that Id() returns a temporary which is stored by reference in the object representing the expression Id(Foo, 4) * v. Thus after the auto statement, autoresult stores a reference to a dead object. If you do not want an abstract expression but the actual result, do not use auto or call eval to enforce evaluation:

auto autoresult = (Id(Foo, 4) * v).eval();

A third option is to make the object returned by Id() available for further computations:

auto id4 = Id(Foo,4);
auto autoresult = id4 * v;

but in this case, anytime you use autoresult then the product will be re-evaluated and the following will output different results:

cout << autoresult;
v.setRandom();
cout << autoresult;