weighted logsumexp

[mileslucas]

• add weighted logsumexp function
• bump version 0.9.5 => 0.9.6

I ran into a use case where I needed a weighted sum in the logsumexp function when copying some python code.

There are no performance regressions in the new version, I also removed the TODO regarding some special-casing on the value u because without it the performance is much worse (30 ns -> 1 us) for cases where dims=:.

Here are the benchmarks I did to make sure there's still a super-fast no allocation for the simplest case

julia> @benchmark logsumexp($([1.0, 2.0]))
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     32.608 ns (0.00% GC)
  median time:      33.277 ns (0.00% GC)
  mean time:        34.576 ns (0.00% GC)
  maximum time:     101.738 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     993

julia> @benchmark logsumexp($([1.0, 2.0]), $([-0.5, 0.5]))
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     36.436 ns (0.00% GC)
  median time:      37.050 ns (0.00% GC)
  mean time:        38.228 ns (0.00% GC)
  maximum time:     150.541 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     992

[bkamins]

in the docstring you say W is an array, but there is no type restriction in the signature. I think it would be cleaner to have it especially that in your code W and X do not have to have matching dimensions, and when do not then zip(X, W) will not iterate every element of X in corner cases.

[mileslucas]

added words and examples

[nalimilan]

Can you just restrict this to ::AbstractArray? I don't see the point of allowing any other type (and we don't in StatsBase). That also allows simplifying the docs. I also don't like allowing broadcasting an array of any dimensions, as currently the weighted sum over dimensions in StatsBase does something different (it applies the weights over the dimension which is being summed over). And undefined behavior is not something we generally accept in Julia packages. So better require the two arrays to have the same size for now.

[mileslucas]

@nalimilan I didn't want to over-specialize because then I couldn't use a call like logsumexp([1, 2], (0.5, 0.3)). This is useful for me because my weights are static in the calculations I need so I prefer to use a Tuple.

As for the making the sizes the same, that would disallow broadcasting like
logsumexp(randn(3, 4), rand(1, 4), dims=1). By not adding a size check it's more flexible.

Perhaps your concerns could be alleviated by thinking about combining some of David's work in #97 here. The only reason I latch onto this array method when my use case will be tuples is that the non-array method currently is pretty slow for my use.

[nalimilan]

@nalimilan I didn't want to over-specialize because then I couldn't use a call like logsumexp([1, 2], (0.5, 0.3)). This is useful for me because my weights are static in the calculations I need so I prefer to use a Tuple.

Could you use StaticArrays instead of tuples? I tend to prefer restrictive APIs unless we have a strong reason to be more flexible, as we can easily end up allowing things that weren't intended, and which can even return incorrect results.

As for the making the sizes the same, that would disallow broadcasting like
logsumexp(randn(3, 4), rand(1, 4), dims=1). By not adding a size check it's more flexible.

As I said the problem is that the definition would be inconsistent with StatsBase, which we should really avoid. With StatsBase conventions, logsumexp(randn(3, 4), rand(1, 4), dims=1) would be written as just logsumexp(randn(3, 4), rand(4), dims=1), but it's a bit more complicated to implement. So again unless there's a strong reason

Perhaps your concerns could be alleviated by thinking about combining some of David's work in #97 here. The only reason I latch onto this array method when my use case will be tuples is that the non-array method currently is pretty slow for my use.

That wouldn't really change anything AFAICT. With non-AbstractArray weights you can't check the size of the input to make sure it makes sense and that there's a one-to-one correspondence between values and weights.

[mileslucas]

Okay, well I don't really agree that this is better, but I've fully restricted the function to only work with abstract arrays which have exactly the same size. I've updated the tests and docstring, too.

[andreasnoack]

Should now be against https://github.com/JuliaStats/LogExpFunctions.jl since logsumexp has moved there.