Hello,

I'm trying to compute Tr(AB) efficiently where A, B are Hermitian (self-adjoint) matrices and Tr is the matrix trace.

The current code I've got is:

Code: Select all

real((a.array() * b.array()).sum());


Noting that Tr(AB) = AB_{ii} = A_{ij}B_{ji} = A_{ij}B_{ij}, via Einstein's index convention. I'd like to use this trick since it doesn't use matrix multiplication. Furthermore, I'd like to use the triangularView mechanism to exploit the fact that A and B are Hermitian. I want to sum the lower off-diagonal terms, multiply that by 2, and add the sum of the diagonal terms. How does one go about this?

In theory you could write something like:

Code: Select all

F.cwiseProduct(G).triangularView<StrictlyUpper>().sum()*2 + F.cwiseProduct(G).diagonal().sum()


but reductions on triangular-views are not implemented yet.

Note that because of the overhead to vectorize a reduction on a triangular part, speedup would be expected for large matrices only, like larger than 100x100.

To get an idea of what could be the speed up, you can try a manual for loop:

Code: Select all

tr = 0;
for(Index j=1; j<A.cols(); ++j)
  tr += A.cwiseProduct(B).cols(j).head(j);
tr = tr*2 + A.cwiseProduct(B).diagonal().sum()


Let us know if you get a significant speed up. If so, this might motivate us to implement triangular reductions!