# Efficient Linear Operators on SparseArrays

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.

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...

• 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. – VolatileStorm Sep 19 '19 at 8:06
• 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. – Henrik Schumacher Sep 19 '19 at 14:04