25
$\begingroup$

I am writing an algorithm where, given a model, I compute likelihoods for a list of datasets and then need to normalize (to probability) each one of the likelihood. So something like [0.00043, 0.00004, 0.00321] might be converted to may be like [0.2, 0.03, 0.77].

My problem is that the log likelihoods, I am working with, are quite small (for instance, in log space, the values are like -269647.432, -231444.981 etc). In my C++ code, when I try to add two of them (by taking their exponent) I get an answer of "Inf". I tried to add them in log-space (Summation/Subtraction of log), but again stumbled upon the same problem.

Can anybody share his/her expert opinion on this?

$\endgroup$
  • $\begingroup$ When you used the functions you pointed to involving $\log(1+\dotsc)$, did you use the log1p function in your language? This uses the Taylor expansion around 1. $\endgroup$ – Neil G Aug 6 '13 at 4:46
  • 1
    $\begingroup$ See also some earlier, related discussion here $\endgroup$ – Glen_b May 3 '15 at 23:43
36
$\begingroup$

Subtract the maximum logarithm from all logs. Throw away all results that are so negative they will underflow the exponential. (Their likelihoods are, for all practical purposes, zero.)

Indeed, if you want a relative precision of $\epsilon$ (such as $\epsilon = 10^{-d}$ for $d$ digits of precision) and you have $n$ likelihoods, throw away any result less than the logarithm of $\epsilon/n$. Then proceed as usual to exponentiate the resulting values and divide each one by the sum of all the exponentials.

For those who like formulas, let the logarithms be $\lambda_1, \lambda_2, \ldots, \lambda_n$ with $\lambda_n = \max(\lambda_i)$. For logarithms to the base $b\gt 1$, define

$$\alpha_i = \cases{ b^{\lambda_i - \lambda_n}, \lambda_i - \lambda_n \ge \log(\epsilon)-\log(n) \\ 0\quad \text{otherwise}.}$$

The normalized likelihoods equal $\alpha_i / \sum_{j=1}^n \alpha_j$, $i = 1, 2, \ldots, n.$ This works because replacing all of the otherwise underflowing $\alpha_i$ by zero makes a total error of at most $(n-1)\epsilon/n\lt \epsilon$ whereas, because $\alpha_n=b^{\lambda_n-\lambda_n}=b^0=1$ and all $\alpha_i$ are non-negative, the denominator $A = \sum_j \alpha_j \ge 1$, whence the total relative error due to the zero-replacement rule is strictly smaller than $\left((n-1)\epsilon/n \right) / A \lt \epsilon$, as desired.

To avoid too much rounding error, compute the sum starting with the smallest values of the $\alpha_i$. This will be done automatically when the $\lambda_i$ are first sorted in increasing order. This is a consideration only for very large $n$.

BTW, this prescription assumed the base of the logs is greater than $1$. For bases $b$ less than $1$, first negate all the logs and proceed as if the base were equal to $1/b$.


Example

Let there be three values with logarithms (natural logs, say) equal to $-269647.432,$ $-231444.981,$ and $-231444.699.$ The last is the largest; subtracting it from each value gives $-38202.733,$ $-0.282,$ and $0.$

Suppose you would like precision comparable to IEEE doubles (about 16 decimal digits), so that $\epsilon=10^{-16}$ and $n=3$. (You can't actually achieve this precision, because $-0.282$ is given only to three significant figures, but that's ok: we're only throwing away values that are guaranteed not to affect the better of the precision you want and the precision you actually have.) Compute $\log(\epsilon/n)$ = $\log(10^{-16}) - \log(3)$ = $-37.93997.$ The first of the three differences, $-38202.733,$ is less than this, so throw it away, leaving just $-0.282$ and $0.$ Exponentiating them gives $\exp(-0.282) = 0.754$ and $\exp(0)=1$ (of course). The normalized values are--in order--$0$ for the one you threw away, $0.754 / (1 + 0.754) = 0.430$, and $1/(1+0.754)=0.570$.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ This is brilliant - so simple, and so obvious in hindsight. @Ikram, please mark this as the correct answer! (unless of course you have something better, in which case please do share) $\endgroup$ – zelanix Jul 14 '16 at 12:29
  • 2
    $\begingroup$ @whuber do we even need to throw away $-38202.733$? Exponentiating that gives us zero anyway, and so it won't contribute to the sum. $\endgroup$ – Taylor Sep 15 '17 at 2:33
  • $\begingroup$ @Taylor exponentiating such a number yield underflow problems. $\endgroup$ – Aku Aug 5 at 17:29
  • $\begingroup$ @ArturJ I suspect Taylor might have been doing the computation in R or some other platform that replaces such underflows with zeros. That behavior is convenient but dangerous: you really need to understand what's going on in your algorithm before you have any right to rely on it. $\endgroup$ – whuber Aug 5 at 18:31
  • 1
    $\begingroup$ good point whuber! @ArturJ as far as I can tell, IEEE 754 doesn't say anything about what to do in the case of underflow, just that it has to be signaled: "[t]he occurrence of an underflow may set a ("sticky") status bit, raise an exception, at the hardware level generate an interrupt, or may cause some combination of these effects." I shouldn't have been relying on the rounding $\endgroup$ – Taylor Aug 6 at 14:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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