Memory leak when creating a SparseArray

Question

Problem: I have a list e={u,v,w} and wish to construct a sparse matrix with entries {u[i],v[i]}->w[i]. Note that the same position can appear several times (duplicate entries), and their values should be summed. For testing purposes, an example:

m=n=10^6; r=5*10^8; 
e={RandomInteger[{1,m},r],RandomInteger[{1,n},r],RandomInteger[{-9,9},r]}; 
ByteCount[e]

12000000496

Solution:

AbsoluteTiming@MaxMemoryUsed[ 
e=Transpose[e[[1;;2]]]->e[[3]]; b=SparseArray[e,{m,n}]; ]

{129.894, 18961378648}

This is fairly efficient, but gives the wrong result, since repeated entries are discarded instead of summed. If I first run the code SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries"->1}]; (I learned about it in this post) and then the above solution, I get the correct result, but also a memory leak:

{153.46, 53065409208}

53GB from 12GB is too much : (. I would really appreciate any suggestions on how to improve this.

Answer 1

Not an answer but too long for a comment.

Well, this additional memory consumption is actually not as bad as you might think, and I am not so sure whether this should be considered a memory leak. The additive assembler has

1. to sort the entries of Transpose[e[[1;;2]]];

2. to reorder the entries of e[[3]];

3. to compute the arrays for the row pointers;

4. to compute the arrays for the column indices; and

5. to accumulate the values from (the reordered) array e[[3]] into the list of nonzeroes.

I guess one would fuse steps 3. to 5. together. Then one has to hold at the same time: arrays of the Dimensions {r,2} (for e[[1;;2]]), {r} (for e[[3]]), {m+1} (for the row pointers), {nnz,1} (for the column indices) and {nnz} (for the list of nonzero entries). Here nnz is the number of nonzeroes in the final matrix. Since this is not known in the beginning, probable r is used as upper bound so that memory of {r,1} (for the column indices) and {r} (for the list of nonzero entries) is allocated for the computations (and the tails are discarded when computations are done).

Also the sorting needs some scratch space - and here one can basically trade runtime for memory consumption (e.g. bucket sort has complexity $O(r)$ but requires quite a lot of memory while quicksort is in-place but has complexity $O(r \, \log(r))$. I don't know which particular sorting algorithm is used here. The memory consumption for the sorting should actually not add up to the one for the rest of the assembly. But maybe the scratch space for the sort was not deallocated directly after the sort but only later? (I find it unlikely, though.)

All this does probably not account for all the memory insufficiency, but for quite a lot of it. Further memory consumption may be due to the behavior of LibraryLink: Typically it has to form a copy for returning the results. But that pretty much depends on the implementation details into which I have to insight.

Moreover, some parts of the assembly might be parallelized, e.g. with OpenMP. Then each thread probably needs some further scratch space to avoid write locks into shared memory. And that may also blow up the memory consumption dramatically.

Anyways, all this suggetes that one should try to figure out how to compute the column indices and row pointers directly and to generate the SparseArray with the undocumented syntax

SparseArray @@ {Automatic, {m, n}, 0, {1, {rp, ci}, vals}}

SparseArray @@ {Automatic, {m, n}, 0., {1, {rp, ci}, vals}}

SparseArray @@ {Automatic, {m, n}, 0.+0.I, {1, {rp, ci}, vals}}

(depending on which scalar type you want), where rp corresponds to b["RowPointers"], ci corresponds to b["ColumnIndices"], and vals corresponds to b["NonzeroValues"], provided that matrix b is in proper CRS format (i.e., if SparseArray`SparseArraySortedQ[b] returns True).

Finally, I'd like to say that one should be careful when interpreting the result of MaxMemoryUsed. Actually, the documentation is super vague about what it does and the sentence

On most computer systems, MaxMemoryUsed[] will give results close to those obtained from external process status requests.

does not really raise any confidence.