operations on subtypes of AbstractSparseVector

[benlorenz]

We are using a custom sparse vector type in Polymake.jl (backed by a C++ library) that is a subtype of AbstractSparseVector. With the recent merge of #175 and the stdlib bump in julia nightly our tests started to fail because scalar*vector multiplication now converts the vector to a julia SparseArrays.SparseVector instead of using the generic AbstractArray operation as a fallback:

Current nightly:

julia> versioninfo()
Julia Version 1.9.0-DEV.1101
Commit 5e081d6101a (2022-08-09 04:39 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.5 (ORCJIT, tigerlake)
  Threads: 1 on 8 virtual cores

julia> using Polymake
[...]

julia> import SparseArrays

julia> v = Polymake.SparseVector{Int64}([1,0,0,0,0,3])
pm::SparseVector<long>
(6) (0 1) (5 3)

julia> v isa SparseArrays.AbstractSparseVector
true

julia> 2*v
6-element SparseVector{Int64, Int64} with 2 stored entries:
  [1]  =  2
  [6]  =  6

julia> @which 2*v
*(a::Number, x::Union{AbstractSparseVector{Tv, Ti}, SubArray{Tv, 1, <:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false}, SubArray{Tv, 1, <:AbstractSparseVector{Tv, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}}, false}} where {Tv, Ti})
     @ SparseArrays ~/software/polymake/julia/julia/julia-5e081d6101/share/julia/stdlib/v1.9/SparseArrays/src/sparsevector.jl:1541

With an older julia nightly we had:

julia> versioninfo()
Julia Version 1.9.0-DEV.926
Commit 51bb96857d (2022-07-08 00:05 UTC)
[...]

julia> v = Polymake.SparseVector{Int64}([1,0,0,0,0,3])
pm::SparseVector<long>
(6) (0 1) (5 3)

julia> 2*v
pm::SparseVector<long>
(6) (0 2) (5 6)

julia> @which 2*v
*(A::Number, B::AbstractArray)
     @ Base arraymath.jl:21

The relevant operations are defined here:

SparseArrays.jl/src/sparsevector.jl

Lines 1540 to 1542 in c4b46a9

	(*)(x::SparseVectorUnion, a::Number) = SparseVector(length(x), copy(nonzeroinds(x)), nonzeros(x) * a)
	(*)(a::Number, x::SparseVectorUnion) = SparseVector(length(x), copy(nonzeroinds(x)), a * nonzeros(x))
	(/)(x::SparseVectorUnion, a::Number) = SparseVector(length(x), copy(nonzeroinds(x)), nonzeros(x) / a)

and the SparseVectorUnion was extended to the abstract type here.

For the negate operation this does not happen and

julia> -v
pm::SparseVector<long>
(6) (0 -1) (5 -3)

still works as expected in the current nightly.

I am not really sure if this is a bug or by design and we need to implement all those operations explicitly?

[SobhanMP]

maybe the function could convert back to the type of x after the operation? or limit it to AbstractSparseArrayC (after #201). the problem with the abstractarray fallback is that it does not understand sparsity, hence does too many operations.

on a side note the speed is really bad at the moment. 50~100 times worst.

julia> a = sprandn(1000, 0.05);

julia> v = Polymake.SparseVector{Float64}(Vector(a));

julia> @btime 2 * $a;@btime 2 .* $a;@btime 2 * $v;@btime 2 .* $v;
  62.529 ns (2 allocations: 896 bytes)
  112.026 ns (4 allocations: 928 bytes)
  11.710 μs (1 allocation: 16 bytes)
  11.581 μs (1 allocation: 16 bytes)

[benlorenz]

maybe the function could convert back to the type of x after the operation?

That will probably be difficult as x might be just a SparseVectorView? Also I don't think any AbstractSparseVector needs to have a constructor with (length, indices, values) as this is rather specific to the implementation.

or limit it to AbstractSparseArrayC (after #210).

I think limiting these operations would be reasonable yes.

the problem with the abstractarray fallback is that it does not understand sparsity, hence does too many operations.
on a side note the speed is really bad at the moment. 50~100 times worst.

This is expected and not really important, the type is almost exclusively used for passing data to back and forth and not for computations in julia.
Within the C++ library there is a proper scalar multiplication operation that is just not exposed to julia yet (and there are various kinds of lazy operations which would delay any multiplications until the elements of the new vector are used).

[SobhanMP]

seems like i've linked the wrong issue, meant to say #201

[SobhanMP]

Also I don't think any AbstractSparseVector needs to have a constructor with (length, indices, values) as this is rather specific to the implementation.

i think a move constructor (one that does not reallocate) from SparseVector can be a necessity? this way packages can define it as they see fit

[dkarrasch]

Perhaps we should redo parts of the "offending" #175? Including the abstract type in the union has side effects beyond what was to be fixed in that PR, but I don't think we need to include it.

[dkarrasch]

I believe the issue is already resolved. We currently have

const SparseVectorUnion{Tv,Ti} = Union{AbstractCompressedVector{Tv,Ti}, SparseColumnView{Tv,Ti}, SparseVectorView{Tv,Ti}}

where

abstract type AbstractCompressedVector{Tv,Ti} <: AbstractSparseVector{Tv,Ti} end

So, if your custom sparse vector type subtypes AbstractSparseVector and not AbstractCompressedVector, then it is no longer included in SparseVectorUnion and all should be well. It may be that the recent change has not yet propagated to Julia Base nightly.

[benlorenz]

Thanks for the update! We only subtype AbstractSparseVector so this should be fine.

It is indeed not yet in julia nightly, but I will try it again after the next stdlib update and re-open the issue if necessary.