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}] [![enter image description here][1]][1] As expected `sum` has 25 terms. We get the same result using [`Total`](https://reference.wolfram.com/language/ref/Total.html): 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, [1]: https://i.sstatic.net/qRR3Y.png