Fastest way to sample from MVN, changing parameters

Specific Domains Julia at Scale
, , ,
5

Looks like using the package is more efficient:

using Distributions, BenchmarkTools, Random, LinearAlgebra
Random.seed!(42)

function with_intermediate(μΣ, N, M, d)
    sum_out = 0.
    for i in 1:M
        μ = rand(d)
        Σ = rand(μΣ)
        X = rand(MvNormal(μ, Σ), N)
        sum_out += X[1,1]
    end
    return sum_out
end

function without_intermediate(μΣ, N, M, d)
    sum_out = 0.
    for i in 1:M
        μ = rand(d)
        S = cholesky(rand(μΣ)).L
        X = muladd(S, randn(d,N), μ)
        sum_out += X[1,1]
    end
    return sum_out
end


function test(d, N, M)
    S = cholesky(I(d))
    μΣ = Wishart(d, S)
    with_intermediate(μΣ, N, M, d)
    without_intermediate(μΣ, N, M, d)
    display(@benchmark with_intermediate($μΣ, $N, $M, $d))
    display(@benchmark without_intermediate($μΣ, $N, $M, $d))
end
test(4, 100, 10)

gives

BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  14.792 μs …  53.984 ms  ┊ GC (min … max):  0.00% … 99.91%
 Time  (median):     17.541 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   23.199 μs ± 539.668 μs  ┊ GC (mean ± σ):  23.25% ±  1.00%

               ▃▇▅█▆▄▁                                          
  ▁▁▁▂▂▃▄▃▃▃▃▄▆████████▆▅▄▃▃▂▂▂▁▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  14.8 μs         Histogram: frequency by time         24.6 μs <

 Memory estimate: 40.47 KiB, allocs estimate: 60.
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  17.791 μs …  57.846 ms  ┊ GC (min … max):  0.00% … 99.90%
 Time  (median):     24.458 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   33.310 μs ± 586.611 μs  ┊ GC (mean ± σ):  26.21% ±  3.13%

                 ▂▄▆▅▇██▇▆▂                                     
  ▁▂▂▂▂▅▆▅▄▃▃▂▃▅▇███████████▇▅▄▃▃▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  17.8 μs         Histogram: frequency by time         37.6 μs <

 Memory estimate: 109.84 KiB, allocs estimate: 90.
1 Like