2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – KraZug Feb 13 '18 at 10:42
3
$\begingroup$

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

Addendum

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]

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

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – KraZug Feb 13 '18 at 10:17

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.