4

Say I want to efficiently evaluate $\sum_{kl}A_{ikjl}B_{kl}$ where $A$, $B$ are numerical tensors. This has been discussed before but with no focus on efficiency. A straightforward way as mentioned there would be

TensorContract[TensorProduct[A, B], {{2, 5}, {4, 6}}]

but this is extremely inefficient both in terms of memory and time.

I came up with

TensorProdContract = 
  Function[{A, B, dims}, 
   Transpose[#, RotateRight@Range@ArrayDepth[#]] &@
     Flatten[A, List@Transpose[dims][[1, All]]].Flatten[B, 
     List@Transpose[dims][[2, All]]]];

which flattens the contracted dimensions and rearranges the arrays so that the fast Dot function can be used. It would be called like

TensorProdContract[A, B, {{2, 1}, {4, 2}}]

But even so, for all dimensions sized at 50, this takes 0.15 seconds on my computer, whereas the tprod library for Matlab (basically compiled code) can do the same thing in 10 times shorter time (called as tprod(A, [1 -1 2 -2], B, [-1 -2])). For size 100, it's 3.5 seconds vs. 0.15 seconds. Given that this should work for arrays of general rank, the Compile approach cannot be really used as far as I know.

Any suggestions?

Community
  • 1
jhrmnn
  • 265
  • 1
  • 8
  • From profiling the code, most time (basically all the difference between Mathematica and Matlab) is spent on the Flatten/Transpose of the 4-D array. – jhrmnn Mar 17 '15 at 21:11
  • My method is equivalent to what @WReach came up with in there. It suffers in terms of speed from having to do all the list manipulation though. The answer by @jose below deals with the inefficiency of the TensorContract/TensorProduct approach. – jhrmnn Mar 18 '15 at 12:16
  • [I've mentioned now the very related question which does not discuss efficiency at all though. The answer below by @jose is the best for both my and the original question. What should happen now?] – jhrmnn Mar 18 '15 at 14:53
  • Have you tried WReach's Transpose method? It appears to be quite fast. Please include code to generate sample data and show that the available methods are still slow. – Mr.Wizard Mar 19 '15 at 10:56

1 Answers1

5

Take two arrays like you describe:

n = 100;
A = RandomReal[1, {n, n, n, n}];
B = RandomReal[1, {n, n}];

With your solution

AbsoluteTiming[
  Dimensions[result1 = TensorProdContract[A, B, {{2, 1}, {4, 2}}]]]

{2.22103, {100, 100}}

To avoid the inefficiency you refer to above, inactivate the tensor product and activate back at the end:

AbsoluteTiming[
  Dimensions[
    result2 = 
      Activate @ TensorContract[
        Inactive[TensorProduct][A, B], {{2, 5}, {4, 6}}]]]

{1.2759, {100, 100}}

result1 === result2

True
m_goldberg
  • 107,779
  • 16
  • 103
  • 257
jose
  • 6,328
  • 1
  • 14
  • 24
  • 1
    Wow, that's some use of Inactivate. Could you elaborate on how exactly does this work? I see that TensorProduct is left with literally nothing to do in this case. Does TensorContract basically crunch through the arguments of TensorProduct, figures out the result symbolically and replaces the arguments? I also see that if I replace TensorProduct with Inactive[someHead], this does not happen, does it mean that this behaviour is explicitly built into TensorContract to work with TensorProduct? It's a shame this is not documented if that's the case. – jhrmnn Mar 18 '15 at 11:55
  • 1
    I see this is still much slower than Dot for cases when Dot can be used. Do you see any way to combine my solution with Inactivate to harness the speed of Dot without doing all the list manipulation? Or would that be impossible because Dot is using some native-code library which directly crunches the blocks of reals as they are in the memory? – jhrmnn Mar 18 '15 at 11:58