struct SparseMatrixCSC cannot be constructed correctly with "repeating" entries

[FangXieTHU]

In many use cases, constructing CSC matrices with "repeating" entries in each column is usually convenient and memory efficient.
In particular, when the same entry appears multiple times, we usually expect that the values for these entries will be added together.
However, it seems that directly using SparseMatrixCSC does not work as intended.

In the following minimal example, I tried to construct an identity matrix using both the IJV format (sparse matrix A) and using CSC format (sparse matrix B).

The entry at [1, 1] appeared twice, each time with value 0.5. Therefore, I would expect the matrix should behave as an 3 dimensional identity matrix.

Using IJV format, the matrix behaves exactly as what I expected. However, the direct construction B using CSC format does not always work properly in the following ways:

1. It cannot be displayed properly. The matrix element at [1, 1] shows 0.5.
2. Converting it into dense matrix by Matrix(B) also returns a matrix with 0.5 at [1, 1].
3. Interestingly, the matrix-vector multiplication operators involving sparse matrix B works correctly.

I am wondering if this behavior of CSC matrix construction is intended, or if this is a glitch in the implementation of display and Matrix functions?

If it is not intended, I believe it would be helpful to get it patched.
If it is intentional, I would like to ask whether it might be possible to introduce an alternative method for constructing CSC-format sparse matrices with repeating entries.

Many thanks to the SparseArrays.jl team for their help and support.

using SparseArrays

col_idx::Vector{Int64} = [1, 1, 2, 3]
row_idx::Vector{Int64} = [1, 1, 2, 3]

colptr::Vector{Int64} = [1, 3, 4, 5]
val::Vector{Float64} = [0.5, 0.5, 1.0, 1.0]

A = sparse(col_idx, row_idx, val, 3, 3)
B = SparseMatrixCSC(3, 3, colptr, row_idx, val)

display(A)
display(B)
display(Matrix(B))

x::Vector{Float64} = [1.0, 1.0, 1.0]

display(A * x)
display(B * x)
display(Matrix(B) * x)

Output:

3×3 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
 1.0   ⋅    ⋅
  ⋅   1.0   ⋅
  ⋅    ⋅   1.0
3×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 0.5   ⋅    ⋅
  ⋅   1.0   ⋅
  ⋅    ⋅   1.0
3×3 Matrix{Float64}:
 0.5  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
3-element Vector{Float64}:
 1.0
 1.0
 1.0
3-element Vector{Float64}:
 1.0
 1.0
 1.0
3-element Vector{Float64}:
 0.5
 1.0
 1.0

[dkarrasch]

The SparseMatrixCSC docstring states that

The standard way of constructing SparseMatrixCSC is through the sparse function.

The 5-arg SparseMatrixCSC constructor appears to be an "expert" function, that doesn't do repeated entry checking, and accepts the arguments as is. Matrix vector multiplication then seems to treat the double entries just as what they are, whereas the Matrix constructor jumps to the next entry once an entry has been filled with a nonzero entry.

So, out of curiosity, why not use the recommended constructor, which also merges repeated entries into a single one?

[FangXieTHU]

So, out of curiosity, why not use the recommended constructor, which also merges repeated entries into a single one?

From what I understand about the documents, the standard sparse function uses IJV format. When the matrix gets very large, for example a few millions dimensions with thousands of nonzero entries in each column, the vector row_idx::Vector{Int64} can be thousands of times longer than constructing colptr::Vector{Int64}. Directly using the CSC format is a lot more memory efficient than the IJV format.

Besides, different columns can be constructed separately using CSC format, which means matrix * vector operations can be parallelized, and will be very helpful in Lanczos algorithms. This is quite commonly used in numerical quantum many-body physics.

[dkarrasch]

Oh, I missed the slight difference. So, what needs to be fixed is the Matrix constructor, right? Detecting and merging repeated entries is more of an improvement. I do wonder, though, if we want to keep checking mechanisms away from the actual constructor and put it in some sparse method, with a different signature.

[j-fu]

You might try out ExtendableSparse.jl (I am the main author). This packages internally uses a linked list intermediate format for a purpose similar to yours.

EDIT: Ah, and: you should be aware of DifferentiationInterface.jl and SparseConnectivityTracer.jl, which show great performance for setting up sparsity pattenrns from nonlinear operators, see JuliaDiff/DifferentiationInterface.jl#774 (comment)

[gdalle]

Perhaps https://sparsearrays.juliasparse.org/dev/#SparseArrays.sparse! is more relevant here?

[ViralBShah]

The SparseMatrixCSC constructor by definition is lightweight and does not deal with such things. As folks have pointed out, sparse is what you use for all those features.

Perhaps we can do a better job of documenting that.

[FangXieTHU]

As folks have pointed out, sparse is what you use for all those features.

It seems that the input that the function sparse needs is IJV format instead of CSC or CSR. Unfortunately, IJV format is a lot less memory efficient. I am not yet sure whether sparse! can be useful in this case.

[dkarrasch]

I think the proposal is to have a new method of sparse that takes CSC-type arguments and does what sparse does with IJV-type arguments. The point is to not put the load on the SparseMatrixCSC constructor, but keep it leightweight.

[j-fu]

I think a SparseMatrixCSC with several occurrences of the same row index in a column would'nt make it much easier to add new entries. Adding new entries still would trigger shifting portions of rowval and nzval.

However, there is a linked list sparse matrix format described in a whitepaper by Y. Saad. See also exercise P.3-16 in his book and the sparskit source. This format implements a linked list in the rowval field, and therefore new entries are added at the end of nzval, and existing entries are easily updated. Unlike the IJV format, this avoids to have repeated entries for the same index.

I have a variant of an implementation in ExtendableSparse.jl: https://github.com/WIAS-PDELib/ExtendableSparse.jl/blob/master/src/matrix/sparsematrixlnk.jl .
(Indeed there are several variants in an attempt to support multithreaded assembly)

Not sure if it makes sense to add this here as an additional format. If the wish occurs, I would try to help.