SparseArrays.SparseMatrixCSC support

[odow]

Adding #4003 made me wonder if we should do something like:

julia> using JuMP, SparseArrays

julia> _check_rectangular(::Any) = error("Unsupported iterator type")
_check_rectangular (generic function with 1 method)

julia> function _check_rectangular(indices::Containers.NestedIterator)
           if length(indices.iterators) != 2
               error(
                   "Unable to construct a SparseMatrixCSC: the set must be " *
                   "two-dimensional. Got: $(length(indices.iterators))",
               )
           end
           I = indices.iterators[1]()
           if I != Base.OneTo(length(I))
               error(
                   "Unable to construct a SparseMatrixCSC: the first dimension must " *
                   "be of the form `1:M`. Got: $I",
               )
           end
           J = indices.iterators[2](1)
           if J != Base.OneTo(length(J))
               error(
                   "Unable to construct a SparseMatrixCSC: the second dimension " *
                   "must be of the form `1:N`. Got: $J",
               )
           end
           for i in I  # Check that there is no triangular indexing
               if indices.iterators[2](i) != J
                   error(
                       "Unable to construct a SparseMatrixCSC: the second dimension " *
                       "must not depend on the value of the first dimension.",
                   )
               end
           end
           return length(I), length(J)
       end
_check_rectangular (generic function with 2 methods)

julia> function Containers.container(
           f::Function,
           indices,
           ::Type{SparseArrays.SparseMatrixCSC},
       )
           m, n = _check_rectangular(indices)
           x = map(args -> args => f(args...), indices)
           I, J, V = first.(first.(x)), last.(first.(x)), last.(x)
           return SparseArrays.sparse(I, J, V, m, n)
       end

julia> model = Model()
A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> @variable(model, x[i in 1:4, j in 1:4; isodd(i + j)], container = SparseMatrixCSC)
4×4 SparseMatrixCSC{VariableRef, Int64} with 8 stored entries:
 ⋅       x[1,2]  ⋅       x[1,4]
 x[2,1]  ⋅       x[2,3]  ⋅
 ⋅       x[3,2]  ⋅       x[3,4]
 x[4,1]  ⋅       x[4,3]  ⋅

The downside is that it is pretty finicky with what shapes we can detect, and iterating over every element isn't very efficient.

It might be better to leave the suggestion in the docs?

[blegat]

The reason we don't use a sparse array and we use our own type, even if the axes are OneTo is that they have different semantics for when an entry is missing. In our sparse array, accessing a missing entry is an error while for sparse arrays it gives zero.

[odow]

I don't mean it should be the default. I mean should we make it opt-in for people that want it.

[blegat]

Yes, I think it would be a good idea, we just need to warn of the potentially surprising semantic difference in the docs.

[odow]

Find a concrete example of the method errors

[odow]

I remember now.

One thing is that broadcasting does not skip nulls (because they default to 0, and your function might return a non-zero) (x-ref #1876):

julia> using JuMP, SparseArrays

julia> model = Model()
A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> @variable(model, x[i in 1:5] >= i)
5-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
 x[4]
 x[5]

julia> A = SparseArrays.sparse([1, 1, 2, 2, 3], [1, 2, 2, 3, 3], x)
3×3 SparseMatrixCSC{VariableRef, Int64} with 5 stored entries:
 x[1]  x[2]  ⋅
 ⋅     x[3]  x[4]
 ⋅     ⋅     x[5]

julia> lower_bound.(A)
ERROR: Cannot call lower_bound with 0 because it is not an affine expression of one variable.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] _eval_as_variable(::typeof(lower_bound), ::AffExpr)
   @ JuMP ~/.julia/packages/JuMP/RGIK3/src/aff_expr.jl:975
 [3] lower_bound
   @ ~/.julia/packages/JuMP/RGIK3/src/aff_expr.jl:1003 [inlined]
 [4] _noshapecheck_map(::typeof(lower_bound), ::SparseMatrixCSC{VariableRef, Int64})
   @ SparseArrays.HigherOrderFns ~/.julia/juliaup/julia-1.10.9+0.x64.apple.darwin14/share/julia/stdlib/v1.10/SparseArrays/src/higherorderfns.jl:180
 [5] copy
   @ ~/.julia/juliaup/julia-1.10.9+0.x64.apple.darwin14/share/julia/stdlib/v1.10/SparseArrays/src/higherorderfns.jl:195 [inlined]
 [6] materialize(bc::Base.Broadcast.Broadcasted{SparseArrays.HigherOrderFns.SparseMatStyle, Nothing, typeof(lower_bound), Tuple{SparseMatrixCSC{VariableRef, Int64}}})
   @ Base.Broadcast ./broadcast.jl:903
 [7] top-level scope
   @ REPL[217]:1

The error message is okay, but somewhat confusing.

Other broadcasts work:

julia> A.^2
3×3 SparseMatrixCSC{QuadExpr, Int64} with 5 stored entries:
 x[1]²  x[2]²  ⋅
 ⋅      x[3]²  x[4]²
 ⋅      ⋅      x[5]²

I think we've fixed most of the bugs I was thinking about, e.g., #2309, #2162, #1950

[odow]

This discourse thread is a fairly good summary of all the reasons why we shouldn't add SparseMatrixCSC support: https://discourse.julialang.org/t/jump-and-sparsearrays/129264

[WalterMadelim]

I think SparseMatrixCSC is mainly developed as a front API of SuiteSparse, where optimized algorithms (including hardware acceleration) exist. Therefore, if users cannot take advantage of this data format in their particular non-pure-numeric applications, then why do we need them?

[odow]

I think I'm going to close this issue as won't fix. User's can manually construct sparse arrays of JuMP types if they want to. But there are enough small issues that if we pretend we have first-class support we'll just get a stream of discourse questions for features that we can't/won't fix.