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.

Douce

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;

Walke

Thanks much! One last question: Isn't the temporary object stored via a const reference binding to the temporary, so the lifetime of the temporary is the same as the lifetime of the expression? – Douce

Yes the temporary is stored via a const reference, but this does not extent its lifetime because there is one indirection. This lifetime extension mechanism is not transitive. – Walke

Answer 2

It probably has a lazy evaluation type that is only safe to evaluate once. You could capture it with:

auto autoresultmatrix = autoresult.eval()

Bronchiectasis

Thanks, I know that if I force evaluation it works (when doing a cast the evaluation is forced), however this makes the code much uglier in general, especially when dealing with complicated expressions. I would have liked to write something like myfunc(A*A+B-C), where myfunc is some function that takes an Eigen expression and spits out the result of the computation as a dynamic matrix, and I though that it will be perfectly fine to use it later in expressions like myfunc(A*A+B-C) * D, without the need for eval. Especially when calling methods like norm, which SHOULD force evaluation. – Douce

@Douce The transpose + print is evaluating and consuming the result. If you did the norm first it would work. I suggested the eval because I thought your concern was being unable to use auto. – Bronchiectasis

And even if I first do std::cout << result.norm(), still get the same behaviour, i.e. junk is being displayed. – Douce

Basically my question is how to write a function that takes an arbitrary Eigen expression as a parameter, then computes the result ALWAYS (so when I use cout<<f(expression) I do not have lazy evaluation, but want f(expression) to be evaluated). – Douce

And as a last comment, if instead of the autoresult I use directly (Id(Foo, 4) * v).transpose(), then it works. I am really a bit puzzled, as I see no difference between (Id(Foo, 4) * v).transpose() and autoresult.transpose(), except that in the first case one has a temporary. – Douce