Inconsistent behavior in assignment to sparse matrices of Dual numbers

[moyner]

trafficstars

Let us say that we have a sparse matrix where the entries are Dual numbers.

using ForwardDiff, SparseArrays
d = ForwardDiff.Dual(1.0, 2.0)
a = sparse([0 d])

Let us define two zero values with different partials. One with zero partials (true zero) and one with non-zero partials (false zero):

# 0, with ∂ = 1
d_false_zero = ForwardDiff.Dual(0.0, 1.0)
# 0, with ∂ = 0
d_true_zero = ForwardDiff.Dual(0.0, 0.0)

If we replace the existing entry at [2], everything is well:

## Insertion into existing element keeps partials
a[2] = d_false_zero
println("Existing element:")
println(a[2].partials == d_false_zero.partials)

This gives the expected output:

Existing element:
true

Now let us try to insert our false zero into the first position, which we know isn't yet stored:

## Insertion into unstored element drops partials
a[1] = d_false_zero
println("Zero element:")
println(a[1].partials == d_false_zero.partials)

Zero element:
false

The partials were lost, due to the iszero check here in SparseArrays: https://github.com/JuliaLang/julia/blob/cb30aa7a089f4d5f1c9f834d96c3788f2e93ebcc/stdlib/SparseArrays/src/sparsematrix.jl#L2652

Some quick and ugly type piracy can fix this, but may have other unintended consequences:

import Base.iszero
Base.iszero(d::ForwardDiff.Dual) = iszero(d.value) && iszero(d.partials)

I'm not sure what the intended behavior should be, since the entry is technically zero, but not a neutral element for addition of Duals.

[mcabbott]

Fixed by #481, I think:

julia> println(a[1].partials == d_false_zero.partials)
true

the entry is technically zero

That's the question. Opinions vary but the most mathematical literature seems to regard Dual as inheriting equality from Tuple, and hence only being zero when the partials are zero. Computer implementations appear to have been written more or less by taking a coin toss.

The alternative solution to this would be for sparse array setindex! always to create an entry, even if you write a zero. If you think of them as efficient storage for dense matrices this is a bad idea. If you think of them as graphs then perhaps the absence of an edge isn't the same as having an edge whose weight is currently zero. But Julia's SparseArrays seems to take the former view.

[SobhanMP]

I think this still an issue on 1.8.5

julia> using SparseArrays, ForwardDiff


julia> ForwardDiff.gradient([1.0, 1.0]) do x
           a = spzeros(eltype(x), 2, 2)
           a[1, 1] = x[1] - x[2]
           a[1, 1]
       end # [0.0, 0.0]
2-element Vector{Float64}:
 0.0
 0.0

julia> ForwardDiff.gradient([1.0, 1.0]) do x
           a = zeros(eltype(x), 2, 2)
           a[1, 1] = x[1] - x[2]    
           a[1, 1]
       end # [1.0, -1.0]
2-element Vector{Float64}:
  1.0
 -1.0

The alternative solution to this would be for sparse array setindex! always to create an entry, even if you write a zero.

@LilithHafner was arguing for something like that... didn't think i'll run into it.

I guess https://github.com/JuliaSparse/SparseArrays.jl/pull/296 fixes it that way.

[mcabbott]

These agree on ForwardDiff#master, i.e. v0.11-dev.

After we merged https://github.com/JuliaDiff/ForwardDiff.jl/pull/481 some people decided they liked the bugs, and so it was pulled from 0.10.x ... and we're stalled on deciding whether to tag 0.11 or 1.0.

[SobhanMP]

but the new is zero check is v !== zero(eltype(A)). wouldn't it fix the problem?

d = ForwardDiff.Dual(1.0, 2.0)
d - 1 !== zero(d) # true

[LilithHafner]

This issue is fixed on Julia master. It is still present on 1.9.0-rc2 because https://github.com/JuliaSparse/SparseArrays.jl/pull/296 hasn't landed there but should be fixed on 1.9.0.

julia> using SparseArrays, ForwardDiff

julia> ForwardDiff.gradient([1.0, 1.0]) do x
           a = spzeros(eltype(x), 2, 2)
           a[1, 1] = x[1] - x[2]
           a[1, 1]
       end # [0.0, 0.0]
2-element Vector{Float64}:
  1.0
 -1.0

julia> ForwardDiff.gradient([1.0, 1.0]) do x
           a = zeros(eltype(x), 2, 2)
           a[1, 1] = x[1] - x[2]    
           a[1, 1]
       end # [1.0, -1.0]
2-element Vector{Float64}:
  1.0
 -1.0

julia> VERSION
v"1.10.0-DEV.1150"

(jl_mKouP7) pkg> st
Status `/private/var/folders/hc/fn82kz1j5vl8w7lwd4l079y80000gn/T/jl_mKouP7/Project.toml`
  [f6369f11] ForwardDiff v0.10.35

[LilithHafner]

(fixed on Julia 1.9.0-rc3)