# Derivatives at infinity issue

**URL:** <https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445>\
**Category:** General Usage\
**Created:** [March 11, 2024, 6:44am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445 "2024-03-11T06:44:15Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 11, 2024, 6:44am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/1 "2024-03-11T06:44:15Z")

</div>

I am trying to implement L’Hôpital’s rule, and my first naive version looks like this:

```julia
using ForwardDiff: derivative
∂(f) = x -> derivative(f,x)
function Hôpital(f,g,a)
    fₐ, gₐ = f(a), g(a)
    any(isnan.((fₐ,gₐ))) && error("fₐ = $fₐ, gₐ = $gₐ")
    r = fₐ / gₐ
    @show fₐ, gₐ, r
    isnan(r) ? Hôpital(∂(f), ∂(g), a) : r
end

Hôpital(x -> x^10 - 1, x -> x^2 - 1,1.0) # 5.0
Hôpital(exp, x -> x,Inf) # Inf
Hôpital(sin, x -> x^2, 0.0) # Inf
Hôpital(log, x -> x, 0) # -Inf, perfect. 
Hôpital(exp, x -> x^2,Inf) # Weirdly doesnt work. 

```

It looks like the second derivative of `exp` at infinity is `NaN`, while it should be `Inf`. Same for `x -> x^2`, while it should be `2`.

is there a way to get these derivatives at infinity correct ?

Note: i am not interested in using a package like Scipy to compute limits, but rather to see how far I can go with standard Julia.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [March 11, 2024, 7:07am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/2 "2024-03-11T07:07:45Z")

</div>

Doesn’t this require symbolic differentiation?

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 11, 2024, 7:09am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/3 "2024-03-11T07:09:31Z")

</div>

Well, as soon as ForwardDiff evaluates the second derivative of `x -> exp(x)` at `x = Inf` to be `Inf` and not `NaN`, this will work. But maybe it will be more efficient with symbolic diff ? I do not know.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 11, 2024, 8:08am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/4 "2024-03-11T08:08:27Z")

</div>

So indeed, this version works a bit better on the examples i could try:

```julia
using Symbolics
@variables t
∂(f) = Symbolics.derivative(f,t,simplify=true)
build(f) = build_function(simplify(f), t, expression=Val{false})
function Hôpital(f,g,a)
    fₜ, gₜ = f(t), g(t)
    L = build(fₜ/gₜ)(a)
    @show fₜ/gₜ
    while isnan(L)
        fₜ, gₜ = ∂(fₜ), ∂(gₜ)
        L = build(fₜ/gₜ)(a)
        @show fₜ/gₜ
    end
    return L
end

Hôpital(x -> x^10 - 1, x -> x^2 - 1,1.0) == 5.0
Hôpital(exp, x -> x^20,Inf) == Inf
Hôpital(sin, x -> x^2, 0.0) == Inf
Hôpital(log, x -> x, 0) == -Inf 
Hôpital(x -> cos(x)-1, x -> x^2, 0) == -1/2
Hôpital(x -> sin(5*x), x -> sin(2*x), 0) == 5/2

```

I wonder if there will be falling cases or if it matches the theorems assumptions correctly…

Also this is pretty slow 😕

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 11, 2024, 9:49am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/5 "2024-03-11T09:49:44Z")

</div>

> [@lrnv](#):
>
> It looks like the second derivative of `exp` at infinity is `NaN`, while it should be `Inf`. Same for `x -> x^2`, while it should be `2`.

I think the correct version of that sentence is

> The derivative at infinity of _the code we use to approximate the derivative of `exp`_ is `NaN`.

Because indeed if we were trying to just differentiate `exp` we would get `Inf` again. So the issue is that the `ForwardDiff.derivative` of `exp` is not exactly the function `exp`, which means ForwardDiff dual numbers presumably don’t work as well with it (due to the lack of custom rule).

This is related to a blog post by @oxinabox I discovered yesterday: [Automatic Differentiation Does Incur Truncation Errors (kinda)](https://www.oxinabox.net/2021/02/08/AD-truncation-error.html)

I don’t have a ready made solution but perhaps it can help you think about the problem

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 11, 2024, 10:02am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/6 "2024-03-11T10:02:34Z")

</div>

You are right on the formulation: the twiced-autodiffed version of julia’s `exp(::Float64)` is not behaving correctly for `Inf`. Also tried `BigFloat`, same issue.

I remember trying out ForwardDiff.jl but also TaylorSeries.jl and TaylorDiff.jl, all three have the same issue.

Reading your link, I think it might be fixable by telling the AD system that `exp'` is `exp`. I wander why `ForwardDiff` is not aware of that already.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 11, 2024, 10:08am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/7 "2024-03-11T10:08:33Z")

</div>

ForwardDiff has special rules for the exponential function applied to a dual number. But it doesn’t have rules for transforming functions into functions (that’s more the realm of symbolic computing).

```julia
julia> using ForwardDiff

julia> ∂(f) = x -> derivative(f, x)
∂ (generic function with 1 method)

julia> exp
exp (generic function with 14 methods)

julia> ∂(exp)
#7 (generic function with 1 method)

```

How is ForwardDiff supposed to know that `#7` has specific derivative behavior?

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 11, 2024, 10:15am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/8 "2024-03-11T10:15:30Z")

</div>

So we are doomed to have

```julia
julia> using ForwardDiff: derivative

julia> derivative(exp,Inf)
Inf

julia> derivative(t -> derivative(exp,t),Inf)
NaN

julia> 

```

? If I get you correctly, this is a dispatch issue : the path for derivating exp correctly exists, but cannot be reached because `t -> derivative(exp,t)` is not the same function.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 11, 2024, 10:17am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/9 "2024-03-11T10:17:31Z")

</div>

Which is why the symbolic approach makes sense.

- a symbolic differentiation system knows that the derivative of exp (as a function) is still exp
- an algorithmic differentiation system only knows how to compute this derivative at a given point

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 11, 2024, 10:19am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/10 "2024-03-11T10:19:16Z")

</div>

But what I do not understand is that there is no chain rule for `t -> derivative(exp,t)`, right, so forwardDiff should just plug a `Dual` into it, and then stepping through the code will encounter the _inside_ `exp(::Dual)` and should do the right thing… But clearly this is not the case.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 11, 2024, 10:38am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/11 "2024-03-11T10:38:58Z")

</div>

Here’s an example that shows you how a `NaN` might arise from this type of code. Not sure this is exactly what goes on inside ForwardDiff but it’s probably close:

```julia
struct MyDual{T}
    x::T
    d::T
end

Base.show(io::IO, d::MyDual) = print(io, "D($(d.x), $(d.d))")

function Base.one(d::MyDual)
    res = MyDual(one(d.x), zero(d.d))
    @info "one($d) = $res"
    return res
end

function Base.:*(d1::MyDual, d2::MyDual)
    res = MyDual(d1.x * d2.x, d1.x * d2.d + d2.x * d1.d)
    @info "$d1 * $d2 = $res"
    return res
end

function Base.exp(d::MyDual)
    res = MyDual(exp(d.x), exp(d.x) * d.d)
    @info "e^$d = $res"
    return res
end

myderivative(f, x) = f(MyDual(x, one(x))).d

exp_der(x) = myderivative(exp, x)

```

```julia
julia> myderivative(exp, Inf)
[ Info: e^D(Inf, 1.0) = D(Inf, Inf)
Inf

julia> myderivative(exp_der, Inf)
[ Info: one(D(Inf, 1.0)) = D(1.0, 0.0)
[ Info: e^D(Inf, 1.0) = D(Inf, Inf)
[ Info: e^D(Inf, 1.0) = D(Inf, Inf)
[ Info: D(Inf, Inf) * D(1.0, 0.0) = D(Inf, NaN)
[ Info: e^D(D(Inf, 1.0), D(1.0, 0.0)) = D(D(Inf, Inf), D(Inf, NaN))
NaN

```

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [March 11, 2024, 10:40am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/12 "2024-03-11T10:40:47Z")

</div>

It is encountering a rule for `exp`, but the problem is that it’s subtracting a term from it multiplied by zero. The actual thing being computed for `exp''(t)` is

```julia
exp_t = exp(t)
1.0 * exp_t + 0.0 * exp_t^2 

```

and when `t = Inf`, then `0.0 * exp_t^2` is `NaN`.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 11, 2024, 10:41am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/13 "2024-03-11T10:41:44Z")

</div>

I think it’s the same multiplication `Inf * 0.0` which you see in my code. I didn’t have the mathematical explanation though!

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [March 11, 2024, 10:52am UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/14 "2024-03-11T10:52:36Z")

</div>

It’s because the AD system is calculating the derivative of the constant value `1.0` to be `0.0`. If you use Zygote.jl though, it has “strong” zeros for constants and so it gets `Inf` here:

```julia-repl
julia> using Zygote: var"'"

julia> exp''(Inf)
Inf

```

and

```julia-repl
julia> function Hôpital(f,g,a)
           fₐ, gₐ = f(a), g(a)
           any(isnan.((fₐ,gₐ))) && error("fₐ = $fₐ, gₐ = $gₐ")
           r = fₐ / gₐ
           @show fₐ, gₐ, r
           isnan(r) ? Hôpital(f', g', a) : r
       end;

julia> Hôpital(x -> x^10 - 1, x -> x^2 - 1,1.0)
(fₐ, gₐ, r) = (0.0, 0.0, NaN)
(fₐ, gₐ, r) = (10.0, 2.0, 5.0)
5.0

julia> Hôpital(exp, x -> x,Inf)
(fₐ, gₐ, r) = (Inf, Inf, NaN)
(fₐ, gₐ, r) = (Inf, 1.0, Inf)
Inf

julia> Hôpital(sin, x -> x^2, 0.0)
(fₐ, gₐ, r) = (0.0, 0.0, NaN)
(fₐ, gₐ, r) = (1.0, 0.0, Inf)
Inf

julia> Hôpital(log, x -> x, 0)
(fₐ, gₐ, r) = (-Inf, 0, -Inf)
-Inf

julia> Hôpital(exp, x -> x^2,Inf) 
(fₐ, gₐ, r) = (Inf, Inf, NaN)
(fₐ, gₐ, r) = (Inf, Inf, NaN)
(fₐ, gₐ, r) = (Inf, 2.0, Inf)
Inf

```

This is something Diffractor.jl _should_ support but currently doesn’t, so I opened an issue here: [Not using "strong zeros" for deriatives of constant values · Issue #278 · JuliaDiff/Diffractor.jl · GitHub](https://github.com/JuliaDiff/Diffractor.jl/issues/278)

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 11, 2024, 6:30pm UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/15 "2024-03-11T18:30:28Z")

</div>

This gives a really short answer to my problem:

```julia
using Zygote: var"'"
function Hôpital(f,g,a)
    r = f(a) / g(a)
    isnan(r) ? Hôpital(f', g', a) : r
end

# All these are fine: 
Hôpital(x -> cos(x)-1 + sin(x^2), x -> x^3, 0) == Inf
Hôpital(x -> x^10 - 1, x -> x^2 - 1,1.0) == 5.0
Hôpital(exp, x -> x^2,Inf) == Inf
Hôpital(sin, x -> x^2, 0.0) == Inf
Hôpital(log, x -> x, 0) == -Inf 
Hôpital(sin, x-> x^2, 0) == Inf
Hôpital(x -> sin(5*x), x -> sin(2*x), 0) == 5/2

```

But, weirdly, it is REAAAAALY slow when the second or third derivative is hit. Might be related to what @oxinabox said on your diffractor issue there [Not using "strong zeros" for deriatives of constant values · Issue #278 · JuliaDiff/Diffractor.jl · GitHub](https://github.com/JuliaDiff/Diffractor.jl/issues/278#issuecomment-1988595891)

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [March 11, 2024, 10:07pm UTC](https://discourse.julialang.org/t/derivatives-at-infinity-issue/111445/16 "2024-03-11T22:07:39Z")

</div>

> [@lrnv](#):
>
> But, weirdly, it is REAAAAALY slow when the second or third derivative is hit.

Not weird at all, Zygote is very much a bad system for higher order AD, so bad that Diffractor.jl was designed specifically so that it’d compose better.

Basically if you do `n`th order derivatives, then Zygote will produce an amount of code that goes like the exponential or factorial of `n`, most of which is trivial junk, and making an insane amount of work for the compiler to sift through.