19

The goal is to assemble a SparseArray in an additive fashion.

Let us assume we have a large List of indices (some will be repeated). We generate a simple test list of indices and values

ind = RandomInteger[{1, 4}, {10, 2}];
val = RandomReal[{-1,1}, Length[ind]];

where each value corresponds to an index from ind. I would like to build a SparseArray in a way such that the repeated index values are summed into the array.

If we simply use:

SparseArray[ind -> val, {4,4}]

only the first index encounter is written into the SparseArray, all repeated indices are ignored.

Current Solution (slow + ugly)

This solution is slow and is only shown to make precise what exactly I am trying to accomplish. We pre-allocate a sparse array of the correct size and use Do to accumulate the values at each index:

n = 5;
ind = RandomInteger[{1, n}, {3*n, 2}];
val = RandomReal[{1, 1}, Length[ind]];
A = SparseArray[{1, 1} -> 0, {n, n}];
Do[
   A[[Sequence @@ ind[[i]]]] += val[[i]]
,{i, 1, Length[val]}
]

There are some great tips on working with SparseArrays in Efficient by-element updates to SparseArrays and SparseArray row operations. A clever combination of GatherBy, Sort, etc. operations on ind and val may be good path to head down. I just can't see it yet.

Community
  • 1
leibs
  • 758
  • 3
  • 11

2 Answers2

21

There is actually an undocumented System Option that tells Mathematica to do this automatically. The default behavior:

ind = {{3, 1}, {3, 3}, {1, 3}, {2, 1}, {3, 2}, {3, 1}, {3, 2}, {3, 3}, {1, 3}, {3, 1}};
val = {1, 1, 3, 0, 3, 4, 3, 1, 1, 1};

SparseArray[ind -> val] // Grid

$ \begin{matrix} 0 & 0 & 3 \\ 0 & 0 & 0 \\ 1 & 3 & 1 \end{matrix} $

Now with our "magic" Option (learned from Oliver Ruebenkoenig):

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
(* equivalently:
   SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}] *)
SparseArray[ind -> val] // Grid

$ \begin{matrix} 0 & 0 & 4 \\ 0 & 0 & 0 \\ 6 & 6 & 2 \end{matrix} $

Arbitrary functions are accepted by version 9 and later; see: Optimising 2D binning code

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> (# - +##2 &)}];

SparseArray[ind -> val] // Grid

$\begin{matrix} 0 & 0 & 2 \\ 0 & 0 & 0 \\ -4 & 0 & 0 \\ \end{matrix}$

To restore the default behavior set:

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
(* equivalently:
   SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}] *)

To encapsulate this setting, use Internal`WithLocalSettings: SetOptions locally?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
  • 4
    Another reference for "TreatRepeatedEntries": What are some useful, undocumented Mathematica functions? – Michael E2 Nov 22 '13 at 13:17
  • @MichaelE2 What a gem. I never would have found that. @Mr.Wizard Any hidden SparseArray tricks for setting all values in a row to zero except the diagonal, which should be set to one. – leibs Nov 22 '13 at 18:11
  • 2
    As far as zeroing a row, I found that I can just take the SparseArray, A and indices of zero rows rows and set A[[rows,All]] = 0.0'. Adding1back to the diagonal of the deleted rows can be done withA = A + SparseArray[Thread[{#, #} &[rows]] -> ConstantArray[1.0, Length[rows]],{nDOF, nDOF}, 0.0];WherenDOFis the dimension ofA`. – leibs Nov 22 '13 at 18:48
7
n=4;    
ind = RandomInteger[{1, n}, {10, 2}]
val = RandomReal[{-1, 1}, Length[ind]]

#[[1, 1]] -> Total[#[[;; , 2]]] & /@ GatherBy[Thread[{ind, val}], First]

a = SparseArray[%, {n, n}]

or quite ugly but compact:

b = SparseArray[#, {n, n}] & /@ Thread[ind -> val] // Total

a == b

True
Kuba
  • 136,707
  • 13
  • 279
  • 740
  • 1
    Also possible :(#[[1, 1]] -> Total[#[[;; , 2]]]) & /@ GatherBy[Thread[{ind, val}], First] – andre314 Nov 21 '13 at 21:40
  • @andre good point, should be faster. – Kuba Nov 21 '13 at 21:40
  • Great solutions. The second method SparseArray[#, {n, n}] & /@ Thread[ind -> val] // Total seems to really hang. The first listed and the solution by @andre are both very fast. THANKS! – leibs Nov 21 '13 at 21:49