Problem with `abs`

[briochemc]

Hi there,

I'm not sure this is intentional or if it is for numerical efficiency reasons, and I may be out of my league discussing this, but I figured I should point this out. When trying to copy the DualNumbers.jl structure into the HyperDualNumbers.jl structure, I noticed that there is a problem with how abs is dealt with. I believe the problem is that abs(x) returns x if x >= 0 (in DualNumbers.jl, HyperDualNumbers.jl, and I believe in ForwardDiff.jl's dual implementation).

Here is a small example where this causes problems. I use ForwardDiff.hessian to evaluate the second derivative of abs(-x^2), which should be the same as the second derivative of x^2, which is always 2:

julia> using ForwardDiff

julia> ForwardDiff.hessian(x -> abs(-x[1]^2), [0.0])
1×1 Array{Float64,2}:
 -2.0

(I used hessian and a single element Vector because I did not know how to get to second derivatives of scalars using ForwardDiff.)

I come with a suggestion: What about dealing with abs in the "physical" way, that is thinking of the absolute value of some x plus or minus some very small term ε. In other words, treating it recursively via something like:

function myabs(d::Dual)
    x, y = value(d), epsilon(d)
    if x > 0
        return d
    elseif x < 0
        return -d
    else
        if y ≥ 0
            return d
        else
            return -d
        end
    end
end

For hyper duals, something like

function myabs(h::Hyper)
    x, y, z, w = value(h), epsilon1(h), epsilon2(h), epsilon12(h)
    if x > 0
        return h
    elseif x < 0
        return -h
    else
        if y+z > 0
            return h
        elseif y+z < 0
            return -h
        else
            if z ≥ 0
                return h
            else z < 0
                return -h
            end
        end
    end
end

Now I know this is not DualNumbers.jl or HyperDualNumbers.jl, to which the above example suggestions apply, but I haven't really figured out how ForwardDiff.jl implements hyperduals so this was the best way I could convey my suggestions. I also believe that it is very likely that one can code the suggested abs in a much more efficient way than my if statements allow. But regardless of that, is that something worse paying some attention?

[KristofferC]

julia> using ForwardDiff

julia> ForwardDiff.hessian(x -> abs(-x[1]^2), [0.0])
1×1 Matrix{Float64}:
 2.0

So I guess this is fixed.