We recently held a competition among the students in a Julia course we are teaching. The goal of the competition was to profile and optimize a provided particle-filter simulation. After all the obvious things were taken care of, most of the time was spent calculating log(sum(exp(w))) of a weight vec…

From what I can see, even though the llvm intrinsic @llvm.exp.v4f64 is output by SIMD.jl, in the end, LLVM still generates four separate calls: julia> @code_native exp(a) .text ... vzeroupper callq *%rdi <--------------- vmovapd %xmm0, %xmm7 vpermilpd …

Haha, so much for my manual SIMD efforts :stuck_out_tongue: Thanks for the explanation, I would not have figured that out. Do you have a feeling for where the issue lies? Is it SIMD.jl, Julia or LLVM that is to blame? @Elrod , thanks for you efforts on SIMDPirates, do you perhaps have any insight int…

LLVM. GCC will vectorize log/exp/sin/etc with the appropriate optimization flags, but LLVM needs those in addition to -fveclib=SVML or some other vector library. using SIMDPirates, SLEEFPirates, LoopVectorization function logsumexp_simdpirates!(w::Vector{T},we) where T offset = maximum(w) …

Thanks for this! On my machines, I still cannot beat Yeppp on Float64, but with SLEEFPirates and LoopVectorization I get significantly closer. julia> w = randn(Float64, N); we = similar(w); julia> @btime logsumexp!($w,$we); 10.358 μs (0 allocations: 0 bytes) julia> @btime logsumexp_yeppp!($w,$we)…

@generated function logsumexp_loopvectorization!(w::Vector{T},we) where T quote offset = maximum(w) N = length(w) s = zero(T) @vectorize $T for i = 1:N wl = w[i] wel = exp(wl-offset) we[i] = wel s += wel end …

I should comment that the offset algorithm that you are using for logsumexp is potentially problematic — you can easily contrive a case where it is off by a factor of 2. In particular, try [1e-20, log(1e-20)] with your logsumexp algorithm function f(x) X = maximum(x) return X + log(sum(exp…

Thanks for the link! I had not considered it so carefully before, but will sure make use of this trick in the future.

It looks like you’re mostly using logsumexp! for sampling from a categorical, is that right? Have you tried using the Gumbel max trick for this? It’s pretty great for sampling a categorical with unnormalized log-probability weights.

That trick does not return the log likelihood like the calculation above, right? Also, they mention that all the log(log(U)) can be precalculated, doesn’t this cause some bias?

Fast logsumexp

General Usage Performance

baggepinnen April 22, 2019, 1:19pm 9

For future reference, here is my implementation
https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/4f9b688d298157dc24a5b0a518d971221fbe15dd/src/resampling.jl#L10

Topic		Replies	Views
Fast LogSumExp over 4th dimension Performance tullio , loopvectorization	36	2593	October 21, 2022
Is this function well optimized for speed? General Usage	13	1026	November 13, 2019
Speeding up my logsumexp function Performance performance , array , tullio , loopvectorization	35	5426	October 21, 2022
Comparing exp() performance on Julia versus numpy Performance	38	2278	August 8, 2022
Reduce(logaddexp,x) has arbitrarily better accuracy than logsumexp(x) Numerics	11	800	January 2, 2024

Fast logsumexp

Related topics