Tensor multiplication

Question

I have eight Tensors to multiply as follows,

$$P=\sum_{all indices}M_{ijkl}M_{mjkl}M_{inkl}M_{mnkl}X_{kl}Y_{kl}X_{kl}Y_{kl}$$

Each M Matrix is say $2^7\times2^7 \times2^7 \times 2^7$ size. Is there any efficient way to perform this multiplication? Tensor contract even with Active and Inactive fails due to the exceptionally large Tensor(rank 10-12 ) that gets generated in between which requires huge amount of memory. The greater than 2 repeated indices is not a mistake, it is what it is.

Answer 1

I think, it is quite feasible. Let $N$ be the tensor dimension, $N=2^7$ in your case. I claim that the computational cost is $\mathcal{O}(N^5)$, which is around $3.2\times 10^{10}$, i.e., the tensor contraction can be computed within 1 minute on a laptop. Observe the following:

1. $M$ can be contracted with $X$ or $Y$ beforehand yielding $A_{ijkl}$, $B_{mjkl}$, $A_{inkl}$ and $B_{mnkl}$. This is just a $\mathcal{O}(N^4)$ operation, use Table for that.

2. Consider first the inner sum over $i,j,m,n$. This can be done sequentially via matrix multiplication (use . and Tr) as follows at the $\mathcal{O}(3 N^3)$ cost

$$T= A.B,\\ T= T.A,\\ T=T.B,\\ x_{kl}=\mathrm{Tr}(T). $$

3. Finally, one performs the sum $\sum_{kl}x_{kl}$ with the $\mathcal{O}(N^2)$ cost (Sum or ParallelSum).

Total computational cost is $\mathcal{O}(N^5)$. The space requirements are also very modest: one needs to store only 2 additional tensors $A$ and $B$ and a matrix $T$. Total additional storage $N^2(2 N^2+1)$, i.e., $\mathcal{O}(N^4)$.

The mathematica code could be as simple as A=Table[..]; B=Table[..]; Sum[a=A[[All,All,k,l]]; b=B[[All,All,k,l]]; Tr[a.b.a.b],{k,N},{l,N}]

Answer 2

Making the ideas in @yarchik's solution explicit:

using smaller versions of your matrices,

t = 4;
M = Array[mm, {t, t, t, t}];
X = Array[xx, {t, t}];
Y = Array[yy, {t, t}];

the exact sum you're looking for is

S = Sum[M[[i, j, k, l]] M[[m, j, k, l]] M[[i, n, k, l]] M[[m, n, k, l]]
        X[[k, l]] Y[[k, l]] X[[k, l]] Y[[k, l]],
        {i, t}, {j, t}, {k, t}, {l, t}, {m, t}, {n, t}];

define intermediates MX and MY:

MX = Transpose[Transpose[M, {3, 4, 1, 2}]*X, {3, 4, 1, 2}];
MY = Transpose[Transpose[M, {3, 4, 1, 2}]*Y, {3, 4, 1, 2}];

define an intermediate A: this step can probably also be done with a list-processing (linear algebra) operation instead of Table; but I can't figure it out right now,

A = Table[MX[[All, j, k, l]] . MY[[All, n, k, l]],
          {j, t}, {n, t}, {k, t}, {l, t}];

Now the sum is a scalar product:

S == Flatten[Transpose[A]] . Flatten[A] // Expand
(*    True    *)