CODEWITHSUNDEEP
juliasparse-matrix
Shep Bryan
Shep Bryan

Reputation: 669

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?

Upvotes: 1

Views: 1148

Answers (1)

Benoit Pasquier
Benoit Pasquier

Reputation: 2995

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.

Upvotes: 1

Related Questions