Search code examples
juliasparse-matrix

Best way to fill a sparse matrix


What is the most efficient way to fill a sparse matrix? I know that sparse matrixes are CSC, so I expected it to be fast to fill them column by column like

using SparseArrays

M = 100
N = 1000

sparr = spzeros(Float64, M, N)

for n = 1:N

  # do some math here

  idx = <<boolean index array of nonzero elements in the nth column>>
  vals =  <<values of nonzero elements in the nth column>>

  sparr[idx, n] = vals

end

However, I find that this scales very poorly with N. Is there a better way to fill the array? Or perhaps, I should not bother with filling the array and instead initialize the matrix differently?


Solution

  • You can do sparse(I, J, V, M, N) directly:

    julia> using SparseArrays
    
    julia> M = 100;
    
    julia> N = 1000;
    
    julia> nz = 2000; # number of nonzeros
    
    julia> I = rand(1:M, nz); # dummy I indices
    
    julia> J = rand(1:N, nz); # dummy J indices
    
    julia> V = randn(nz); # dummy matrix values
    
    julia> sparse(I, J, V, M, N)
    100×1000 SparseMatrixCSC{Float64, Int64} with 1982 stored entries:
    ⣻⣿⣿⣿⣿⡿⣾⣿⣿⣿⣿⣿⣿⣷⣾⣽⣿⢿⢿⣿⣿⣿⢿⣿⣾⣿⣽⣿⣿⣾⣿⣿⣿⣿⣿⣿⣿⣿⣻⣿
    ⣼⣿⣿⡿⣿⣿⡽⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣻⣿⡿⣿⣿⣿⡿⣿⡿⣯⢿⣿⠾⣿⣿⡿⢿⣿⣻⡿⣾
    

    which should scale decently with size. For more expert use, you could directly construct the SparseMatrixCSC object.


    EDIT:

    Note that if you have to stick with the pseudo code you gave with a for-loop for column indices and values, you could simply concatenate them and create I, J, and V:

    I = Int[]
    J = Int[]
    V = Float64[]
    
    for n = 1:N
    
      # do some math here
    
      idx = <<boolean index array of nonzero elements in the nth column>>
      vals =  <<values of nonzero elements in the nth column>>
    
      I = [I; idx]
      J = [J; fill(n, length(I))]
      V = [V; vals]
    
    end
    

    but that'd be slower I think.