I have some high dimensional high rank tensors, let's say $$F_{ijkl}$$ and I need to find $$F^{abcd}=g^{ai}g^{bj}g^{ck}g^{dl}F_{ijkl}.$$ Here $g^{ij}$ is the contravariant metric.

Simple summation in Mathematica takes way to much time


But I can sum over first and last index using matrix multiplication, so first I calculate $F^a{}_{ij}{}^d$ and the do summation over last two indices:

Fuddu=g.F.g; Do[F^{abcd}=Sum[Fuddu[[a,j,k,d]]g[[b,j]]g[[c,k]],{j,1,dim},{k,1,dim}],{a,1,dim},{b,1,dim},{c,1,dim},{d,1,dim}]

This way is much faster but still takes a lot of time, I need smth even faster, any ideas guys?

Edit: I cannot give you my info, but you can choose some big dim like 10, fill square matrix g of dimension dim and rank-4 tensor F of dimension dim with some random functions/numbers:

dim=10;g=RandomReal[{0, 1}, {dim, dim}];F=RandomReal[{0, 1}, {dim, dim,dim,dim}];
  • $\begingroup$ I ran your code as is, but I get a bunch of errors. dim is not specified and Fup, g is not initialized as an array, etc. Please add the additional information. $\endgroup$ – QuantumDot Jul 25 '15 at 18:09
  • 1
    $\begingroup$ I think the built-in function TensorContract will be helpful to you. $\endgroup$ – QuantumDot Jul 25 '15 at 18:22
  • $\begingroup$ It is only in Mathematica 10, and it is only symbolic, no? $\endgroup$ – Yuri Jul 25 '15 at 18:26
  • $\begingroup$ @g3n1uss Notice there is a type in the first computation: it should be ...Sum[F[[i,j,k,l]]... $\endgroup$ – Federico Aug 9 '15 at 23:34

My keyboard is broken. So here is fast answer (on Mathematica 9); more later...

Here is your input:

dim = 3; g = RandomReal[{0, 1}, {dim, dim}];
F = RandomReal[{0, 1}, {dim, dim, dim, dim}];

Now multiply four g's and the F. Use TensorProduct[g, g, g, g, F] (don't run this yet--it's slow) to generate the rank 12 tensor (unrepeated indices).

Now contract the 2nd and 9th indices, 4th and 10th indices, 6th and 11th indices, and also the 8th and 12th indices:

TensorContract[TensorProduct[g, g, g, g, F], {{2, 9}, {4, 10}, {6, 11}, {8, 12}}]

This is fastest I can get (0.3 sec on my machine).

  • $\begingroup$ Thank you very much, but I do not have TensorContract (Mathematica 8). $\endgroup$ – Yuri Jul 25 '15 at 21:27
  • $\begingroup$ aw :( buy upgrade $\endgroup$ – QuantumDot Jul 25 '15 at 21:28
  • 4
    $\begingroup$ When using TensorContract in a situation like this, it is important to inactivate the TensorProduct to avoid a very large intermediate expression. The timing of Activate[TensorContract[Inactive[TensorProduct][g, g, g, g, F], {{2, 9}, {4, 10}, {6, 11}, {8, 12}}]] is about 10^4 times faster than that without Activate/Inactive for dim=5. And the gain is even larger in higher dimensions. $\endgroup$ – jose Jul 26 '15 at 15:39


My solution is faster than the accepted one, even with the Activate/Inactive trick suggested in the comments.

Original answer

You can define your tensor contraction routine using the builtins Dot and Transpose. Here is an example:

DotAt[T_?TensorQ, U_?TensorQ, m_Integer?Positive, n_Integer?Positive] := 
    With[{dimT = ArrayDepth@T, dimU = ArrayDepth@U},
       Dot[Transpose[T, Insert[Range[dimT - 1], dimT, m]], 
           Transpose[U, Insert[Range[2, dimU], 1, n]]]]

DotAt[T, U, m, n] contracts the $m$-th index of $T$ with the $n$-th index of $U$. With this definition you have for example that Dot[T, U] is equivalent to DotAt[T, U, Length@Dimensions@T, 1].

From this, you can go on and define the equivalent of TensorContract, with an easier syntax for multiple contractions, but I'm leaving that to you :)

Your problem can now be solved by

myFup = DotAt[g, DotAt[g, DotAt[g, DotAt[g, F, 2, 4], 2, 4], 2, 4], 2, 4];

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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