Eigen and C++11 type inference fails for Cholesky of matrix product

Question

I am trying to take the cholesky decomposition of the product of a matrix with its transpose, using Eigen and C++11 "auto" type. The problem comes when I try to do

auto c = a * b
auto cTc = c.tranpose() * c;
auto chol = cTc.llt();

I am using XCode 6.1, Eigen 3.2.2. The type error I get is here.

This minimal example shows the problem on my machine. Change the type of c from auto to MatrixXd to see it work.

#include <iostream>
#include <Eigen/Eigen>
using namespace std;
using namespace Eigen;


int main(int argc, const char * argv[]) {
    MatrixXd a = MatrixXd::Random(100, 3);
    MatrixXd b = MatrixXd::Random(3, 100);
    auto c = a * b;
    auto cTc = c.transpose() * c;
    auto chol = cTc.llt();
    return 0;
}

Is there a way to make this work while still using auto?

As a side question, is there a performance reason to not assert the matrix is a MatrixXd at each stage? Using auto would allow Eigen to keep the answer as whatever weird template expression it fancies. I'm not sure if typing it as MatrixXd would cause problems or not.

Answer 1

Let me summarize what's is going on and why it's wrong. First of all, let's instantiate the auto keywords with the types they are taking:

typedef GeneralProduct<MatrixXd,MatrixXd> Prod;
Prod c = a * b;
GeneralProduct<Transpose<Prod>,Prod> cTc = c.transpose() * c;

Note that Eigen is an expression template library. Here, GeneralProduct<> is an abstract type representing the product. No computation are performed. Therefore, if you copy cTc to a MatrixXdas:

MatrixXd d = cTc;

which is equivalent to:

MatrixXd d = c.transpose() * c;

then the product a*b will be carried out twice! So in any case it is much preferable to evaluate a*b within an explicit temporary, and same for c^T*c:

MatrixXd c = a * b;
MatrixXd cTc = c.transpose() * c;

The last line:

auto chol = cTc.llt();

is also rather wrong. If cTc is an abstract product type, then it tries to instantiate a Cholesky factorization working on a an abstract product type which is not possible. Now, if cTc is a MatrixXd, then your code should work but this still not the preferred way as the method llt() is rather to implement one-liner expression like:

VectorXd b = ...;
VectorXd x = cTc.llt().solve(b);

If you want a named Cholesky object, then rather use its constructor:

LLT<MatrixXd> chol(cTc);

or even:

LLT chol(c.transpose() * c);

which is equivalent unless you have to use c.transpose() * c in other computations.

Finally, depending of the sizes of a and b, it might be preferable to compute cTc as:

MatrixXd cTc = b.transpose() * (a.transpose() * a) * b;

In the future (i.e., Eigen 3.3), Eigen will be able to see:

auto c = a * b;
MatrixXd cTc = c.transpose() * c;

as a product of four matrices m0.transpose() * m1.transpose() * m2 * m3 and put the parenthesis at the right place. However, it cannot know that m0==m3 and m1==m2, and therefore if the preferred way is to evaluate a*b in a temporary, then you will still have to do it yourself.