18
$\begingroup$

One way to evaluate the following sums is combining Table and Sum:

$u_{abcd} = \sum_{e=1}^3 v_{aeb}w_{ced}$

$q_{ab} = \sum_{d,e=1}^3 v_{d e a}w_{deb}$

It will look like

v = Table[Times[i, j, k], {i, 3}, {j, 3}, {k, 3}]
w = Table[Times[i+2, j-1, k+1], {i, 3}, {j, 3}, {k, 3}]

u = Table[Sum[v[[a, e, b]] w[[c, e, d]], {e, 3}], {a, 3}, {b, 3}, {c, 3}, {d,  3}]
q = Table[Sum[v[[d, e, a]] w[[d, e, b]], {d, 3}, {e, 3}], {a, 3}, {b, 3}]

Is there a more elegant way to evaluate sums like these?

$\endgroup$
5
  • 7
    $\begingroup$ Just replacing the sum with Dot as in u2 = Outer[Dot[v[[#1, All, #2]], w[[#3, All, #4]]] &, Range[3], Range[3], Range[3], Range[3]]. $\endgroup$ Commented Jan 14, 2013 at 15:27
  • $\begingroup$ @b.gatessucks I'm curious why you didn't post that as an answer. I took the liberty of posting it for you, but I encourage you to post such things yourself in the future. If you wish to post it now, I'll be glad to delete mine. $\endgroup$
    – Mr.Wizard
    Commented Jan 14, 2013 at 15:33
  • 1
    $\begingroup$ @Mr.Wizard It's not gonna be more helpful in an answer and it's quicker to add a comment. Thanks for taking the time. $\endgroup$ Commented Jan 14, 2013 at 15:40
  • 2
    $\begingroup$ @b.gatessucks actually, it is: comments cannot be Accepted and voting is limited. It takes hardly a moment longer to make an answer than a comment, and you should get the credit for your own answers/ideas. Nevertheless, it's your prerogative. $\endgroup$
    – Mr.Wizard
    Commented Jan 14, 2013 at 15:42
  • 3
    $\begingroup$ From the tutorial on tensors: "You can think of Inner as performing a "contraction" of the last index of one tensor with the first index of another. If you want to perform contractions across other pairs of indices, you can do so by first transposing the appropriate indices into the first or last position, then applying Inner, and then transposing the result back." For multiple contractions, as in the second example, transpose all involved indexes to the end, Flatten them into a single index, and then perform a contraction. $\endgroup$
    – whuber
    Commented Jan 14, 2013 at 15:58

5 Answers 5

19
$\begingroup$

If we transpose the indices of $v$ and $w$ so that $v'_{abe} = v_{aeb}$ and $w'_{ecd} = w_{ced}$, then we can compute $u = v' \cdot w'$:

u === Transpose[v, {1, 3, 2}] . Transpose[w, {2, 1, 3}]
(* True *)

We can use a similar trick to compute $q$ if we reorder $v_{dea} \to v''_{ade}$, except that this time the $d$ and $e$ indices in $v_{ade}''$ and $w_{deb}$ must treated as if they comprised a single index to be contracted. Flatten can do this for us:

q === Flatten[v, {{3}, {1, 2}}] . Flatten[w, {{1, 2}, {3}}]
(* True *)

For details about the Flatten syntax used here, see Flatten command: matrix as second argument.

$\endgroup$
2
  • 2
    $\begingroup$ Much better. I was trying to figure this out myself earlier but failed. Big +1. $\endgroup$
    – Mr.Wizard
    Commented Jan 15, 2013 at 1:49
  • 1
    $\begingroup$ For anyone coming here in more modern times (i.e. after the advent of TensorContract, TensorProduct, and crucially Inactive) there is a later answer by Bren that reformulates this answer by Jose Martin-Garcia and which is much more convenient and much faster for large tensors. $\endgroup$
    – b3m2a1
    Commented Jun 30, 2019 at 5:32
19
$\begingroup$

In Mathematica version 9, you can do these kinds of things much more naturally. Here are the two quantities that you wanted:

u = TensorContract[TensorProduct[v, w], {2, 5}];

q = TensorContract[TensorProduct[v, w], {{1, 4}, {2, 5}}];

The contraction is performed on the tensor product in which the first three indices belong to the factor v and the last three indices label w. Therefore, the indices corresponding to $e$ in your sum for u are in slots {2, 5}, and the two summations for q run over slots {1, 4} (for the variable $d$) and {2,5} (as in v).

$\endgroup$
8
  • $\begingroup$ No love for this answer? If it works it's nice. (I can't test it.) $\endgroup$
    – Mr.Wizard
    Commented Jan 15, 2013 at 6:27
  • 1
    $\begingroup$ @Mr.Wizard Trust me, it works, it's the best. I consider this to be one of the most important news in ver.9, much more interesting than auto-completion tosh. $\endgroup$
    – Artes
    Commented Jan 15, 2013 at 12:56
  • $\begingroup$ @Artes Thanks - and I still haven't gotten used to autocompletion; when it appears and I want to just navigate away by down-arrow, it annoyingly snags the cursor... I guess we have to take the good with the bad. $\endgroup$
    – Jens
    Commented Jan 15, 2013 at 15:01
  • $\begingroup$ @Jens Since I've been getting used to autocompletion for 2 months, it seems quite natural and sometimes I lack it when using ver.8. Nevertheless I find it's a toy unlike the new tensor capabilites which are definitely still underestimated. $\endgroup$
    – Artes
    Commented Jan 17, 2013 at 16:21
  • 3
    $\begingroup$ @Jens Doesn't this waste memory in a huge way because the calculation stores the tensor product as an intermediate value? Or does TensorProduct have a special UpValue or something... $\endgroup$
    – Ian Hincks
    Commented Oct 19, 2013 at 0:36
12
$\begingroup$

Since the question "Efficient tensor product followed by contraction" asking for an efficient solution to this problem has been marked as a duplicate, I find it appropriate to add here an encapsulated version of the answer to that question by @m_goldberg. Note that this works efficiently for the contraction of any number of index pairs. The notation follows Mathematica's TensorContract.

TensorMultiply[ A, B, {{$a_1$,$b_1$},{$a_2$,$b_2$},...} ]

$\hspace{7.5mm}$yields the tensor product of tensors A and B in the pairs {$a_i$,$b_i$} of slots.

TensorMultiply[A_, B_, pairs_] := 
  Activate@TensorContract[
    Inactive[TensorProduct][A, B], (# + {0, TensorRank[A]}) & /@ pairs];
$\endgroup$
4
  • $\begingroup$ I must have overlooked this earlier - it works very well (+1). $\endgroup$
    – Jens
    Commented May 15, 2016 at 0:07
  • $\begingroup$ This is both quite clever and very efficient! (It makes sense that this works as well as it does as it was proposed by Jose Martin-Garcia). For high dimensional tensors this vastly outperforms the Transpose approach. It is also vastly more convenient. $\endgroup$
    – b3m2a1
    Commented Jun 30, 2019 at 5:30
  • $\begingroup$ Can this be generalized to more than two tensors? I thought about if for a minute, but the generalization isn't completely trivial because you have to worry about index pairs that refer to different tensor pairs. $\endgroup$ Commented May 12 at 17:54
  • $\begingroup$ OK, I made the generalization and submitted it as a suggested edit addition. $\endgroup$ Commented May 30 at 20:07
5
$\begingroup$

Incorporating b.gatessucks's recommendation:

u === Outer[v[[#1, All, #2]].w[[#3, All, #4]] &, ##] & @@ Range@{3, 3, 3, 3}

True

$\endgroup$
1
  • $\begingroup$ How would you calculate q with Outer[] and Dot[]? $\endgroup$
    – sjdh
    Commented Jan 15, 2013 at 10:43
4
$\begingroup$

Edit

This is the fastest solution:

Do[TensorContract[TensorProduct[v, w], {2, 5}], {10000}] // Timing // AbsoluteTiming
Do[TensorMultiply[v, w, {{2, 2}}], {10000}] // Timing // AbsoluteTiming
Do[Activate@TensorContract[Inactive[TensorProduct][v, w], {2, 5}], {10000}] // Timing // AbsoluteTiming
Do[DotAt[v, w, 2, 2], {10000}] // Timing // AbsoluteTiming
{5.36688, {5.36012, Null}}
{0.753744, {0.752581, Null}}
{0.57401, {0.573111, Null}}
{0.143519, {0.143405, Null}}

Original answer

As used in a recent answer recent answer, one can easily define a function which contracts two tensors:

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

DotAt contracts the index m of T with the index n of U.

This way, u = TensorContract[TensorProduct[v, w], {2, 5}]; would become u = DotAt[v, w, 2, 2]. This just shows that Dot and Transpose can be combined to perform all the operations and that it is actually quite easy to implement a form of TensorContract with the syntax that one prefers.

$\endgroup$

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.