Convenience type for assembling sparse matrices

Internals & Design
1

I find the IJV-approach to assembling sparse matrices quite cumbersome as it forces me to write three lines to set a single entry A[i,j] = v:

push!(I,i)
push!(J,j)
push!(V,v)

The situation becomes even worse if instead of a single entry I want to add a whole block of entries A[i1:i2, j1:j2] = v since then the entries of i1:i2 and j1:j2 have to be repeated.

I think it would be nice to have a convenience SparseMatrixAssembler type which wraps the IJV vectors and allows to use setindex! to achieve the above. Once you are done assembling, you simply call SparseMatrixCSC(assembler[,m,n]) to convert to a proper sparse matrix. For example:

julia> A = SparseMatrixAssembler();

julia> A[2,2] = π;

julia> SparseMatrixCSC(A,3,3)
3Ă—3 SparseMatrixCSC{Float64,Int64} with 1 stored entry:
  [2, 2]  =  3.14159

A rough prototype implementation can be found at GitHub - ettersi/SparseMatrixAssemblers.jl: Assemble sparse matrices using block assignment. Any chance a more carefully designed version of this could make it into the standard SparseArrays package?

4 Likes
2

I’m replying just to say that I have the same issue in 2023. Scipy has the coo format which serves the same purpose; Julia is lagging behind here.

Another advantage of this solution is that it would make easier to switch between assembling a dense and a sparse matrix, since the same syntax would work in both cases.

3 Likes
3

a) I think SparseArrayInitializer or something of that nature would be more consistent with existing names
b) Is there any particular reason to restrict this to matrices? If you built it to be generalizable to other-dimensional-arrays, how would you do it differently?

4

Unlike SciPy, there’s no reason for this to be built in, and adding more baggage to the SparseArray stdlib doesn’t make sense, so if you’re waiting for this to be solved in base Julia, it won’t. If you need such a type why not use the implementation above or something based on it?

5

An oldie: FinEtools.jl/AssemblyModule.jl at main · PetrKryslUCSD/FinEtools.jl · GitHub

And a “more Julian” solution Elfel.jl/Assemblers.jl at master · PetrKryslUCSD/Elfel.jl · GitHub

1 Like
6

I’ll confess that I find the sparse(i,j,v,m,n) constructor rather non-ergonomic. It’s adopted verbatim from MATLAB (where almost anything other than Array{Float64,N} is a pain to use) and I consider it quite clumsy in Julia.

I would have found a constructor based on an iterable (e.g., Vector) of Pair{CartesianIndex{N},T} to be much easier to use, for example. Then you could add entries by push!ing onto a single vector instead of 3.

Sometimes I grit my teeth and build the i/j/v vectors separately. Other times I end up building my single vector (with eltype Tuple{Int,Int,T} or the above Pair) and then decompose it into the 3 to pass to sparse. As suggested above, this isn’t a very hard transformation to do. It’s also something easily done by a package without the need to modify SparseArrays, so perhaps we accept this ship as “sailed”.

4 Likes