Logarithm of a sum without underflow

logsum.py

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()

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

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.

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

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

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.