Julia: AbstractSparseArray Interface

[ChrisRackauckas]

It would be nice to have a clearly defined AbstractSparseArray <: AbstractArray with an interface that defines where the non-zero values are, and that way generic looping over sparse arrays can be defined. Special matrix types like Tridiagonal should fall under AbstractSparseArray to allow for generic sparse iteration, and then it would make it easier for things like @dlfivefifty 's BandedMatrices.jl to plug into a more generic system of sparsity support.

[ChrisRackauckas]

@timholy 's ArrayIteration.jl has some ideas here:

https://github.com/timholy/ArrayIteration.jl

But to really make that work well I think we need a standard interface in Base for everyone's array types to conform to in order to allow generic codes to work on it. The problem of course is to make generic sparse iteration efficient, if possible.

[dlfivefifty]

In theory this is a good idea. In practice, I'm doubtful structured sparse matrices (BandedMatrix, SymTridiagonal, etc.) will benefit that much, as they probably need specialized linear algebra implementations to achieve top speeds. I don't see a general purpose interface that doesn't know the structure being fast enough.

[KristofferC]

Why do you want to loop over sparse array instead of using a higher level function like map, mapreduce etc?

[dlfivefifty]

map and broadcast are almost useless without a fill-in value. (exp.(A) is filled in with 1s)

[ChrisRackauckas]

I want to easily get information about sparsity patterns to be able to do sparse differentiation described here:

https://github.com/JuliaDiffEq/DiffEqDiffTools.jl/issues/29

I assume that this kind of stuff will be needed on a lot of different types, not only for our finite difference uses, but also in eventual AD implementations as well. I assume there must be some better way to handle this than to specialize to every matrix type, though I could be wrong there.

[ChrisRackauckas]

Note that I am not saying anything about creating a generic sparse matrix multiplication or anything like that, that would be silly. Of course that requires utilizing the representation in order to be efficient at all. But I am not sure that the query for the sparsity structure has to be different for each sparse type, and as @dlfivefifty noted sometimes a simple loop over non-zero elements is required. If there was an easy way to get an efficient iterator for the right way to walk through the non-zero elements (along with gathering the indices) then that would be helpful.

[dlfivefifty]

I think at the basic level we would want something like the following to work for general sparse matrices:

convert(::Type{SparseMatrixCSC}, A::AbstractSparseMatrix) = sparse(nonzerorows(A), nonzerocolumns(A), nonzerovalues(A))

But it raises some questions on whether these should be separate iterators, or a single iterator whose items are (i,j,A[i,j]) for all i,j such that A[i,j] ≠ 0.

[mbauman]

We've done quite a bit of brainstorming about this in https://github.com/JuliaLang/julia/issues/15648 — which is precisely what led to Tim's ArrayIteration.jl package.

[KristofferC]

Indeed, this falls very well in scope of #15648, so closing.

[KristofferC]

FWIW, I posted an example iterator that gives (i, j, A[i,j]) for sparse matrices at https://discourse.julialang.org/t/how-to-get-the-row-and-column-value-for-sparsematrixcsc-elements/8613/16.

[mbauman]

I do think there could be half-way solutions here that could be quicker and easier to implement that might satisfy the root issue here. They won't be the whole ArrayIteration thing, but it seems clearer in my mind now that SparseArrays is a package. Heck, SparseArrays could just export a new findstored (and maybe even an eachstored iterator) function that does the "right thing" for dense arrays.

I imagine that doesn't solve your real problem, though, since you want "not-all-indices-are-stored" for dispatch purposes. That gets into the array storage trait that Sheehan's been working on — I've been imagining that as a 1.1 feature.

[Sacha0]

Indeed, this falls very well in scope of #15648, so closing.

Though indeed the part of this issue concerning iteration falls under #15648, the broader issue of clarification of the AbstractSparse(Array|Vector|Matrix) interfaces seems to stand alone. I have long thought of opening a similar issue, as coding against the former abstractions isn't possible at present due to lack of well-defined interfaces. So I might advocate reopening this issue with focus on its distinct part. Best! :)

[mbauman]

Yes, I agree. Key question here is _what does it mean to be an AbstractSparseArray?_ (akin to #10064). This does dovetail with the implementation discussions in #15648, but let's try to focus on the key properties that make an array an AbstractSparseArray. Is it that the non-stored values are zeros? Does it have something to do with the complexity of indexing? Or the typical/ideal proportion of stored to non stored? What should this say?

help?> AbstractSparseArray
search: AbstractSparseArray AbstractSparseMatrix AbstractSparseVector

  No documentation found.

[StefanKarpinski]

AbstractSparseArray to me means that the storage size of the array is << O(m*n) – that's what sparse data structures are all about.

[ViralBShah]

Is this the right place for this discussion? Developing a new clean interface is something that should happen in JuliaSparse ideally.

[ViralBShah]

Saying the obvious here - J. H. Wilkinson's informal working definition of a sparse matrix was "any matrix with enough zeros that it pays to take advantage of them". This basically translates to having algorithms that can operate in time proportional to nnz(S), for a sparse matrix S.

I am closing this - as there isn't much actionable here. If people want to keep it around for discussion, we can reopen.

[mbauman]

Here's the action — and why I re-opened this last time:

help?> AbstractSparseArray
search: AbstractSparseArray AbstractSparseMatrix AbstractSparseVector

  No documentation found.

I'd love for this to go one step further and define a handful of core helper functions that we then can use with all arrays that are particularly efficient for sparse implementations. Of course, that's more than just a doc change, and somewhat more targeted by https://github.com/JuliaLang/julia/issues/26613 (which IMO is also a partial duplicate).

[ViralBShah]

Maybe we can continue that in #26613. I suspect that it isn't easy to come up with a data structure agnostic API that is fast.