Skip to content

Instantly share code, notes, and snippets.

Created April 13, 2010 06:39
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mblondel/364369 to your computer and use it in GitHub Desktop.
Save mblondel/364369 to your computer and use it in GitHub Desktop.
Logarithm of a sum without underflow
import numpy as np
def _logsum(logx, logy):
Return log(x+y), avoiding arithmetic underflow/overflow.
logx: log(x)
logy: log(y)
x + y = e^logx + e^logy
= e^logx (1 + e^(logy-logx))
log(x+y) = logx + log(1 + e^(logy-logx)) (1)
log(x+y) = logy + log(1 + e^(logx-logy)) (2)
The computation of the exponential overflows earlier and is less precise
for big values than for small values. Due to the presence of logy-logx
(resp. logx-logy), (1) is preferred when logx > logy and (2) is preferred
if logx > logy:
return logx + np.log(1 + np.exp(logy-logx))
return logy + np.log(1 + np.exp(logx-logy))
logsum_ufunc is a numpy ufunc (universal function) and as a result contains
reduce, accumulate, reduceat and outer.
logsum_ufunc = np.frompyfunc(_logsum, 2, 1)
logsum(loga, axis=0, dtype=None, out=None)
Take the log of the sum of array elements over a given axis.
For example, for an array a=[a_1,...,a_N], it returns \log \sum_n a_n.
loga: numpy.log(a)
axis: The axis along which to apply the log sum.
dtype: The type used to represent the intermediate results.
out: A location into which the result is stored.
logsum(1darray) => scalar
logsum(2darray, axis=0) => 1darray
logsum = logsum_ufunc.reduce
# Unit-tests...
if __name__ == "__main__":
import unittest
class Test(unittest.TestCase):
def test_1d(self):
inp = np.arange(1,10)
out = logsum(np.log(inp))
expected = np.log(np.sum(inp))
self.assertAlmostEquals(out, expected)
def test_2d(self):
for a in (0,1):
inp = np.arange(1,10).reshape(3,3)
out = logsum(np.log(inp), axis=a)
expected = np.log(np.sum(inp, axis=a))
#self.assertTrue(np.allclose(out, expected))
for i in range(3):
self.assertAlmostEquals(out[i], expected[i])
Copy link

What's the advantage of your _logsum over numpy.logaddexp?

Copy link

I didn't know that numpy had this function at that time. I know two tricks to compute log sums but I don't remember which one does numpy use.

Copy link

Ok. I was summing a lot of logprobs today, so I wondered if your method would buy me anything that logaddexp wouldn't :)

Copy link

Protip: use log1p(x) instead of log(1+x), but it only makes a difference for small x.

Copy link

timvieira commented Jun 3, 2018

If you really want to go down the rabbit hole, check out this article. Despite the title, it does cover log(1+exp(x)), which is the crux of logaddexp.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment