getindex(A::SparseMatrixCSC, ...) creates dense column vector: out of memory

[cschwarzbach]

A[I,J] for a sparse matrix A and index vectors I and J is implemented in sparse/sparsematrix.jl. The computational core consists of three different algorithms which are implemented in Base.SparseMatrix.getindex_I_sorted_bsearch_A, Base.SparseMatrix.getindex_I_sorted_bsearch_I and Base.SparseMatrix.getindex_I_sorted_linear. A heuristics is employed to choose the best algorithm for a given problem size. The second and third algorithm internally use a dense vector of size size(A,1) as cache. For moderately sized matrices, the memory footprint is not significant. For very large sparse matrices of, e.g., size(A,1) > 10^9, the memory footprint becomes significant. At the same time, the heuristics that is used to pick the best algorithm fails. See the following code example.

mA  = 2^27
nA  = 4
nzA = 2^15 # number of non-zeros per column
A   = sprand(mA,nA,nzA/mA*nA)

println("Input: ", size(A,1), " x ", size(A,2), " sparse matrix with ", nnz(A), " non-zero entries")

nI = 2^18
I  = sort(rand(1:mA, nI))
J  = 1:nA

println("Slicing ", length(I), " x ", length(J), " submatrix")
B = A[I,J]
println("Output: ", size(B,1), " x ", size(B,2), " sparse matrix with ", nnz(B), " non-zero entries")
println()

f = [ Base.SparseMatrix.getindex_I_sorted
      Base.SparseMatrix.getindex_I_sorted_bsearch_A
      Base.SparseMatrix.getindex_I_sorted_bsearch_I
      Base.SparseMatrix.getindex_I_sorted_linear ]

for k = 1:length(f)
    println("B = ", f[k], "(A,I,J)")
    @time B = f[k](A,I,J)
    @time B = f[k](A,I,J)
    @time B = f[k](A,I,J)
    println()
end

The generated screen output looks as follows:

Input: 134217728 x 4 sparse matrix with 525820 non-zero entries
Slicing 262144 x 4 submatrix
Output: 262144 x 4 sparse matrix with 1010 non-zero entries

B = getindex_I_sorted(A,I,J)
  0.914429 seconds (6 allocations: 1.000 GB, 9.74% gc time)
  0.926163 seconds (6 allocations: 1.000 GB, 6.86% gc time)
  0.754537 seconds (6 allocations: 1.000 GB, 8.42% gc time)

B = getindex_I_sorted_bsearch_A(A,I,J)
  0.023197 seconds (4 allocations: 16.094 KB)
  0.023501 seconds (4 allocations: 16.094 KB)
  0.025746 seconds (4 allocations: 16.094 KB)

B = getindex_I_sorted_bsearch_I(A,I,J)
  0.728229 seconds (6 allocations: 1.000 GB, 9.55% gc time)
  0.701622 seconds (6 allocations: 1.000 GB, 8.67% gc time)
  0.725397 seconds (6 allocations: 1.000 GB, 8.03% gc time)

B = getindex_I_sorted_linear(A,I,J)
  0.453750 seconds (6 allocations: 1.000 GB, 13.28% gc time)
  0.454082 seconds (6 allocations: 1.000 GB, 13.58% gc time)
  0.477592 seconds (6 allocations: 1.000 GB, 13.85% gc time)

This is with a very recent Julia 0.4:

Julia Version 0.4.0-dev+7002
Commit f42b222* (2015-08-26 20:27 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-4980HQ CPU @ 2.80GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

Note that the heuristic choice (getindex_I_sorted(A,I,J) equivalent to A[I,J]) allocates 1 GB memory and has the longest runtime of the three algorithms.

I think that the conversion of even just a column of a sparse matrix to a dense vector should be avoided.

[jiahao]

Please see #11324 - support for sparse vectors is planned for 0.5.

[eldadHaber]

The problem is not only for sparse vectors. It's for all sparse matrices

In principle you should not convert a sparse array into a dense one. We as users take our data types seriously and do not like unexpected conversions.

[IainNZ]

In my opinion, the conversion of even just a column of a sparse matrix defies the purpose of sparse storage. As a user, I carefully pick my tools and select sparse storage because I cannot afford dense storage. I expect the same from the algorithms operating on my sparse object. In fact, the internal conversion to a dense column vector has caused my actual application to run out of memory and it took me a lot of time and effort to track the problem down to getindex.

@cschwarzbach I'm not sure if you are new to participating in open source communities, but the tone and the way this paragraph was written is not necessary. It would have sufficed to simply say "I believe take dense slices shouldn't be produced from sparse matrices", and did not require you explaining how careful you are when using selecting a data structure or how much difficulty this implementation weakness in a free, open-source software caused you - it can be read as insulting to people who have done the work that exists so far that you are otherwise benefiting from.

@eldadHaber: @jiahao understands that very well, he was referring to that issue because ideally slicing a column from a sparse matrix would give you a sparse vector. Same comment goes for your tone as for @cschwarzbach - there doesn't really exist such a clear user-developer distinction here, you are not paying for a product, but participating in the development of a pre-release open source software project.

[eldadHaber]

To @IainNZ and @jiahao

No intention to be insulting and I appologize if it seems so. We certainly understand the effort that is invested here.
We are "heavy" users that develop a large software library using Julia and invest a lot of effort in developing and then teaching and promoting the language. We want and need this effort to be successful. Your success is ours as well

So please, can we get back to the point.
I would suggest to fix this and not convert sparse to dense.

[nalimilan]

So please, can we get back to the point.

I would suggest to fix this and not convert sparse to dense.

I think what @jiahao and @IainNZ meant was that the point is well taken. We're just waiting for sparse vectors support to be merged into Base to switch getindex to returning this new structure instead of a dense vector. That should happen for the next 0.5 release (it's too late to include this breaking change in 0.4), but if you really need it right now you could probably code something using the SparseVectors.jl package.

[timholy]

I would suggest to fix this and not convert sparse to dense.

Just a reminder, it's best if the people who most care about an issue implement the fix themselves; you guys obviously care a lot about this. (That doesn't mean that some hero won't step up to the plate and do your work for you, but you can never count on when/if that will happen.) Want to give it a try? It might not be hard. See https://github.com/JuliaLang/julia/blob/master/CONTRIBUTING.md

[ViralBShah]

@cschwarzbach Can you suggest a heuristic that will work for you?

[ViralBShah]

Cc: @tanmaykm

[cschwarzbach]

@IainNZ I'm sorry if I hurt anybody's feelings; I had no intention to do so. I've deleted the offending section in my post. I appreciate the hard work that the Julia community has done to pull this amazing project together.

[eldadHaber]

If you can put your fix it will be better.

On Aug 29, 2015, at 11:42 AM, Christoph Schwarzbach [email protected] wrote:

@IainNZ I'm sorry if I hurt anybody's feelings; I had no intention to do so. I've deleted the offending section in my post. I appreciate the hard work that the Julia community has done to pull this amazing project together.

—
Reply to this email directly or view it on GitHub.

[cschwarzbach]

@ViralBShah @timholy

Of the three algorithms for getindex that are implemented for sparse matrices, only Base.SparseMatrix.getindex_I_sorted_bsearch_A does not create a dense vector of size size(A,1). Therefore, I would just default to this implementation and discontinue using the other two which create a dense vector as a cache. This is my preferred choice. However, I'm aware of the fact that the other two algorithms are there for a reason. I don't quite know what impact my proposed change will have on the runtime tests that have gone into determining the existing heuristics and if the potential increase in runtime will be acceptable by the community.

Of course, one could come up with a more elaborate heuristic which considers dense vectors up to a certain size acceptable. However, this seems to me very problem and hardware dependent. I'm afraid of opening Pandora's box.

[eldadHaber]

I also suggest to stick to the sparse algorithm. I would add in the help file that if the size of the index set is large then the user should density her array before getting the index. This will give the choice in the hand of the user.

[ViralBShah]

I would really like to tweak the heuristics, and not slow down things for people who benefit from this. Clearly, we just need to put a ceiling on the largest dense vector we create, for which we can figure out a good threshold.

[mauro3]

The relevant performance test are here: https://github.com/JuliaLang/julia/blob/349a4e197728d010d107472d44ba9ccece0876f7/test/perf/sparse/getindex.jl

It would be great to add above test case to it (well, maybe with smaller size).

At least my sparse matrices (and I suspect many others) are pretty square. Then allocating a dense vector of size(A,1) should have a memory footprint a lot smaller than the matrix itself. So, I can imagine that in some situations the allocation is worth it, and presumably that is what was found when implementing the heuristics.