Cancellation with sparse arrays

[gdalle]

Consider a function $f(t)$ such that $f(0)$ is the fully zero sparse vector. Apparently, ForwardDiff cannot compute $f'(0)$.

julia> using ForwardDiff: derivative

julia> using SparseArrays: sparse

julia> x = sparse([1.0]);

julia> dx = [1.0];

julia> derivative(t -> x + t * dx - x, 0)  # incorrect
1-element SparseVector{Float64, Int64} with 0 stored entries

julia> derivative(t -> Vector(x) + t * dx - Vector(x), 0)  # correct
1-element Vector{Float64}:
 1.0

[mcabbott]

The answer here is changed by #481:

julia> derivative(t -> x + t * dx - x, 0)  # was zero above
1-element SparseArrays.SparseVector{Float64, Int64} with 1 stored entry:
  [1]  =  1.0

(jl_PuuLcA) pkg> st ForwardDiff
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_PuuLcA/Project.toml`
  [f6369f11] ForwardDiff v0.11.0-DEV `https://github.com/JuliaDiff/ForwardDiff.jl.git#master`

However, it's not entirely clear to me what the right behaviour is. It depends whether you think the sparse array is just a compact representation of a dense array or whether you think the sparsity information is structural, e.g. a sparse matrix represents a graph.

ChainRules takes the latter view, that the location of zeros is a structural feature like Diagonal. In this case, the bug in this particular example is that x's structure is that it has one nonzero entry, but this is lost by the function t -> x + t * dx - x, which promotes an accidental zero to a structural feature.

julia> let t = 0
         x + t * dx - x
       end
1-element SparseArrays.SparseVector{Float64, Int64} with 0 stored entries

julia> let t = eps()
         x + t * dx - x
       end
1-element SparseArrays.SparseVector{Float64, Int64} with 1 stored entry:
  [1]  =  2.22045e-16

Edit: The reason 481 works here is that it alters the check SparseArrays is using to decide whether to drop the zero. Even if we think the let t = 0 answer is a bug, the Dual case seems correct:

julia> let t = ForwardDiff.Dual(0, 1)
         x + t * dx - x
       end
1-element SparseArrays.SparseVector{ForwardDiff.Dual{Nothing, Float64, 1}, Int64} with 1 stored entry:
  [1]  =  Dual{Nothing}(0.0,1.0)

[gdalle]

Your explanation makes sense, and I also think that the zeros should be structural (so I actuall end up disagreeing with the upcoming release 🤣)

Perhaps my real gripe is that x - x doesn't have the same sparsity structure as x. But I imagine there is a good reason?

[mcabbott]

I am not sure there's a good reason. Somewhere there was an issue complaining about this kind of behaviour, and my memory is that most were horrified that values could alter sparsity in something like this:

julia> x = sparse([1.0]);

julia> x .*= 0;

julia> x
1-element SparseArrays.SparseVector{Float64, Int64} with 0 stored entries

(May have had some mistakes above at first, sorry.)

[gdalle]

I'll try to dig up the issue, wanna add to the outrage

[gdalle]

For reference: JuliaSparse/SparseArrays.jl#190