SparseMatrix is very slow

General Usage
1

Hi I have the following simple code:

import SparseArrays

function func0(K::Array{Float64,2}, M::Array{Float64,2}, r::Array{Float64,2}, x::Array{Float64,2}, row::Array{Int32,1}, col::Array{Int32,1}, data::Array{Float64,1}, m::Int64, n::Int64)
    c = SparseArrays.sparse(row, col, data, m, n)
    u = 1.0 ./ x    # dense
    ku = (1.0 ./ (transpose(K) * u))  # dense
    v = c .* ku   # sparse * dense element-wise product (c is sparse and ku is dense), very very slow
end

Somehow line v = c .* ku is super slow, 10x slower than when c is represented as dense. And c is very sparse. Am I missing anything? Thanks!

2

My guess is that the broadcast machinery here doesn’t know that the element-wise product of sparse and dense is sparse, so is looping over everything (not the nonzeroes) of the sparse matrix. I think your best bet will be re-writing this as a loop over non-zero elements of c.

1 Like
3

I see. Thanks a lot for the suggestion! But the slowdown seems beyond the overhead of looping all elements in the sparse matrix. I tried to conver sparse matrix c to dense, and then do the element-wise product, and that was much faster although it also loops over all elements.

Is it possible that the sparse matrix is CSC format and the dense array is in row-major order, so the iteration is causing a lot of cache misses?

Thanks!

1 Like
4

Iterating over sparse arrays is much more expensive than iterating over dense arrays.
Dense arrays in Julia are also column major.

5

Oh oh I see. Thanks!

6

A good iteration scheme when using sparse matrices is the following (from the documentation: Sparse Arrays · The Julia Language)

A = sparse(I,J,V)
rows = rowvals(A)
vals = nonzeros(A)
m, n = size(A)
for j = 1:n
   for i in nzrange(A, j)
      row = rows[i]
      val = vals[i]
      # perform sparse wizardry...
   end
end

You could probably adapt it to do your multiplication.

And welcome to the Julia community! :slight_smile:

3 Likes
7

Thanks a lot! That example is very helpful!!

Btw, are such methods on SparseMatrix documented somewhere? like rowvals(A) and nonzeros(A)?

I tried to do methodswith(SparseMatrixCSC), but it doesn’t have any comments for a method.

Thanks!

8

I would start reading here Sparse Arrays · The Julia Language if I were you.

3 Likes
9

yes sure, just scroll down in the link I sent you. For example, rowvals would be here.

What you can also do: In the Julia REPL enter “?” to get in the help mode. Then you can also enter e.g. “rowvals” to get documentation info:

julia> ?
help?> rowvals
search: rowvals

  rowvals(A::SparseMatrixCSC)

  Return a vector of the row indices of A. Any modifications to the returned vector will mutate A as well. Providing access to how the row indices are stored internally can be useful in conjunction with iterating over structural nonzero values. See also nonzeros and nzrange.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> A = sparse(2I, 3, 3)
  3Ă—3 SparseMatrixCSC{Int64,Int64} with 3 stored entries:
    [1, 1]  =  2
    [2, 2]  =  2
    [3, 3]  =  2
  
  julia> rowvals(A)
  3-element Array{Int64,1}:
   1
   2
   3
10

I see. Thanks a lot for the detailed example! Appreciate it!