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.


3 Answers 3


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}];
(*{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}];
ttt = t.tt;
{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.

  • $\begingroup$ Thank you! This was very helpful. $\endgroup$
    – Alehud
    Mar 18, 2021 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 "SymbolicTensors`ArrayContract".

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];
SymbolicTensors`ArrayContract[TP, {{3, 5, 7}}]

This creates a large output of the form:

enter image description here


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]

enter image description here

Have fun!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.