Inequality, ε-balls, and accidental structural zeros

[mcabbott]

Consider this function:

sq(x) = x==1 ? one(x) : x^2

@test FiniteDifferences.central_fdm(5, 1)(sq, 1) ≈ 2.0
@test_broken ForwardDiff.derivative(sq, 1.0) == 2.0

Here ForwardDiff gets the wrong answer, according to your first calculus class: The derivative is defined by taking limits, evaluating sq(x + ε) for some small ε, and these always see the continuum x^2 not the special point.

One to think about this is to say that x==1 really means abs(x-1) < ζ with some tiny ζ, which we keep finite until we are sure we aren't confused. The calculus class assumption is that ζ << ε.

The assumption of ForwardDiff is the opposite. Its Dual(x,1) encodes a perturbation x + 1ε with ε smaller than everything else around, and in particular ε << ζ. Or in other words, sq is viewed as being piecewise continuous, with a small flat area of width 2ζ, which is still large enough for us to see that its slope is zero.

Of course nobody really writes contrived examples like sq. But they do write things like this:

function prod1(xs::Vector)
    p = one(eltype(xs))
    for x in xs
        p = p * x
    end
    p
end

function prod2(xs::Vector)
    p = one(eltype(xs))
    for x in xs
        p = p * x
        p == 0 && break  # exit early once you know the answer
    end
    p
end

@test ForwardDiff.gradient(prod1, [1,2,0,4,0,6]) == zeros(6)
@test_broken ForwardDiff.gradient(prod2, [1,2,0,4,0,6]) == zeros(6)

This has almost the same problem as #197, where det(A) tests for istriu(A) || istril(A) before calling a simpler branch. The fact that f(x,y) == g(x,y) when y==0 does not imply that df/dy == dg/dy. So it seems AD ought not to take that branch.

In which case we want something like this:

Base.:(==)(x::Dual, y::Int) = x.value == y && iszero(x.partials)
Base.:(!=)(x::Dual, y::Int) = x.value != y || !iszero(x.partials)

This fixes the tests above, and (a slightly more careful version) fixes #197 and #407.

However, it means that fun(Dual(x,1)).value need not be equal to fun(x), on a discontinuous function. Although fun(Dual(x,0)).value should still be equal, @assert zero(x) == 0 isn't broken, and there should be no problems where functions use things like zero(eltype(xs)) for type-stability.

The idea that the forward evaluation is unchanged is often thought of as an axiom of AD, but for discontinuous functions, I think that's another way of saying ε << ζ. Which is a choice. And one that your calculus teacher would disapprove of. The point of evaluating a function with dual numbers is, presumably, to find derivatives, so finding them correctly ought to have a higher priority.

There are other comparisons to think about, for example:

sq2(x) = x>1 ? x^2 : x<1 ? x^2 : one(x)

clamp2(x, lo=0, hi=1) = x>hi ? oftype(x,hi) : x<lo ? oftype(x,lo) : x
clamp3(x, lo=0, hi=1) = x>=hi ? oftype(x,hi) : x<=lo ? oftype(x,lo) : x

[ForwardDiff.derivative(cl, 1.0) for cl in [x->clamp(x,0,1), clamp2, clamp3]] == [1,1,0]
[central_fdm(5, 1)(cl, 1.0) for cl in [x->clamp(x,0,1), clamp2, clamp3]] ≈ [0.5, 0.5, 0.5]

I'm not sure how often simulating x==1 as in sq2(x) happens in the wild. Perhaps from some combination like f(x) = relu(x) + 0.1*relu(-x)?

But clamping parameters to some range is routine. Your calculus teacher would throw an error here, but that's probably not the most helpful response for the computer.

Returning a nonzero derivative here is useful because, if this is some parameter being optimised, it means gradient descent won't get stuck against the wall, when the gradient is away from it. So you can argue that the ability to choose which sub-gradient ForwardDiff will use is a feature. The 0.5 gradient alla FiniteDifferences would also be fine for getting un-stuck, but it's very difficult to picture how ForwardDiff could evaluate both branches, and easy to picture doing so having awful side-effects.

Here is one way to relate the present rule for >(::Dual, ::Real) and >=(::Dual, ::Real) to the finite-everything ζ << ε story. We can say that while the ε-ball overlaps with both sides, the vote from the longer side (longer by about 2ζ) always wins by a hair:

----------(==========1==========)----------  abs(x-1) < ε
---------------------1-(===================  x > 1+ζ
          +++++++++++++..........            gradient votes, clamp2(1.0)

Trying out the above ==(::Dual, ::Real) rule, it looks like the tests of this package all pass, except for the ones explicitly testing such rules. It would be interesting to know if this breaks any other uses in the wild. It would also be interesting to think up other pathological examples, maybe I've missed something important.

Also:

• Another way to talk about this: The problem with prod2 above, and det in support for the determinant function #197, is that they promote accidental zeros of the input to structural zeros. And AD then respects these, and gets the wrong answer. What looked like a simple optimisation when writing for real numbers, has been unintentionally promoted to a constraint on what derivatives are allowed. This is the reverse of the discussion about preserving structural zeros in things like Zygote.gradient(sum∘exp, Diagonal(ones(3)))[1].

• Some other packages get this right, such as TaylorSeries.jl, this C++ code, and this Ruby. Some get it wrong (according to me) like this Mathematica code, this paper with Matlab, and this blog post, although he changed his mind from right to wrong. More mathematical treatments seem to regard the tuple (x,δx) as inheriting == from tuples, i.e. they get it right.

• Similar things were also discussed in Problem with abs #377, where the example is this:

sq3(z) = abs(-z^2)
sq4(z) = abs2(z)
sq5(z) = z^2

[ForwardDiff.derivative(f, 0.0) for f in [sq, sq2, sq3, sq4, sq5]] == [0,0,0,0,0]
[ForwardDiff.hessian(x -> f(x[1]), [0.0])[1] for f in [sq, sq2, sq3, sq4, sq5]] == [2,2,-2,2,2]

A rule was suggested there in which x > y behaves differently for x.value == y.value, breaking such ties by comparing x.partials > y.partials. In the clamp2 example, whether you get stuck against the wall presumably shouldn't depend on whether you minimise loss(x) or maximise -loss(x), so we probably don't want to compare x.partials .> 0 when only x is a dual number. But the rule when both x and y are dual might be worth some more thought.

[mcabbott]

There's another example here: TuringLang/DistributionsAD.jl#23 (comment). Behaviour on tagged ForwardDiff:

using ForwardDiff, LinearAlgebra
A = Matrix(I, 2,2)

ForwardDiff.gradient(x -> sum(x \ [1,3]), A)           # [-1  0; 0 -3]
ForwardDiff.gradient(x -> sum(cholesky(x) \ [1,3]), A) # [-1 -4; 0 -3] 

Aplus2 = [1 0; 0.001 1]  # perturb A[2]
sum(Aplus2 \ [1,3]) - sum(A \ [1,3])                      # -0.001
sum(cholesky(Aplus2) \ [1,3]) - sum(cholesky(A) \ [1,3])  # PosDefException: matrix is not Hermitian

Aplus23 = [1 0.001; 0.001 1]  # perturb A[2] and A[3] together
sum(Aplus23 \ [1,3]) - sum(A \ [1,3])                     # -0.004
sum(cholesky(Aplus23) \ [1,3]) - sum(cholesky(A) \ [1,3]) # -0.004

Perhaps one could claim the cholesky result is particular sub-gradient. But will the algorithm always behave this way? I'm don't know. The one for simple \ is silently wrong --- or at least, it has promoted A to be structurally diagonal. The Jacobians are:

ForwardDiff.jacobian(x -> x \ [1,3], A)
 # -1.0  -0.0  -0.0  -0.0
 # -0.0  -0.0  -0.0  -3.0

ForwardDiff.jacobian(x -> cholesky(x) \ [1,3], A)
 # -1.0  0.0  -3.0   0.0
 #  0.0  0.0  -1.0  -3.0

With the proposed change to ==, the results instead match finite differences. Here \ may choose a different algorithm, but this shouldn't cause a discontinuity of the forward pass.

ForwardDiff.gradient(x -> sum(x \ [1,3]), A)           # [-1 -3; -1 -3]
ForwardDiff.gradient(x -> sum(cholesky(x) \ [1,3]), A) # ERROR: PosDefException

ForwardDiff.jacobian(x -> x \ [1,3], A)
 # -1.0   0.0  -3.0   0.0
 #  0.0  -1.0   0.0  -3.0

Aplus2 = [1 0; 0.001 1]
Aplus2 \ [1, 3] - A \ [1, 3]  # [0, -0.001]

Aplus3 = [1 0.001; 0 1]
Aplus3 \ [1, 3] - A \ [1, 3]  # [-0.003, 0]

The example just above that, in TuringLang/DistributionsAD.jl#23 (comment), also gives an error on this branch, also from cholesky. As does a finite perturbation.

[mateuszbaran]

I didn't fully think it through but wouldn't it break functions like f(x) = x == 0.0 ? 1.0 : sin(x)/x? What if I replaced == 0.0 with iszero or ≈ 0.0?

[mateuszbaran]

Right now if I wanted correct derivatives of such functions up to a certain order I'd just put in the "measure zero branch" a Taylor approximation of the function.

[mcabbott]

This f would indeed change. However, it has problems on a region bigger than one point, and so it would be better anyway for the smoothed branch to be used within some finite interval. For example:

julia> fpi(x) = x==0 ? 1.0 : sin(pi*x) / (pi*x);

julia> fpi(1e-40)
1.0

julia> ForwardDiff.derivative(fpi, 1e-40)
1.2089258196146292e24

This is a variant of #466, for which the solution is

julia> cosc(1e-40) 
-3.2898681336964526e-40

Even if you don't have an exact closed-form derivative, I think that you would usually want to replace f(x) with a constant / polynomial within some small interval, not at a single point. Perhaps always, are there exceptions?

Since the fallback is iszero(x) = x == zero(x) I'd be hesitant to mess with that, but isapprox could be given a different behaviour for dual numbers. At zero it doesn't by default widen the interval:

julia> ForwardDiff.Dual(1e-40,0) ≈ 0
false

julia> nextfloat(0.0) ≈ 0.0
false

[mcabbott]

One more data point is that https://github.com/JuliaPhysics/Measurements.jl gets this right according to my argument above. While its error bars aren't exactly dual numbers, they are a related species. Here's a version of #536 in which mul! does not take the shortcut:

julia> using Measurements

julia> λ = measurement(0, 0.1)
0.0 ± 0.1

julia> iszero(λ)
false

julia> A = measurement.([1 2; 3 4], 0.1);

julia> B = measurement.([5 6; 7 8], 0.1);

julia> A * (B * λ)
2×2 Matrix{Measurement{Float64}}:
 0.0±1.9  0.0±2.2
 0.0±4.3  0.0±5.0

julia> mul!(similar(A), A, B, λ, 0)  # like issue 536
2×2 Matrix{Measurement{Float64}}:
 0.0±1.9  0.0±2.2
 0.0±4.3  0.0±5.0

[ChrisRackauckas]

The idea that the forward evaluation is unchanged is often thought of as an axiom of AD

And it's definitely a false one. It's an axiom for computer scientists, but not for numerical analysts 😅. A nice example is in the space of ODEs. Automatic differentiation is equivalent to solving the expanded ODE known as the forward sensitivity equations, essentially:

u' = f(u,p,t)
d/dt (du/dp) = df/du du/dp + df/dp

Straight automatic differentiation is equivalent to solving the expanded ODE with the adaptive error controls only applying to the first part of the equation u' = f(u,p,t). Are there ODEs for which the second part is unstable when adaptivity is only applied to the first part? Yes. SciML/SciMLSensitivity.jl#273 is a real-world example where this came up. The solution was SciML/SciMLSensitivity.jl#273, i.e. the norm used in the ODE solver has to account for the pseudo-ODEs if you want it to be stable, and so the default norm that is used adds the partials to the primal part.

https://github.com/SciML/DiffEqBase.jl/blob/v6.83.1/src/forwarddiff.jl#L31-L34

This means that solving with ForwardDiff gives different stepping behavior, but if you don't do that, then there will be cases where you have "infinite" derivative because of numerical instability in the derivative calculation even when the primal is stable. So definitely, this axiom does not hold for the realities of numerical computing.

So back to the core of the thread, I definitely agree with you. In fact, DiffEq specializes its interpolation computation in order to work around this kind of issue. Normally it would just pull sol.u[i] if the interpolant matches sol.t[i] directly, but it needs to still use the interpolation if it's a dual number since otherwise the derivative is zero. If it had this epsilon ball definition, the workaround to force Dual numbers to not take the sol.t[i] == t branch would be eliminated.