Potential Performance regression between Julia 1.4.2 and 1.5 Beta 1

[lungben]

Hi,

I observed a strange effect comparing Timings for an Asian Pricer implementation in Julia 1.4.2 and 1.5 Beta 1.
For the vectorized implementation, Julia 1.5 is significantly faster than Julia 1.4.2 (which is great!).
However, for the loop implementation, it is the other way around - Julia 1.5 took nearly twice as long.

Could you please check? Or am I missing something obvious?

Thanks and Best Regards
Benjamin

using BenchmarkTools
using Random
using Statistics 

T, r, K, s, Xo = promote(1.0, 0.05, 55.0, 0.3, 50) # align types for performance
m, mc = 1000, 10_000

# Source: https://gist.github.com/bkamins/7ecb8bf493cac8d6f6f3c141eedbaf4b

function asian_sample_julia(T, r, K, σ, X₀, m::Integer)
    X = X₀
    x̂ = zero(X)
    Δ = T/m
    for i in 1:m
        X *= exp((r-σ^2/2)*Δ + σ*√Δ*randn())
        x̂ += X
    end
    exp(-r*T)*max(x̂/m - K, 0)
end

function asian_sample_vec_julia(T, r, K, σ, X₀, m::Integer; rng=randn)
    Δ = T/m
    X = rng(typeof(X₀), m)
    X .*= σ*√Δ
    X .+= (r-σ^2/2)*Δ
    X .= exp.(cumsum!(X, X))
    exp(-r*T)*max(mean(X)*X₀ - K, 0)
end

asian_loop_julia(T, r, K, σ, X₀, m::Integer, mc:: Integer, f; kwargs...) = mean(f(T, r, K, σ, X₀, m; kwargs...) for _ in 1:mc)

versioninfo()

println("Loop implementation:")
@btime asian_loop_julia($T, $r, $K, $s, $Xo, $m, $mc, $asian_sample_julia)

println("Vectorized implementation:")
@btime asian_loop_julia($T, $r, $K, $s, $Xo, $m, $mc, $asian_sample_vec_julia)

Output:

Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 8
Loop implementation:
  103.978 ms (0 allocations: 0 bytes)
Vectorized implementation:
  165.363 ms (10000 allocations: 77.51 MiB)
2.0182973558431843

Julia Version 1.5.0-beta1.0
Commit 6443f6c95a (2020-05-28 17:42 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 8
Loop implementation:
  191.175 ms (0 allocations: 0 bytes)
Vectorized implementation:
  115.469 ms (10000 allocations: 77.51 MiB)
2.0806637302421787

[KristofferC]

What happens if you take out the random number generation?

[lungben]

When I repace the randn() in the loop implementation by a deterministic function or by a uniform random number rand(), the performance is the same for Julia 1.4.2 and 1.5 Beta.

uniform random generator

Julia Version 1.4.2
Loop implementation:
  82.991 ms (0 allocations: 0 bytes)

Julia Version 1.5.0-beta1.0
Loop implementation:
  84.192 ms (0 allocations: 0 bytes)

The root cause seems to be the randn() function:
Julia 1.4.2:

julia> @btime randn()
  8.407 ns (0 allocations: 0 bytes)
-0.5964650134297568

julia> @btime rand()
  6.999 ns (0 allocations: 0 bytes)
0.7545456334265139

Julia 1.5 Beta:

julia> @btime randn()
  9.409 ns (0 allocations: 0 bytes)
0.5146215154773097

julia> @btime rand()
  6.600 ns (0 allocations: 0 bytes)
0.030946873310495837

The timings for rand() are essentially the same, whereas randn() takes significantly longer in 1.5.
However, when timing randn alone, it is only a 10% difference, whereas it is nearly a factor 2 for the whole function.

Interestingly, when generating randn numbers for a large vector, the picture changes completely:
Julia 1.4.2:

julia> @btime randn(1000);
  6.700 μs (1 allocation: 7.94 KiB)

Julia 1.5:

julia> @btime randn(1000);
  2.256 μs (1 allocation: 7.94 KiB)

[KristofferC]

Interestingly, when generating randn numbers for a large vector, the picture changes completely:

Yeah, performance of bulk generation of normally distributed random numbers was improved in #35078. Perhaps it had an effect on the scalar case?

However, when timing randn alone, it is only a 10% difference, whereas it is nearly a factor 2 for the whole function.

Might be some cache effects?

What happens if you create an RNG outside the loop and use that one?
Something like:

rng = MersenneTwister()
for ....
    x = ... randn(rng)
end

[rfourquet]

Yes for vectors the algorithm was specifically improved.

For scalars, I indeed oberve a 6% slowdown between 1.4 (9.127ns) and 1.5 (9.725). What is surprising to me is that by evaling the following (which is roughly the definition of randn in 1.4, the performance is back (8.932 ns), when it seems that this should be stricitly equivalent to the 1.5. version:

@eval Random @inline function randn(rng::AbstractRNG=default_rng())
          @inbounds begin
              r = rand(rng, UInt52Raw())
              r &= 0x000fffffffffffff
              rabs = Int64(r>>1) # One bit for the sign
              idx = rabs & 0xFF
              x = ifelse(r % Bool, -rabs, rabs)*wi[idx+1]
              rabs < ki[idx+1] && return x # 99.3% of the time we return here 1st try
              return randn_unlikely(rng, idx, rabs, x)
          end
      end

Maybe one @inline in 1.5 is discarded?

[lungben]

What happens if you create an RNG outside the loop and use that one?

When using randn(rng) with a pre-defined rng=MersenneTwister(), I get the following timings for the Asian pricing:

Julia 1.4.2:

Loop implementation:
  77.685 ms (4 allocations: 144 bytes)

julia> @btime randn($rng)
  5.999 ns (0 allocations: 0 bytes)
-0.12093968291229917

julia> @btime randn()
  8.908 ns (0 allocations: 0 bytes)
1.6566234022831083

julia> @btime randn(MersenneTwister())
  10.300 μs (9 allocations: 19.56 KiB)
-0.9996654920302606

Julia 1.5:

Loop implementation:
  193.106 ms (0 allocations: 0 bytes)

julia> @btime randn($rng)
  7.100 ns (0 allocations: 0 bytes)
0.571128035953164

julia> @btime randn()
  9.409 ns (0 allocations: 0 bytes)
-0.8816111967004022

julia> @btime randn(MersenneTwister())
  10.201 μs (8 allocations: 19.55 KiB)
-0.3469181045198116

Thus, the time differences for scalar randn numbers are even more exposed when pre-creating the rng.

Btw. great performance tip pre-creating the rng (I was not aware of this large time saving)!

[lungben]

I tried to compare the assembler generated by @code_native randn() for both versions.
I do not understand most of it, but there are a few differences (1st version is for Julia 1.4.2).

; │││┌ @ RNGs.jl:187 within `gen_rand'
; ││││┌ @ Base.jl:33 within `getproperty'
        movq    8(%rdi), %rcx
        movq    16(%rdi), %rax

; │││┌ @ RNGs.jl:187 within `gen_rand'
; ││││┌ @ Base.jl:33 within `getproperty'
        movq    8(%rdi), %r9
        testq   %r9, %r9
        je      L209
        leaq    48(%rsp), %rcx

This looks like a bounds-check (?)

Furthermore, for V1.4.2 there is a larger code-block (40 lines) at the end missing in the V1.5 assembly. I do not know what the last block is doing.

[lungben]

Replacing the function
@inline randn(rng::AbstractRNG=default_rng()) = _randn(rng, rand(rng, UInt52Raw()))
by the old implementation

@inline function randn(rng::AbstractRNG=default_rng())
          @inbounds begin
              r = rand(rng, UInt52Raw())
              r &= 0x000fffffffffffff
              rabs = Int64(r>>1) # One bit for the sign
              idx = rabs & 0xFF
              x = ifelse(r % Bool, -rabs, rabs)*wi[idx+1]
              rabs < ki[idx+1] && return x # 99.3% of the time we return here 1st try
              return randn_unlikely(rng, idx, rabs, x)
          end
      end

should bring the scalar performance back without loosing the performance optimization for larger arrays.

However, it is not a very nice solution: the code for _randn and randn is duplicated and it still would be good to understand why the new solution is slower for scalars in the first place.

[lungben]

Comparing the LLVM code between the new 1.5 implementation and the old one (by @eval-ling the old version above), I got the following difference:
In the Julia 1.5 version, there is the line
%30 = call double @j__randn_3665(%jl_value_t* nonnull %0, i64 %29) #0
which is not present in the old implementation.
Therefore, in the old implementation there are about 60 additional lines of code for an inlined function call.

Thus, I think @rfourquet is right, it is an inlining issue.
Unfortunately I do not know how to fix it...

[mbauman]

This is now uniformly faster for me on 1.7.

1.4

julia> @btime asian_loop_julia($T, $r, $K, $s, $Xo, $m, $mc, $asian_sample_julia)
  106.939 ms (0 allocations: 0 bytes)
1.9963023110637297

julia> println("Vectorized implementation:")
Vectorized implementation:

julia> @btime asian_loop_julia($T, $r, $K, $s, $Xo, $m, $mc, $asian_sample_vec_julia)
  122.707 ms (10000 allocations: 77.51 MiB)
2.00700254990639

julia> @btime randn(1000);
  5.596 μs (1 allocation: 7.94 KiB)

julia> @btime randn();
  6.967 ns (0 allocations: 0 bytes)

1.7-DEV

julia> @btime asian_loop_julia($T, $r, $K, $s, $Xo, $m, $mc, $asian_sample_julia)
  82.253 ms (0 allocations: 0 bytes)
1.9881008237998168

julia> println("Vectorized implementation:")
Vectorized implementation:

julia> @btime asian_loop_julia($T, $r, $K, $s, $Xo, $m, $mc, $asian_sample_vec_julia)
  102.999 ms (10000 allocations: 77.51 MiB)
2.0080692645661906

julia> @btime randn(1000);
  2.333 μs (1 allocation: 7.94 KiB)

julia> @btime randn();
  4.470 ns (0 allocations: 0 bytes)