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;
A = RandomInteger[{-10, 10}, {n, n}];
B = RandomInteger[{-10, 10}, {n, n, n, 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?

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