Efficient implementation of tensorial Rayleigh product

Question

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

Answer 1

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

Speed up

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!