Apply reduction along specific axes

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.