Consider a tensor $T\in\mathbb{R}^{N\times N\times N\times M}$ and two vectors $x,y\in\mathbb{R}^N$. I want to compute the $N\times M$ vector defined by $X_{ij}=\operatorname{tr}(x^\top T_{:,:,i,j}y)=\operatorname{tr}_{12}(yx^\top T)$ efficiently.
I tried this in two different ways:
TCtable[x_, T_, y_] := ParallelTable[Module[{Tslice},
Tslice = T[[;; , ;; , i, j]];
Tr[x\[Transpose].Tslice.y] ], {i, 1, Length[x]}, {j,1,Last[Dimensions[T]]}];
TCtrace[x_, T_, y_] := TensorContract[y.x\[Transpose].T, {1, 2}];
My tensors have the very nice property that $T_{:,:,i,j}$ is very sparse for all $i,j$ (so I am representing my tensor in Mathematica as a sparse array).
With $N=500,M=3$ the parallel table method takes about 1 second, while the explicit tensor multiplication and partial trace takes about 20 seconds. Are there other clever ways to speed this up? I want to compute this for many different tensors and vectors of the same size so if there's a way to compile the code or amortize the complexity that would also be great!
ParallelTable
, but because the second one needs times-$N^2$ operations. $\endgroup$