Dot operation broadcast for a "StaticSparseVector"

For various reasons, I want to create a static sparse vector (SSV) type. In contrast to what is done (for matrices) in LuxurySparse.jl, I want the sparsity structure to be part of the data type.

I do that by using a recursive datatype SSV parametrized by an index I that has a value for this index, and a next containing the rest of the vector.

struct EmptySSV{L} end
struct SSV{I,V,N} # Static Sparse Vector, recursive, in which the sparsity pattern is part of the type
    next ::N
SSV{I}(value,next) where{I} = SSV{I,typeof(value),typeof(next)}(value,next)

function StaticVector(l,i,v) # general purpose, type unstable constructor
    out     = EmptySSV{l}()
    for j   = reverse(eachindex((i)))
        out = SSV{i[j]}(v[j],out)
    return out

for OP∈(:(+),:(-))
    @eval function Base.$OP(v1::SSV{I1,V1,N1},v2::SSV{I2,V2,N2}) where{I1,V1,N1,I2,V2,N2}
        if     I1==I2 SSV{I1}($OP(v1.value,v2.value),$OP(,
        elseif I1< I2 SSV{I1}(    v1.value          ,$OP(,v2     ))
        elseif I1> I2 SSV{I2}(             v2.value ,$OP(v1     ,
    @eval Base.$OP(v::SSV{I,V,N}  ,e::EmptySSV{L }) where{I,V,N,L} = SSV{I}(v.value,$OP(,e     ))
    @eval Base.$OP(e::EmptySSV{L },v::SSV{I,V,N  }) where{I,V,N,L} = SSV{I}(v.value,$OP(e     ,
    @eval Base.$OP( ::EmptySSV{L }, ::EmptySSV{L }) where{L      } = EmptySSV{L}()
    @eval Base.$OP( ::EmptySSV{L1}, ::EmptySSV{L2}) where{L1,L2  } = error("Attempt to add vectors of ≠ lengths")
# Tests

v = StaticVector(10,[3,5],[3.,5.])
w = StaticVector(10,[2,5],[2.,5.])

@show v[3]
@show v[5]
#@show v[1]   # throws "Not stored"
#@show v[11]  # throws "Out of bounds"
@show v+w===StaticVector(10,[2,3,5],[2.,3.,10.])

Everything goes fine until I try to make it possible to have vectorised operations on these vectors. I tried to overload broadcast:

function Base.Broadcast.broadcast(op,v1::SSV{I1,V1,N1},v2::SSV{I2,V2,N2}) where{I1,V1,N1,I2,V2,N2}
    if     I1==I2 SSV{I1}(broadcast(op,v1.value,v2.value),broadcast(op,,
    elseif I1< I2 SSV{I1}(             v1.value          ,broadcast(op,,v2     ))
    elseif I1> I2 SSV{I2}(                      v2.value ,broadcast(op,v1     ,
Base.Broadcast.broadcast(op,v::SSV{I,V,N}  ,e::EmptySSV{L }) where{I,V,N,L} = SSV{I}(v.value,broadcast(op,,e     ))
Base.Broadcast.broadcast(op,e::EmptySSV{L },v::SSV{I,V,N  }) where{I,V,N,L} = SSV{I}(v.value,broadcast(op,e     ,
Base.Broadcast.broadcast(op, ::EmptySSV{L }, ::EmptySSV{L }) where{L      } = EmptySSV{L}()
Base.Broadcast.broadcast(op, ::EmptySSV{L1}, ::EmptySSV{L2}) where{L1,L2  } = error("Attempt to add vectors of ≠ lengths")

but my methods do not seem to be called, I get an error that no method matching iterate(::SSV{3, Float64, SSV{5, Float64, EmptySSV{10}}}) which indeed I did not define, but did not call from my broadcast. So I was not supposed to overload broadcast maybe?

How are broadcast and iterate connected? I do not believe Julia would be able, with me implementing ìterate alone, to do the sort of recursive braiding I do in my + and - methods, and so I am concerned that I do not have all the bits in place…

Looking for inspiration, I find that SparseArray.jl has functions like _bcast_binary_map. … but called by what?

Anyone knows how things look under the hood?

Maybe I can rephrase my question better: My SSV type has the peculiar property that getindex is not implemented for it. This is because this vector (by design) can combine elements of different types, and hence deciding what zero to return on indexing an element not explicitely stored is not possible. Nonsense? Well, this workhorse type is to be wrapped into more specific types…

What do I need to implement to allow broadcasting on my own vector type?
Is it a requirement for broadcasting that indexing must be possible?

There’s a list of functions required for broadcasting: Interfaces · The Julia Language

But if you can’t know the type of the zeros, how will broadcasting work? What happens with isa.(Integer, x)?

Hi Gustaphe,

I will study the link you gave me in fine details! Give me a few days…

Very good question indeed! Your question guides me: this should return a StaticSparseVector of booleans. A central issue: I do not want broadcasting to result in full vectors (unless… see below)

What if I add a Vector to a StaticSparseVector. I would return the same Vector, except with an additions where relevant, and so I never have to find out “what kind of zero”. Indeed, for efficiency, I do not want to spend time adding zeros.

The problem with that approach is, a sparse vector if booleans is just a representation of a vector where every element is false except the specified. But that doesn’t really make sense for you. How about

f(x) = x isa Integer || iszero(x)

Should this be true or false for a sparse vector with only Integers in it? Because intuitively I’d say true, but with your suggestion it’s either false or undefined.

Touché, excellent!

OK, my StaticSparseVector is a weird customer, and it is hereby renamed HeterogenousStaticSparseVector (HSSV). You point out that sparse and heterogenous put together make it impossible to implement the whole sparse vector interface, including, quite possibly, broadcasting. I’ll chew on this!

And indeed, this is not (or no longer!!! :grin:) my intention to have the whole interface (though I need a.*b, and this operation has a well defined result…): HSSV will be wrapped into various other containers with a well defined eltype, to which we will be able to give a reasonable “AbstractVector” interface.

So why heterogenous? Well, you can do some cool stuff with it. One favorite example is implementing a matrix as sparse vector of sparse vector. How could I go static (stay off the heap) otherwise? A matrix is a row of columns, and each column has a different type (since sparsity pattern is part of the vector type). Of course for a matrix, you can (should!) implement this is other ways, but this flexibility will be valuable for forward automatic differentiation to higher orders with sparse partials, which are themselves forwarddiff objects with sparse partials…

Thank you for this discussion, this helps me a lot!