Performing sparse sum on Mathematica

Question

I want to evaluate a sum in Mathematica of the form

Sum[g[[i,j,k,l,m,n]] x g[[o,p,q,r,s,t]] x (complicated function of the indices i,j etc), {i,0,3}, {j,0,3}, {k,0,3}, {l,0,3}, {m,0,3}, {n,0,3}, {o,0,3}, {p,0,3}, {q,0,3}, {r,0,3}, {s,0,3}, {t,0,3}]

But all these indices range from 0 to 3, so the total number of cases to sum over is 4^12, which will take an unforgiving amount of time. However, barely any elements of the array g[[i,j,k,l,m,n]] are nonzero -- there are probably around 8 nonzero entries -- so I would like to restrict the sum over {i,j,k,l,m,n,o,p,q,r,s,t} to precisely those combinations of indices for which both factors of g are nonzero.

I can't find a way to do this for summation over multiple indices, where the allowed index choices are particular combinations of {i,j,k,l,m,n} as opposed to specific values of each particular index. Any help appreciated!

tried SparseArray and the property "NonzeroPositions"? See also: What are SparseArray Properties? How and when should they be used? — kglr, Mar 04 '22 at 23:04
@kglr those properties are now documented and this one is referred to as "ExplicitPositions". https://reference.wolfram.com/language/ref/SparseArray.html#1948657 — Greg Hurst, Mar 04 '22 at 23:37
Thank you very much @GregHurst (didn't know about the v13.0 update). — kglr, Mar 04 '22 at 23:42
@kglr How does that help me evaluate this sum? I understand that I can use NonzeroPositions to determine which elements of my array g are nonzero, but how do I then restrict the sum to only these values of the indices? — pseudo spin, Mar 04 '22 at 23:43
pseudo spin, by g[[i,j,k,l,m,n]] x g[[o,p,q,r,s,t]] x ... did you mean g[[i,j,k,l,m,n]] + g[[o,p,q,r,s,t]] + ...? — kglr, Mar 04 '22 at 23:45
@kglr No, terribly sorry -- I meant that I have a product of three quantities, which are functions of 12 indices, and I want to sum over those indices. I have edited my question to make this clearer!! — pseudo spin, Mar 04 '22 at 23:47
@pseudospin is there a condition that determines which elements are zero based on the symmetries of the indices? it was not clear to me from the OP that's why I am asking — , Mar 04 '22 at 23:55
one possible way: check if this gives what you need: g = SparseArray[g]; nzp = s["NonzeroPositions"]; Sum[s[[## & @@ i]] s[[## & @@ j]] FOO[i, j], {i, nzp}, {j, nzp}] — kglr, Mar 04 '22 at 23:57
@kcr Not as far as I can tell -- there possibly is but it would be hard to figure it out — pseudo spin, Mar 04 '22 at 23:59
Total[s[[## & @@ #]] s[[## & @@ #2]] FOO[##] & @@@ Tuples[nzp, 2]] should give the same result. — kglr, Mar 05 '22 at 00:06
@kglr Using this on a test case, I think this works fine! Do you want to post this as an answer so I can accept it? — pseudo spin, Mar 05 '22 at 00:29

kglr · Accepted Answer · 2022-03-05T16:49:57.357

complicatedFunctionOfIndices = FOO[FromDigits @ Sort @ #, FromDigits @ ReverseSort @#2]&

SeedRandom[1]
myg = RandomChoice[{3000, 1, 1, 1, 1, 1} -> Range[0, 5], {4, 4, 4, 4, 4, 4}];

Dimensions[myg]

{4, 4, 4, 4, 4, 4}

sa = SparseArray[myg]

nonZeroPositions = sa["NonzeroPositions"]

{{1, 1, 3, 3, 4, 3}, {1, 4, 3, 3, 4, 4}, {2, 2, 1, 2, 2, 4}, 
 {4, 2, 1, 2, 1, 4}, {4, 2, 2, 1, 3, 2}}

nonZeroValues = sa["NonzeroValues"]

 {3, 5, 3, 1, 4}

sum = Sum[sa[[## & @@ i]] sa[[## & @@ j]] complicatedFunctionOfIndices[i, j], 
  {i, nonZeroPositions}, {j, nonZeroPositions}]

As expected sum has 25 terms.

We get the same result using Total:

total = Total[sa[[## & @@ #]] sa[[## & @@ #2]] complicatedFunctionOfIndices[##]  & @@@ 
  Tuples[nonZeroPositions, 2]];
sum == total

True

Note: Per Greg Hurt's comment above, replace "NonzeroPositions" with "ExplicitPositions" if you have version 13.0,

Performing sparse sum on Mathematica

1 Answers1