# Trying to make my robust logistic function differentiable as part of AD, getting errors so far

**URL:** <https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838>\
**Category:** General Usage\
**Tags:** question\
**Created:** [May 5, 2020, 9:46pm UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838 "2020-05-05T21:46:22Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [May 5, 2020, 9:46pm UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/1 "2020-05-05T21:46:22Z")

</div>

I found that `StatsFuns.logistic` occasionally had issues with really extreme values and so have this more robust version of the logistic function defined as below:

```julia
const LOGISTIC_64 = log(Float64(1)/eps(Float64) - Float64(1))
const LOGISTIC_32 = log(Float32(1)/eps(Float32) - Float32(1))
"""
Return the logistic function computed in a numerically stable way:
``logistic(x) = 1/(1+exp(-x))``
"""
function logistic(x::Float64)
	x > LOGISTIC_64 && return one(x)
	x < -LOGISTIC_64 && return zero(x)
	return one(x) / (one(x) + exp(-x))
end
function logistic(x::Float32)
	x > LOGISTIC_32 && return one(x)
	x < -LOGISTIC_32 && return zero(x)
	return one(x) / (one(x) + exp(-x))
end
logistic(x) = logistic(float(x))

```

I am using it as part of some densities for MCMC sampling and Bayesian inference. I’d like to be able to apply AD techniques to some functions that involve the logistic function as implemented above, but when I do end up getting the following error:

```julia
ERROR: StackOverflowError:
Stacktrace:
 [1] logistic(::ForwardDiff.Dual{ForwardDiff.Tag{var"#57#58",Float64},Float64,4}) at /Users/harrisonwilde/Library/Mobile Documents/com~apple~CloudDocs/PhD/Synthetic Data/src/creditcard/distrib_utils.jl:21 (repeats 79984 times)

```

Presumably as it repeatedly calls the bottom line over and over because calling `float(x)` where `x` is of type `ForwardDiff.Dual{ForwardDiff.Tag{var"#57#58",Float64},Float64,4}` must not work. How can I write my functions so they are “differentiable” in this setting, any guidance would be much appreciated!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2020, 9:53pm UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/2 "2020-05-05T21:53:37Z")

</div>

Just write your function like this

```julia
function LOGISTIC(T) 
    log(one(T)/eps(T) - one(T))
end

function logistic(x::T) where {T}
    LOGISTIC_T = LOGISTIC(T)
    x > LOGISTIC_T && return one(x)
    x < -LOGISTIC_T && return zero(x)
    return one(x) / (one(x) + exp(-x)) end

```

This way, your function will work for any type supporting `eps`, `one`, and `exp`. Note that due to specialization, this should be just as fast as your implementation.

---

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [May 5, 2020, 11:39pm UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/3 "2020-05-05T23:39:11Z")

</div>

> [@Oscar\_Smith](#):
>
> function LOGISTIC(T) log(one(T)/eps(T) - one(T)) end function logistic(x::T) LOGISTIC\_T = LOGISTIC(T) x \> LOGISTIC\_T && return one(x) x \< -LOGISTIC\_T && return zero(x) return one(x) / (one(x) + exp(-x)) end

Thank you, that makes a lot of sense, I think you forgot to add a `where {T}` on the function definition though right?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2020, 11:40pm UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/4 "2020-05-05T23:40:14Z")

</div>

Fixed, Good point!

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [May 6, 2020, 1:04am UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/5 "2020-05-06T01:04:35Z")

</div>

This should be fixed in StatsFuns.jl directly. I’ll try to submit a patch tonight since there’s a few different issues worth fixing here:

1. It’s better to use `exp(x) / (exp(x) + 1)` here to handle subnormals properly.

Compare and contrast the implementations in the region that generates subnormal outputs:

```julia
let x = -740.0
    exp(x) / (exp(x) + 1)
end
# 4.2e-322

let x = -740.0
    1 / (1 + exp(-x))
end
# 0.0

```

1. You might want to use `ifelse` instead of `&&` since minimizing branches increases the chances SIMD will work. It may not be sufficient to get SIMD to kick in, but I would guess it will work out better.

2. A potentially bigger issue is that, if you’re ever computing `log(logistic(x))`, you should use the specialized technique from Equation 10 on page 7 of [this paper](https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf) rather than compute `logistic(x)` and then apply `log` to the result.

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [May 6, 2020, 6:13am UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/6 "2020-05-06T06:13:43Z")

</div>

> [@johnmyleswhite](#):
>
> Compare and contrast the implementations in the region that generates subnormal outputs:

I’m not sure I agree with this example, can you clarify why you think the first one is preferable? Also the thresholding here makes sure that the answer would be zero anyway but in a predictable way since the threshold for Float64 is around `±36`.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 6, 2020, 6:26am UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/7 "2020-05-06T06:26:59Z")

</div>

The OP formula doesn’t deal with subnormals properly, but it does have the benefit of producing `1.0` instead of `NaN` for larger positive numbers. Possibly the best of both worlds would be to switch which you use based on sign of x.

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [May 6, 2020, 10:07am UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/8 "2020-05-06T10:07:47Z")

</div>

Run the calculation in higher precision and then truncate to `Float64`:

```julia
Float64(1 / (1 + exp(-big(-740.0))))
# 4.2e-322

```

Over all intervals I’ve tested, the first implementation is closer to the higher precision result – which is to say that it seems **the first strictly dominates the second in precision**.

I think it’s best not to switch based on the sign of `x` since you can simply guard against `NaN`, but the first implementation is more precise for large positive inputs in the pre-overflow region. See the PR I made for what I think is the proper way to do this: [Improve accuracy of logistic by johnmyleswhite · Pull Request #94 · JuliaStats/StatsFuns.jl · GitHub](https://github.com/JuliaStats/StatsFuns.jl/pull/94)

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [May 6, 2020, 10:21am UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/9 "2020-05-06T10:21:13Z")

</div>

I think my question was more to do with “what would you do with that 4.2e-322” ? I think there’s a point where it makes sense to truncate results to avoid nasty propagation of errors. In the context of special funs, use of `big` etc, I could see users willing to use your version but usually this function is used in ML/Stats where it makes little sense to keep track of these values, just like it would make little sense (IMO) to train a neural net in big precision.

Anyway I think it’s just a matter of perspective, I agree with your point about `big` precision, my perspective is that if you truncate (which I believe can often make sense) the point is irrelevant.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [May 6, 2020, 10:45am UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/10 "2020-05-06T10:45:52Z")

</div>

Tangentially related, subnormals can be really slow to compute with ([Subnormal number - Wikipedia](https://en.wikipedia.org/wiki/Denormal_number#Performance_issues), [50x speed difference in gemv for different values in vector](https://discourse.julialang.org/t/50x-speed-difference-in-gemv-for-different-values-in-vector/2755), [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/#Treat-Subnormal-Numbers-as-Zeros-1))

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [May 6, 2020, 11:03am UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/11 "2020-05-06T11:03:28Z")

</div>

That’s fair: there’s a tradeoff here between precision and performance.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 6, 2020, 4:08pm UTC](https://discourse.julialang.org/t/trying-to-make-my-robust-logistic-function-differentiable-as-part-of-ad-getting-errors-so-far/38838/12 "2020-05-06T16:08:26Z")

</div>

> [@HarrisonWilde](#):
>
> ```julia
> const LOGISTIC_64 = log(Float64(1)/eps(Float64) - Float64(1)) 
> const LOGISTIC_32 = log(Float32(1)/eps(Float32) - Float32(1))
> ...
> one(x) / (one(x) + exp(-x))
> 
> ```

I don’t really have any opinion on the numerical issues, but there seems to be a lot of unnecessary typing here. Due to type promotion you can just as well write

```julia
const LOGISTIC_64 = log(1 / eps(Float64) - 1) 
const LOGISTIC_32 = log(1 / eps(Float32) - 1)
...
1 / (1 + exp(-x))

```

If you really need literal `Float64` or `Float32` you can write

```julia
1.0 # Float64
1e0 # also Float64
1f0 # Float32

```

Integers are really nice, you can use them in most places and they will promote to the correct precision. This too, is redundant:

```julia
logistic(x) = logistic(float(x))

```

since the division inside the `logistic` function will promote the values to appropriate floats.