# Contraction of three tensors with the same index

I have a 3-dimensional tensor $$T_{ijk}$$. I need to calculate a tensor $$M_{ijnmpq}=\sum_k T_{ijk}T_{nkm}T_{kpq}$$ Is there a way to do such contractions in Mathematica, avoiding loops? Or maybe there is a convenient way of rewriting this sum in terms of usual tensor contractions where indices of summation appear only twice? Thanks.

Very beautiful question. Indeed, there exist an efficient and concise method.

1. Let us first do the brute force calculation

n = 4;
t = RandomReal[{0, 1}, {n, n, n}];
TTT = Table[
Sum[ t[[i, j, k]] t[[p, k, q]] t[[k, r, s]], {k, n}],
{i, n}, {j, n}, {p, n}, {q, n}, {r, n}, {s, n}];
Dimensions[TTT]
(*{4, 4, 4, 4, 4, 4}*)


2. Now we build an auxiliary rank-5 tensor and take a dot product of the original rank-3 tensor and the auxiliary rank-5 tensor leading to the rank-6 result. The numerical complexity is $$\text{Multiplications: }\mathcal{O}(n^7),$$ which is optimal in this case. But the biggest benefit is that one avoids the storage of a rank-9 tensor as in the TensorProduct[t, t, t] approach. $$\text{Storage: }\mathcal{O}(n^6)\text{ (present) vs. } \mathcal{O}(n^9) \text{ (naive).}$$

tt = Table[TensorProduct[t[[;; , k, ;;]], t[[k, ;; , ;;]]], {k, n}];
Dimensions[tt]
ttt = t.tt;
Dimensions[ttt]
{4, 4, 4, 4, 4}
{4, 4, 4, 4, 4, 4}


3. Now we verify the calculation

Norm[Flatten[ttt - TTT]]
ttt == TTT
(* 4.15065*10^-15 *)
(* True *)


The timing on my rather old laptop and $$n=10$$ is for the fast method $$0.0053\,s$$ vs. $$24.7952\,s$$ for the brute force, the speedup around 5000.

• Thank you! This was very helpful. – Alehud Mar 18 at 1:21

One possibility is to create the tensor product and then do a contraction over the 3., 5. and 7. index.

MMA has 2 high level functions for this: "TensorProduct" and "SymbolicTensorsArrayContract".

Assuming that your tensor: T is given as an array of rank 3 I create a symbolic tensor for this example, but you can take any numeric tensor:

T = Array[#1 #2 #3 &, {3, 3, 3}];


We now create the tensor product and the contraction:

TP = TensorProduct[T, T, T];
SymbolicTensorsArrayContract[TP, {{3, 5, 7}}]


This creates a large output of the form:

Try this:

Subscript[M, i_, j_, n_, m_, p_, q_] :=
Sum[Subscript[T, i, j, k]*Subscript[T, n, k, m]*Subscript[T, k, p,
q], {k, 1, 3}]


Then, for example,

Subscript[M, 1, 1, 1, 1, 1, 1]


Have fun!