Skip to content

Instantly share code, notes, and snippets.

@mblondel mblondel/logsum.py
Created Apr 13, 2010

Embed
What would you like to do?
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

This comment has been minimized.

Copy link

commented May 16, 2011

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

@mblondel

This comment has been minimized.

Copy link
Owner Author

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

This comment has been minimized.

Copy link

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

This comment has been minimized.

Copy link

commented Jun 3, 2018

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

@timvieira

This comment has been minimized.

Copy link

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
You can’t perform that action at this time.