Hi! Hopefully this is the right place for this question, I have asked on stackoverflow and both answers were preceded by "I'm not an expert in Eigen but..." - they were certainly relevant but lacked the depth I'm hoping for here.

Here is the (non) working example:

Code: Select all

#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);
    ProductReturnType<MatrixXd,MatrixXd>::Type c = a * b;
    LLT<MatrixXd> cholcTc = (c.transpose() * c).llt();
    return 0;
}


This gives me the following error: https://gist.github.com/c0g/b8fa6863990a5dc334b0

The problem can be solved by calling eval on c.transpose() or forcing c to be a MatrixXd.

I realise that premature optimisation is a Bad Thing, however I have two questions that I can't work out how to answer by profiling (since it won't compile)

1) Can Eigen work out that c.transpose() * c = b.transpose() * (a.transpose() * a) * b? If so, that should reduce computational complexity considerably (goes from a 100x100 matrix product to one 100x3 matrix product and two 100x3 times 3x3 product). In numpy as a quick test with 1000x3 it goes from 126ms of cpu time to 10.7ms - pretty big difference!

2) If Eigen is this smart, how do I keep it as a expression right up until the end rather than using evals or casting to MatrixXd?

Tom

Detailed answer on SO (http://stackoverflow.com/questions/2711 ... ix-product).