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:

JuliaDiff/FiniteDiff.jl#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 #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! :)