# Efficient implementation of tensorial Rayleigh product

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];][]


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];][]
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];][]
AbsoluteTiming[res2 = rp2[A, B];][]
res1 == res2


0.00277906

0.770904

True

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!

• Thank you very much, I just did not see the transposition! That was the crucial point. Thanks! – Mauricio Fernández Jul 5 '16 at 16:53
• @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! – AccidentalFourierTransform Jan 16 at 16:53
• @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}]; – yarchik Jan 19 at 18:19
• Also possible: Nest[A.Transpose[#, RotateLeft[Range]] &, B, 4]. – Henrik Schumacher Mar 29 at 10:16
• @HenrikSchumacher Yes, it's a very nice concise way of writing, and functional too. – yarchik Mar 29 at 15:12