Inference like Eigen and C ++ 11 fails for Cholesky matrix product

Question

I am trying to take a choleck decomposition of a matrix product with its transposition using Eigen and C ++ 11 "auto". The problem occurs 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 a problem on my machine. Change type c from auto to MatrixXd to see how it works.

 #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 when using auto?

As a side question, is there a performance reason for not claiming that the matrix is MatrixXd at every stage? Using auto would allow Eigen to save the response like any strange template expression that it represents. I'm not sure if typing it, since MatrixXd will cause problems or not.

Answer 1

Let me summarize what happens and why it is wrong. First of all, create auto instances with the types they accept:

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

Please note that Eigen is an expression template library. Here, GeneralProduct<> is an abstract type representing the work. No calculations are performed. Therefore, if you copy cTc to MatrixXd as:

 MatrixXd d = cTc; 

which is equivalent to:

 MatrixXd d = c.transpose() * c; 

then the product a*b will be executed twice! Therefore, in any case, it is much preferable to evaluate a*b within the explicit time value, as well as for c^T*c :

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

Last line:

 auto chol = cTc.llt(); 

also wrong. If cTc is an abstract product type, it tries to create a Cholesky factorization instance that works on an abstract product type that is not possible. Now, if cTc is MatrixXd , then your code should work, but this is still not the preferred way, since the llt() method is more likely to implement a single-line expression, for example:

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

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

 LLT<MatrixXd> chol(cTc); 

or even:

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

which is equivalent if you should not use c.transpose() * c in other calculations.

Finally, depending on the sizes of a and b , it would be preferable to calculate 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 the product of the four matrices m0.transpose() * m1.transpose() * m2 * m3 and put the brackets in the right place. However, he cannot know that m0==m3 and m1==m2 , and therefore, if the preferred method is to estimate a*b in a temporary order, you still have to do it yourself.

Answer 2

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

EDIT: I am not an expert on the Eigen product or casting specifications. I just assumed the answer from the error message and confirmed it in the online documentation. Profiling is always the best choice to test the performance of different parts of your code.

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 them to work. (Exam templates are an extreme example of this.) This avoids unnecessary copying of the complete data matrix in many situations. Unfortunately, auto is quite happy to just create an object like a proxy, which usually library users will never explicitly declare. Since you need to ultimately calculate the numbers, there is no effect associated with the cast on MatrixXd . (Scott Meyers, in Effective Modern C ++, gives a related example of using auto with vector<bool> , where operator[](size_t i) returns the proxy.)

Answer 4

DO NOT use auto with Eigen expressions. I ran into even more “dramatic” problems with this before, see

independent automatic type inference in a common product

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

MatrixXd from an expression to a specific type of type MatrixXd should be extremely fast if you do not want a lazy evaluation (since the result is evaluated during the task).