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.

Breadroot

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.

Danielledaniels

Thanks - it's really great to hear from a developer of the library! My reason for this was to see if Eigen could correctly optimise m0.transpose() * m1.transpose() * m2 * m3 when they have useful shapes - hence I wanted to keep everything in expression-space until the last minute. Is it due to how templates work that I can't take a cholesky decomposition of a GeneralProduct, is it just that no-one has cared enough to add it to Eigen yet or is there a reason why doing it is stupid? – Breadroot

Answer 2

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.

Repine

Is there a performance hit for casting to MatrixXd? Obviously in this constrained example it will be minimal, but in real life? What would you say the solution is here - using ProductReturnType for this would seem to go a bit mad when doing the c.tranpose() * c step. – Breadroot

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.)

Me

Answer 4

DO NOT use auto with Eigen expressions. I bumped into even more "dramatic" issues with this before, see

eigen auto type deduction in general product

and was advised by one of the Eigen creators (Gael) not to use auto with Eigen expressions.

The cast from an expression to a specific type like MatrixXd should be extremely fast, unless you want lazy evaluation (since when doing the cast the result is evaluated).

Glottis

I think this is slightly different to yours. You return an object from Id that is in turn referenced by autoresult - then the object from Id goes away and autoresult references something that doesn't exit. c.transpose() references c, which still exists so my auto cTc shouldn't reference anything that doesn't exist - similar to how ggael offers as one solution to use auto id = Id(Foo, 4). I wanted to test if Eigen can correctly notice that c.transpose() * c = b.transpose() * a.transpose() * a * b which simplifies the math, so ideally it would be kept as an expression. – Breadroot

yes, you are correct, however, if I haven't got the explanation, I would have had absolutely no idea why the code behaved abnormally. And since Eigen has a whole lot of types (basically every expression is a different type), plus indirect referencing, auto makes things complicated as you are not really aware of what goes on under the hood, and how are these expressions evaluated. – Glottis