Efficient matrix transpose matrix multiplication in Eigen

Question

I have access to a number of matrix libraries, but for this project I am using Eigen, due to its compile time definition and its inclusion of SVD.

Now, I am doing the following operation:

Eigen::Matrix<double,M,N> A;     // populated in the code

Eigen::Matrix<double,N,N> B = A.transpose() * A;

As I understand, this makes a copy of A and forms the transpose, which is multiplied by A again. This operation is being performed on relatively small matrices (M=20-30,N=3), but many millions of times per second, meaning it must be as fast as possible.

I read that using the following is faster:

B.noalias() = A.transpose() * A;

I could write my own subroutine that accepts A as an input and fills B, but I was wondering if there is an efficient, existing implementation that uses the least amount of cycles.

Answer 1

First of all, since Eigen relies on template expressions, A.transpose() does not evaluate into a temporary.

Secondly, in:

Matrix<double,N,N> B = A.transpose() * A;

Eigen knows that B cannot appear on the right hand side of the expression (because here the compiler calls the constructor of B), and therefore, no temporary is created at all. This is equivalent to:

Matrix<double,N,N> B;             // declare first
B.noalias() = A.transpose() * A;  // eval later

Finally, for such small matrices, I don't expect that the use of B.selfadjointView().rankUpdate(A) will help (as suggested in kennytm comment).

On the otherhand, with N=3, it might be worth trying the lazy implementation:

B = A.transpose().lazyProduct(A)

just to be sure. Eigen's has built-in heuristics to choose the best product implementation, but since the heuristic has to be simple and fast to evaluate, it might not be 100% right.