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

**URL:** <https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229>\
**Category:** Numerics\
**Created:** [January 1, 2024, 2:16pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229 "2024-01-01T14:16:26Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [January 1, 2024, 2:16pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/1 "2024-01-01T14:16:26Z")

</div>

It occurred to me this morning that `logsumexp` is associative and commutative and can be calculated with a simple reduce. We could write

```julia
logsumexp(x) = reduce(logaddexp,x)

```

Isn’t that pretty!

So how does that compare to the really fancy implementation in LogExpFunctions.jl? The reduction is 5x slower, but surprisingly fast, equivalent to a different 1-pass algorithm than discovered by [Sebastian Nowozin: Streaming Log-sum-exp Computation](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)

```julia
using LogExpFunctions
using BenchmarkTools

x = rand(10^7)

@btime a1 = logsumexp(x)
  70.892 ms (1 allocation: 16 bytes)
16.659565226722435

@btime a2 = reduce(logaddexp, x)
 310.882 ms (1 allocation: 16 bytes)
16.659565226722393

```

But also: Interesting! a bit of numerical difference at 10^(-14)! Who is more accurate? I looked at the same algorithms at Float32:

```julia
y = Float32.(x)
b1 = logsumexp(y)
b2 = reduce(logaddexp, y)
b1-Float32(a1)
0.013198853f0
b2-Float32(a1)
0.0f0

```

The reduction is much more numerically stable! For larger vectors of rand the logsumexp gets arbitrarily bad, albeit at somewhat enormous vectors.

In a sense, the rand case is _particularly_ bad because the maximum in this example had a lot of “multiplicity” many values in `x` near the maximum contribute a non-trivial amount to the final result. So subtracting off the maximum is not enough, it’s actually better to subtract off the larger of the _running sum_ or the next number. That’s what the reduction does.

In fact, if you know a little physics, the log-multiplicity is exactly the entropy in the stat-mech sense (-Σ\_i p\_i log p\_i), the effective number of states contributing to the free energy (logsumexp)

```julia
logp = x .- logsumexp(x)
logmplogp = logp .+ log.(-logp)
logS = logsumexp(logmplogp)
S = exp(logS)
16.0774514705695

```

quite near the logsumexp(x), so logsumexp(x) here is clearly “entropy driven” rather than “softmax driven”

Anyway, there are clearly some tradeoffs due to the vectorial nature of cpus so I’m not sure if this an issue. I wouldn’t neccesarily want to force a 5x slowdown on the most expensive calculation in many algorithms for accuracy that may not be important. Is it possible that the reduction version can be made faster with whatever magic is being used on [logsumexp](https://github.com/JuliaStats/LogExpFunctions.jl/blob/master/src/logsumexp.jl)? Should I open an issue?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 1, 2024, 3:09pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/2 "2024-01-01T15:09:03Z")

</div>

> [@cdawg](#):
>
> `logsumexp(x) = reduce(logaddexp,x)`

Where is `logaddexp` defined?

_Update_ - it’s [here](https://github.com/JuliaStats/LogExpFunctions.jl/blob/8b28b3a3c4e55b420af2442a2e5462dbc1ceae56/src/basicfuns.jl#L350-L363) in LogExpFunctions.jl.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [January 1, 2024, 3:24pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/3 "2024-01-01T15:24:32Z")

</div>

```julia
julia> using BenchmarkTools, LogExpFunctions

julia> x = rand(10^7)

julia> @btime a1 = logsumexp($x)
  36.006 ms (0 allocations: 0 bytes)
16.659561175877272

julia> @btime a2 = log(sum(exp,$x))
  28.252 ms (0 allocations: 0 bytes)
16.659561175877275

julia> f(x) = reduce(logaddexp,x)
f (generic function with 1 method)

julia> @btime a3 = f($x)
  240.178 ms (0 allocations: 0 bytes)
16.659561175877275

```

I think that this simple example must not illustrate the benefits of `logsumexp` since `log(sum(exp,x))` is faster. Likely you need different inputs which might induce overflow/underflow without care.

---

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [January 1, 2024, 3:46pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/4 "2024-01-01T15:46:39Z")

</div>

`x = rand(10^7)` works with log(sum(exp.(x))) because x is confined to the nearly-linear part of exp(x), so there’s nothing to worry about.

But, say x = 1000 .+ rand(10^7). Then `log(sum(exp.(x)))` doesn’t work, but both `reduce(logaddexp, x)` and `logsumexp(x)` do. In the space of algorithms that DO work, reduce(logaddexp,x) seems to be more accurate for large vectors.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [January 1, 2024, 4:29pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/5 "2024-01-01T16:29:47Z")

</div>

> [@cdawg](#):
>
> ```julia
> y = Float32.(x)
> b1 = logsumexp(y)
> b2 = reduce(logaddexp, y)
> b1-Float32(a1)
> 0.013198853f0
> b2-Float32(a1)
> 0.0f0
> 
> ```

But is this the right way to check accuracy? Converting `Float64` to `Float32` changes the values and _should_ lead to some difference. Isn’t it better to generate the random values as `Float32` initially and then compare reduction in both precisions?

---

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [January 1, 2024, 5:09pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/6 "2024-01-01T17:09:29Z")

</div>

sure, you’re right, it’s more intuitive to go the other way and promote and then check, and compare to eps:

```julia
y = rand(Float32,10^7)
x = Float64.(y)
a1 = logsumexp(x)
a2 = reduce(logaddexp, x)
b1 = logsumexp(y)
b2 = reduce(logaddexp, y)
(b1-Float32(a1))/eps(Float32)
111040.0f0
(b2-Float32(a1))/eps(Float32)
16.0f0

```

but it amounts to the same thing, since the precision of Float 64 is so much higher.

It seemed intuitive to me to check for accuracy by checking if the Float32 conversion “commutes” with the operation, but I don’t know what the standard method is as I’m not a numerics person.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 1, 2024, 5:47pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/7 "2024-01-01T17:47:26Z")

</div>

To measure accuracy, the results can be compared against:

```julia
c1 = logsumexp(BigFloat.(x))

```

(slow for sure, but a good baseline for accuracy).

Another faster performance/accuracy tradeoff Pareto frontier sweet-spot is:

```julia
c2 = log(sum(exp, x; init=big"0.0"))

```

Another attempt at increasing accuracy, does so slightly with this `x`, but depends on the distribution of `x` values:

```julia
c3 = log(sum(expm1, x; init=big"0.0")+length(x))

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 1, 2024, 6:46pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/8 "2024-01-01T18:46:21Z")

</div>

> [@Dan](#):
>
> Another faster performance/accuracy tradeoff Pareto frontier sweet-spot is:
> 
> `c2 = log(sum(exp, x; init=big"0.0"))`

This is susceptible to spurious overflow if the `x` elements can exceed 710.

---

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [January 1, 2024, 6:57pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/9 "2024-01-01T18:57:53Z")

</div>

```julia
x = rand(10^9)
a1 = logsumexp(x)
a2 = reduce(logaddexp, x)
(a1-a2)/eps(a1)
  306.0
# is only around 25 for 10^8 and 16 for 10^7

```

I have to say the error at the Float64 level is a lot better behaved, presumably the numerics here are super well tested. I do think the reduce version really is more stable.

I wouldn’t expect better guarantees at this accuracy. But the situation really is weird at Float16 / Float32 level. I guess not the first time a FloatNot64 behaves weird?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 1, 2024, 7:02pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/10 "2024-01-01T19:02:31Z")

</div>

> [@cdawg](#):
>
> In the space of algorithms that DO work, reduce(logaddexp,x) seems to be more accurate for large vectors.

I think that’s because this uses [pairwise summation](https://en.wikipedia.org/wiki/Pairwise_summation), whereas LogExpFunctions.jl uses [a “one-pass” algorithm with an `init` argument](https://github.com/JuliaStats/LogExpFunctions.jl/blob/master/src/logsumexp.jl#L81-L107), 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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [January 1, 2024, 7:22pm UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/11 "2024-01-01T19:22:06Z")

</div>

I see. I think you’re right and the entropy explanation I gave is bogus, this is probably pure pairwise enhancement.

I just wouldn’t condition on AbstractArray since I often give it generators. I also worry about summing say [0, 10^-16], you might run into places where log1p is better and logaddexp does that near the limit of what’s possible with Float64’s.

I just have a few more questions:

- Should I open an issue with LogExpFunctions?
- Any idea what is slowing down the reduce(logaddexp, x) here?
- Does it make sense that the errors are so much worse on the one-pass for Float16 and Float32?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 2, 2024, 4:01am UTC](https://discourse.julialang.org/t/reduce-logaddexp-x-has-arbitrarily-better-accuracy-than-logsumexp-x/108229/12 "2024-01-02T04:01:31Z")

</div>

> [@cdawg](#):
>
> Should I open an issue with LogExpFunctions?

Couldn’t hurt.

> [@cdawg](#):
>
> Does it make sense that the errors are so much worse on the one-pass for Float16 and Float32?

You notice roundoff accumulation much more quickly with lower precisions.

> [@cdawg](#):
>
> I just wouldn’t condition on AbstractArray since I often give it generators.

Doing pairwise summation is a bit trickier when you don’t have random access: [use pairwise order for mapreduce on arbitrary iterators by stevengj · Pull Request #52397 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/52397)

You could also use Kahan/compensated summation, maybe in blocks to improve performance.