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

The problem is that the first multiplication returns a Eigen::GeneralProduct instead of a MatrixXd and auto is picking up the return type. You can implicitly create a MatrixXd from a Eigen::GeneralProduct so when you explicitly declare the type it works correctly. See http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralProduct.html

EDIT: I'm not an expert on the Eigen product or performance characteristics of doing the casting. I just surmised the answer from the error message and confirmed from the online documentation. Profiling is always your best bet for checking the performance of different parts of your code.

Answer 2

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.

Answer 3

I'm not an expert at Eigen, but libraries like this often return proxy objects from operations and then use implicit conversion or constructors to force the actual work. (Expression Templates are an extreme example of this.) This avoids unnecessary copying of the full matrix of data in many situations. Unfortunately, auto is quite happy to just create an object of the proxy type, which normally users of the library would never explicitly declare. Since you need to ultimately have the numbers calculated, there is not a performance hit per se from casting to a MatrixXd. (Scott Meyers, in "Effective Modern C++", gives the related example of using auto with vector<bool>, where operator[](size_t i) returns a proxy.)