Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
I just was about to post a new question asking this same thing when I came across your question. Thanks, leibs and Mr.Wizard! –  DumpsterDoofus Dec 24 '13 at 23:29
add comment

2 Answers

up vote 7 down vote accepted

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{array}{ccc} 0 & 0 & 3 \\ 0 & 0 & 0 \\ 1 & 3 & 1 \end{array} $

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

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];

SparseArray[ind -> val] // Grid

$ \begin{array}{ccc} 0 & 0 & 4 \\ 0 & 0 & 0 \\ 6 & 6 & 2 \end{array} $

To restore the default behavior set:

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];

Additional option values are accepted by Version 9; see: Optimising 2D binning code

share|improve this answer
2  
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
    
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'. Adding 1` back to the diagonal of the deleted rows can be done with A = A + SparseArray[Thread[{#, #} &[rows]] -> ConstantArray[1.0, Length[rows]],{nDOF, nDOF}, 0.0]; Where nDOF is the dimension of A. –  leibs Nov 22 '13 at 18:48
add comment
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
share|improve this answer
1  
Also possible :(#[[1, 1]] -> Total[#[[;; , 2]]]) & /@ GatherBy[Thread[{ind, val}], First] –  andre 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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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