Efficient ways to generate a matrix from multiplying two elements of some matrix

Question

I am working on the single particle entanglement entropy and in the process I have something like

Table[
    Sum[
        states[[eindex, spaceindex]] * states[[eindex, spaceindex2]], 
        {eindex,1, Na/2}(*trace over lower-half filling*)
    ], 
    {spaceindex, 1, Nb}, 
    {spaceindex2, 1, Nb}(*some subsystem*)
]

where states is a Na by Na square matrix, from diagonalizing some other matrix. 1<Nb<Na is some sub-system. The matrix I am trying to construct is essentially $C_{i,j}=Tr[states[k,i]*states[k,j]]$ traced over the first index.

The above code gets the job done but I wonder if there are more efficient ways.

Also, I keep getting myself into this situation where I need to generate a matrix efficiently via similar constructions, for example, projection matrices. I would also like to know some general ways to directly translate the equations into efficient code.

Answer 1

Matrix multiplication is your friend because these operations are run by highly optimized libraries (BLAS).

Read out the matrices first, Transpose one of them and then Dot them together:

n = 1000;
Nb = 500;
Na = 400;
A = RandomReal[{-1, 1}, {n, n}];

ilist = Range[1, Nb];
jlist = Range[1, Nb];
klist = Range[1, Quotient[Na, 2]];

result0 = Table[Sum[A[[k, i]] A[[k, j]], {k, 1, Na/2}], {i, 1, Nb}, {j, 1, Nb}]; // AbsoluteTiming // First

result = Transpose[A[[klist, ilist]]] . A[[klist, jlist]]; // RepeatedTiming // First

Max[Abs[result0 - result]]

28.7556

0.000697215

2.84217*10^-14

As you can see, this is 40000 times faster with negligible error (which stems from then fact that floating point additions is not really associative).