10
$\begingroup$

Most math libraries have a number of versions of logarithm functions. Most of the time we assume them to be perfect, but actually quite a lot of them just offer a certain number of digits of precision.

For some functions, there are numerically more stable variants. For example Fortran, R, Java and C both have Math.log1p, for computing log(1.0+x) (which offers higher precision for small values of x), and the counterpart expm1. Here the numeric problems arise from a loss in precision - if x is really small, 1.0 + x loses digits in order to preserve the 1 at the beginning.

I've seen such functions for increased precision in a number of situations. This seems to be quite common whenever you are implementing distribution functions (Gamma, Beta, Poisson etc.) with high numeric precision. For example the Gamma function seems to be most of the time used as logGamma. In general, going to "logspace" can improve precision a lot, and so R seems to have a "logspace" flag on most functions.

Another example, in R, there exists log1mexp for log(1 - exp(p)): http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf

I've been playing around with entropy and information theoretic measures. A very common term there is

p * -log(p)

where usually, one would want the base of the logarithm to be 2, not e; but just as often this is only a linear factor, and you can as well use the natural logarithm (so this is not of key importance to me). Anyway, do you know if there is a faster / more direct / more precise way of computing this term? I'm having it all over the place, so it could really pay off to make it a bit more precise and fast (save me the usual "premature optimization" stuff, thanks).

I don't see any obvious reason that would cause a loss in precision. So I'm mostly interested if there is any nice trick to speed up this computation. That maybe even saves me treating the p=0 corner case (which sensibly is 0, although log(0) does not exist) or gives me base 2 for free (although a single multiplication with a constant obviously is not killer expensive). Thanks.

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ If you're worried about over/underflow, note that because $m \sim 10^{-308} \leq p \leq M \sim 10^{-308}$ (double precision), $|\log p|$ will be at most $\approx 700$, so it can't ever be very large. $p \log p$ (0 when $p=0$) should be perfectly fine. As a side note, fast formulas are can be different from accurate formulas, so sometimes you can't get both. $\endgroup$
    – Kirill
    Commented Nov 9, 2012 at 18:07
  • $\begingroup$ R comes with a log2 function which depending on your OS can be a simple wrapper around log/log(2) or make use of the fact that C99 added a log2 function. $\endgroup$
    – anonymous
    Commented Aug 31, 2013 at 0:38

1 Answer 1

Reset to default
10
$\begingroup$

$p \log p$ won't suffer from precision loss anywhere in $[0,1]$, and won't suffer from exponent overflow near $0$ either. Thus, the fast, accurate way is

p ? p * -log(p) : 0
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.