Reduce(logaddexp,x) has arbitrarily better accuracy than logsumexp(x)

Specific Domains Numerics
10

I think that’s because this uses pairwise summation, whereas LogExpFunctions.jl uses a “one-pass” algorithm with an init argument, which ends up calling foldl (i.e. naive summation).

Although it should be possible to re-implement the one-pass algorithm in a pairwise fashion, it seems a lot easier just to use the two-pass algorithm, which only seems to be about 20% slower than the one-pass algorithm and benefits from pairwise summation:

function logsumexp_2pass(x::AbstractArray)
    isempty(x) && return log(sum(exp, x)) # -Inf
    max = maximum(x)
    return max + log(sum(x -> exp(x - max), x))
end

The fact that Julia’s mapreduce uses pairwise associativity in certain cases is a two-edged sword — it makes it more accurate without losing performance, but then people are continually surprised when other cases are less accurate.

1 Like