Tensor multiplication

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.

• An explicit $M$ is appreciated, I guess. – Αλέξανδρος Ζεγγ Apr 29 at 7:38
• You can take it to be a completely randomreal matrix – Roopayan Ghosh Apr 29 at 17:28

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).$$

1. 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}]

• This seems easily achievable thanks. I shall do it and check, shall accept it as an answer once I do it. – Roopayan Ghosh Apr 29 at 21:35
• @RoopayanGhosh Have you already checked it? – yarchik May 4 at 18:11
• Oh thanks for reminding, I did and it worked, accepting the answer. Sorry I forgot – Roopayan Ghosh May 5 at 19:32

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    *)

• A is rank-4 and the whole thing scales as yours, as far as I can see. – Roman Apr 29 at 16:47
• This seems to have same complexity to me as well. Thanks for the explicit answer. – Roopayan Ghosh Apr 29 at 21:37
• @Roman Sorry, I misunderstood your method. It became completely clear after you added All. – yarchik Apr 30 at 5:21