log1pexp(x) when x < -37

[cossio]

I was reading R's implementation of log1pexp and they include one additional branch (see https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf, equation 10).

Currently StatsFuns does this:

https://github.com/JuliaStats/StatsFuns.jl/blob/ae40bc1c6ee63a5624ba6cd168da5bf43f473e59/src/basicfuns.jl#L66

Whereas including the additional branch it would be:

log1pexp_v2(x::AbstractFloat) = x ≤ -37. ? exp(x) : x ≤ 18. ? log1p(exp(x)) : x ≤ 33.3 ? x + exp(-x) : x

Apparently including the x < -37. check reduces the function's computation time by 25% (because it saves the call to log1p I guess). But I am no expert doing this kind of performance tests so please can someone else more knowledgeable check this? Here is the test I did:

using BenchmarkTools
@btime for x = -100. : 1. : 100., _ = 1 : 100
    log1pexp_v2(x)
end   #  310.471 μs (0 allocations: 0 bytes)
@btime for x = -100. : 1. : 100., _ = 1 : 100
    log1pexp(x)
end  #   423.671 μs (0 allocations: 0 bytes)

Both log1pexp_v2 and log1pexp give numerically identical results.

If this is worth it I can prepare a PR.

Update: Actually it seems much of the time is due to the oftype(exp(-x), x) bit in the final branch. I think specializing on AbstractFloat would take care of that. Even after this change the x ≤ -37. check is faster.

[dmbates]

I can see how your version will be faster because it avoids the exp(-x) in oftype(exp(-x), x). You may consider adding

log1pexp_v2(x::Number) = log1pexp_v2(float(x))

to get coverage for integer arguments.

[cossio]

I would do log1pexp_v2(x::Real) = log1pexp_v2(float(x)), since we don't want to cover complex arguments.

[dmbates]

You're right. That's better.

[cossio]

Should I prepare a PR on this?

[dmbates]

Yes, please.

[cossio]

@dmbates Please see JuliaStats/StatsFuns.jl#70

[Uroc327]

If I may ask, what's the current state of this?

[devmotion]

I don't know, I just transferred the issue when log1pexp was moved from StatsFuns to LogExpFunctions. @cossio, do you plan to open a PR to LogExpFunctions that preferably also addresses the comments in JuliaStats/StatsFuns.jl#70?

[tpapp]

JuliaStats/StatsFuns.jl#70 was just closed, but a PR here is welcome.

[cossio]

I'll do a PR right after #35 is done.

[tpapp]

Great, thanks!