What about this? You could still de-vectorize the reduction functions (sum
and maximum
), but it seems to perform pretty well:
import StatsFuns
using BenchmarkTools
function logsumexp1(u, axes)
m = maximum(u, axes)
return m .+ log.(sum(exp.(u .- m), axes))
end
logsumexp2(u, axes) = mapslices(StatsFuns.logsumexp, u, axes) # for comparison
julia> u = rand(10,100000);
julia> @btime logsumexp1($u, 2);
22.964 ms (15 allocations: 7.63 MiB)
julia> @btime logsumexp2($u, 2)
33.472 ms (138 allocations: 786.13 KiB)
julia> u = rand(100000,10);
julia> @btime logsumexp1($u, 2);
25.086 ms (18 allocations: 9.92 MiB)
julia> @btime logsumexp2($u, 2);
223.134 ms (1099527 allocations: 25.17 MiB)
You don’t need to add axes with newaxis
since Julia’s reductions preserve the number of dimensions.