Instantly share code, notes, and snippets.

# mblondel/logsum.py Created Apr 13, 2010

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) Rationale: x + y = e^logx + e^logy = e^logx (1 + e^(logy-logx)) log(x+y) = logx + log(1 + e^(logy-logx)) (1) Likewise, 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 otherwise. """ if logx > logy: return logx + np.log(1 + np.exp(logy-logx)) else: 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. Examples: 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]) unittest.main()

### larsmans commented May 16, 2011

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

### mblondel commented May 16, 2011

 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.

### larsmans commented May 16, 2011

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

### timvieira commented Jun 3, 2018

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

### timvieira commented Jun 3, 2018 • edited

 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.