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
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.
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:
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!