4
$\begingroup$

I am interested in the tensor product $\hat{B} = A \star B$ (which at least I know as Rayleigh product), defined with components

\begin{equation} \hat{B}_{i_1 i_2 ... i_ n} = \sum_{j_1 = 1}^d \sum_{j_1 = 1}^d ... \sum_{j_n = 1}^d A_{i_1 j_1} A_{i_2 j_2} ... A_{i_n j_n} B_{j_1 j_2 ... j_n} \ . \end{equation}

How would you implement this efficiently for tensors of arbitrary order? The following rude implementation using TensorProduct and TensorContract is just slow

rp1[A_, B_] := Block[
   {n, temp, ind},
   n = TensorRank[B];
   temp = ConstantArray[A, n];
   AppendTo[temp, B];
   ind = Table[{2*k, 2*n + k}, {k, n}];
   TensorContract[TensorProduct @@ temp, ind]
   ];

For example

n = 4; (*tensor order*)
d = 4; (*dimension*)
A = RandomInteger[{-10, 10}, {d, d}];
B = RandomInteger[{-10, 10}, ConstantArray[d,n]];
AbsoluteTiming[res1 = rp1[A, B];][[1]]

0.146019

For fourth-order tensors you can do this as

rp2[A_, B_] := Block[
   {dim, i1, i2, i3, i4},
   dim = Length[A];
   Table[
    B.A[[i4]].A[[i3]].A[[i2]].A[[i1]]
    , {i1, dim}, {i2, dim}, {i3, dim}, {i4, dim}]
   ];

which is a lot faster

AbsoluteTiming[res2 = rp2[A, B];][[1]]
res1 == res2

0.0023153

True

I am sure, I just dont have the right perspective on this, but how would you implement this for tensors of arbitrary order?

EDIT: This is just based on yarchik's answer (please see his answer), just to give the code for general tensor order/rank

rp[A_, B_] := Block[
   {n, it, t1},
   n = TensorRank[B];
   it = RotateLeft@Range[n];
   t1 = B;
   Do[t1 = TensorTranspose[A.t1, it], {i, n}];
   t1
   ];

Testing

n = 4; (*tensor order*)
d = 10; (*dimension*)
A = RandomInteger[10, {d, d}];
B = RandomInteger[10, ConstantArray[d, n]];
AbsoluteTiming[res1 = rp[A, B];][[1]]
AbsoluteTiming[res2 = rp2[A, B];][[1]]
res1 == res2

0.00277906

0.770904

True

$\endgroup$
6
$\begingroup$

One can still speed up your code a lot

rp3[A_, B_] := Module[{a,c},
  a = Transpose[A];
  c = B;
  Do[c = Transpose[c.a, {2, 3, 4, 1}], {i, 4}];
  c
  ]

Now timing subroutine:

test[n_] := Module[{t0, t1, nrm, A, B,c0,c1},
  A = RandomInteger[{-10, 10}, {n, n}];
  B = RandomInteger[{-10, 10}, {n, n, n, n}];
  t0 = Timing[c0 = rp2[A, B];] // First;
  t1 = Timing[c1 = rp3[A, B];] // First;
  nrm = Norm[Flatten[c1 - c0]];
  {n, t0, t1, nrm}
  ]

Now results:

p = Table[test[n], {n, 4, 14}]
(*{{4, 0.002887, 0.000146, 0}, {5, 0.007295, 0.000156, 0}, {6, 0.023915,
   0.000331, 0}, {7, 0.065420, 0.000586, 0}, {8, 0.180312, 0.001182, 
  0}, {9, 0.437847, 0.002028, 0}, {10, 0.942596, 0.003337, 0}, {11, 
  1.941590, 0.005264, 0}, {12, 3.863945, 0.008290, 0}, {13, 7.321387, 
  0.012194, 0}, {14, 12.766735, 0.019249, 0}}*)

Now plots

ListLogPlot[{p[[All, 1 ;; 2]], p[[All, 1 ;; 3 ;; 2]]}, Joined -> True]
speedup = Transpose[{p[[All, 1]], p[[All, 2]]/p[[All, 3]]}]
(*{{4, 20.}, {5, 47.}, {6, 72.}, {7, 112.}, {8, 153.}, {9, 216.}, {10, 
  282.5}, {11, 368.8}, {12, 466.1}, {13, 600.4}, {14, 663.2}}*)
ListLogPlot[speedup, Joined -> True]

Comparison of run-times

enter image description here

Speed up

enter image description here

For sizes in the range of 4 to 14 we have a speedup in the range of 20 to 663 !

Explanation

Matrix-matrix multiplication is much faster than a set of corresponding matrix-vector multiplications!

|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ Thank you very much, I just did not see the transposition! That was the crucial point. Thanks! $\endgroup$ – Mauricio Fernández Jul 5 '16 at 16:53
  • 1
    $\begingroup$ @yarchik can you please explain {2, 3, 4, 1} in the Transpose? In particular, how to generalise this part of the function to arbitrary n and d (in the notation of the OP). Thanks! $\endgroup$ – AccidentalFourierTransform Jan 16 at 16:53
  • 2
    $\begingroup$ @AccidentalFourierTransform The idea is to use matrix multiplication, which is a very efficient operation, for the computation of each sum. To this end, we need to change the order of indices of matrix Bin such a way as to have the active one as the last one. Now concerning the generalization: d is the dimension of the inner sum in each matrix multiplication. Therefore, it is already taken care of. n is the number of sums= number of dimensions of B. To make it more general, we can write Do[c = Transpose[c.a, RotateLeft[Range[n]]], {i, n}]; $\endgroup$ – yarchik Jan 19 at 18:19
  • 1
    $\begingroup$ Also possible: Nest[A.Transpose[#, RotateLeft[Range[4]]] &, B, 4]. $\endgroup$ – Henrik Schumacher Mar 29 at 10:16
  • $\begingroup$ @HenrikSchumacher Yes, it's a very nice concise way of writing, and functional too. $\endgroup$ – yarchik Mar 29 at 15:12

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.