# Nested sum in replacement rule of unknown size

I have a function which takes a list as its argument, and I want to transform it by multiplying by a matrix for each element of a list. To be specific, I want the following transformation: $$\hat{f}_i = \sum_p^n A_{ip}f_p, \quad \hat{f}_{ij} =\sum_i^n \sum_j^n A_{ip}A_{jq} f_{pq}, \quad \hat{f}_{ijk} =\sum_i^n \sum_j^n \sum_k^n A_{ip} A_{jq} A_{kr} f_{pqr}$$ and so on.

I can write some replacement rules to do this explicitly for each size of the list,

n = 4;
reps = {f[{i_}] :> Sum[A[i, p] f[{p}], {p, 1, n}],
f[{i_, j_}] :> Sum[A[i, p] A[j, q] f[{p, q}], {p, 1, n}, {q, 1, n}],
f[{i_, j_, k_}] :> Sum[A[i, p] A[j, q] A[k, r] f[{p, q, r}], {p,1,n}, {q,1,n}, {r,1,n}]}


For the first order, I get

f[{1}] /. reps
(* A[1, 1] f[{1}] + A[1, 2] f[{2}] + A[1, 3] f[{3}] + A[1, 4] f[{4}] *)


with increasingly large expressions for the longer lists. The large expressions are fine, and speed is not a concern. I want to manipulate them symbolically, without assuming values for the matrix A at this point.

I need a solution that takes any size list (up to say length 12), and figure there must be a better way than explicitly listing them out.

• Note that speed is not my primary goal here, and I don't have an explicit matrix of numbers, I want to do this symbolically. Feb 13, 2018 at 10:42

Using TensorContract with TensorProduct won't be as fast as Compile or Dot solutions, but might be more understandable:

transform[A_][tensor_] := With[{rank = TensorRank[tensor]},
Activate @ TensorContract[
Inactive[TensorProduct] @@ Append[ConstantArray[A, rank], tensor],
Table[{2 i, 2 rank+i}, {i, rank}]
]
]


(see this answer for an explanation of the use of Inactive here)

Here are examples for ranks 1, 2, and 3:

b1 = transform[A][f1]; //RepeatedTiming
b2 = transform[A][f2]; //RepeatedTiming
b3 = transform[A][f3]; //RepeatedTiming


{0.000080, Null}

{0.00011, Null}

{0.00030, Null}

and a check with @HenrikSchumacher's answer:

a1 == b1
a2 == b2
Block[{Internal\$EqualTolerance = 4}, a3 == b3]


True

True

True

Based on your comment, perhaps the following variation of my answer is more like what you're looking for:

transform[A_, t_, indices_] := TensorContract[
TensorProduct @@ Append[A[[#]]& /@ indices, t],
Table[{i, i+Length[indices]}, {i, Length[indices]}]
]


Examples:

transform[Array[A, {4,4}], Array[f, 4], {1}]
transform[Array[A, {4,4}], Array[f, {4,4}], {1,3}]


A[1, 1] f[1] + A[1, 2] f[2] + A[1, 3] f[3] + A[1, 4] f[4]

A[1, 1] A[3, 1] f[1, 1] + A[1, 1] A[3, 2] f[1, 2] + A[1, 1] A[3, 3] f[1, 3] + A[1, 1] A[3, 4] f[1, 4] + A[1, 2] A[3, 1] f[2, 1] + A[1, 2] A[3, 2] f[2, 2] + A[1, 2] A[3, 3] f[2, 3] + A[1, 2] A[3, 4] f[2, 4] + A[1, 3] A[3, 1] f[3, 1] + A[1, 3] A[3, 2] f[3, 2] + A[1, 3] A[3, 3] f[3, 3] + A[1, 3] A[3, 4] f[3, 4] + A[1, 4] A[3, 1] f[4, 1] + A[1, 4] A[3, 2] f[4, 2] + A[1, 4] A[3, 3] f[4, 3] + A[1, 4] A[3, 4] f[4, 4]

• I do not have an explicit tensor, I just have a list of indices. I want to manipulate it symbolically before giving any numbers. Feb 13, 2018 at 10:20
• The addendum version does what I need, thanks. Mar 16, 2018 at 14:47

Better store the f and A as vectors and matrices and use Dot. That's easier and much, much faster.

n = 10;
A = RandomReal[{-1, 1}, {n, n}];
f1 = RandomReal[{-1, 1}, {n}];
f2 = RandomReal[{-1, 1}, {n, n}];
f3 = RandomReal[{-1, 1}, {n, n, n}];

a1 = A.f1; // AbsoluteTiming
b1 = Table[
Sum[A[[i, p]] f1[[p]], {p, 1, n}],
{i, 1, n}]; // AbsoluteTiming
Max[Abs[a1 - b1]]


{8.*10^-6, Null}

{0.000152, Null}

0.

a2 = A.f2.A\[Transpose]; // AbsoluteTiming
b2 = Table[
Sum[A[[i, p]] A[[j, q]] f2[[p, q]], {p, 1, n}, {q, 1, n}],
{i, 1, n}, {j, 1, n}]; // AbsoluteTiming
Max[Abs[a2 - b2]]


{0.000014, Null}

{0.019346, Null}

1.33227*10^-15

The third one is a bit more involved:

b3 = Table[
Sum[A[[i, p]] A[[j, q]] A[[k, r]] f3[[p, q, r]], {p, 1, n}, {q, 1,
n}, {r, 1, n}],
{i, 1, n}, {j, 1, n}, {k, 1, n}]; // AbsoluteTiming
a3 = Compile[{{A, _Real, 2}, {B, _Real, 2}}, A.B,
RuntimeAttributes -> {Listable},
Parallelization -> True
][A, A.f3.A\[Transpose]]; // AbsoluteTiming
Max[Abs[a3 - b3]]
`

{2.66485, Null}

{0.000152, Null}

5.32907*10^-15

• I want a replacement rule that will work for any size list. I can't see how to generalise what you have written to do that. Feb 13, 2018 at 10:17