2
$\begingroup$

I have a rank four SparseArray of dimension {40, 40, 40, 40}. I wish to apply a linear operator to it, whose form I know abstractly. That is, I have a simple closed form for the action of the operator on each basis vector. Currently, I am implementing the action of this operator as:

operator[s_SparseArray] := DeleteCases[ArrayRules[s], HoldPattern[{_, _} -> 0]] // Map[operator] // Total // SparseArray;

operator[{a_, b_, c_, d_} -> e_] := 
    Total[SparseArray[#, {40,40,40,40}] & /@ {
         {a, b, c, d} -> e/4, 
         {b, a, c, d} -> -e/4, 
         {b, c, a, d} -> -e/4, 
         {b, c, d, a} -> -e/4, 
         {c, b, a, d} -> e/4, 
         {c, b, d, a} -> e/4, 
         {c, d, b, a} -> e/4, 
         {d, c, b, a} -> -e/4
     }];

This approach, with a matrix with 73695 non-zero elements is taking ~ 6 seconds on my laptop. Are there any better ideas for implementing something like this? I do this thousands of times, so any order of magnitude games would be important.

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
3
$\begingroup$

Using Sum and Transpose is about 50 times faster on my machine:

n = 40;
m = 73695;

weights = 0.25 {1, -1, -1, -1, 1, 1, 1, -1};
perms = InversePermutation /@ {{1, 2, 3, 4}, {2, 1, 3, 4}, {2, 3, 1,4}, {2, 3, 4, 1}, {3, 2, 1, 4}, {3, 2, 4, 1}, {3, 4, 2, 1}, {4, 3, 2, 1}};

A = SparseArray[RandomInteger[{1, n}, {m, 4}] -> RandomReal[{-1, 1}, m], {n, n, n, n}];
B = operator[A]; // AbsoluteTiming // First
B2 = Sum[weights[[i]] Transpose[A, perms[[i]]], {i, 1, Length[perms]}]; // AbsoluteTiming // First
Max[Abs[B - B2]]

4.73512

0.097895

0.

Moreover, one could probably safe quite an amount of time if one knew that all the input matrices share the same sparsity pattern...

Improve this answer
$\endgroup$
2
  • $\begingroup$ There are some nice tricks in only a few lines of code! I was unaware of how Transpose can be used to rearrange the "axes" of higher rank tensors. I questioned the use of Sum and explicit indexing, but naive use of Map/Dot breaks it out of SparseArray format. Note that in the question I am using exact arithmetic, and this answer uses floating point arithmetic. This change alone is responsible for a factor of three in the speedup. I hadn't realized the importance, but my input arrays are rational valued. $\endgroup$
    – V.S.
    Commented Sep 19, 2019 at 8:06
  • $\begingroup$ Indeed, I usually try to avoid Sum and to replace it by Dot in some way. So I was a bit surprised that the Sum method worked that well. Originally, I tried something really fancy with A["NonzeroValues"] and A["NonzeroPositions"] and to reassemble the tensor again. But that was (only slightly) slower than Sum+Transpose. As a general advice: Avoid to use ArrayRules[A] because it unpacks arrays; use A["NonzeroValues"] and A["NonzeroPositions"] instead. $\endgroup$
    – Henrik Schumacher
    Commented Sep 19, 2019 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.