# Fastest way to sample from MVN, changing parameters

**URL:** <https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145>\
**Category:** Julia at Scale\
**Tags:** statistics, distributions, gaussian-process, probablistic\
**Created:** [April 17, 2025, 2:15am UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145 "2025-04-17T02:15:12Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Daniel\_Johansson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_johansson/32/216620_2.png) [@Daniel\_Johansson](https://discourse.julialang.org/u/Daniel_Johansson)\
**Post date:** [April 17, 2025, 2:15am UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/1 "2025-04-17T02:15:12Z")

</div>

Hey guys!  
Does anybody know what is the most efficient way of sampling from a multivariate normal distribution whose parameters are constantly updated. I’m constructing a Gibbs sampler for a MVN whose parameters depend on matrix multiplication and inversion, thus I want to make things as efficient as possible.

---

<div class="post-metadata">

**Author:** ![Daniel\_Johansson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_johansson/32/216620_2.png) [@Daniel\_Johansson](https://discourse.julialang.org/u/Daniel_Johansson)\
**Post date:** [April 17, 2025, 2:18am UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/2 "2025-04-17T02:18:49Z")

</div>

For now, I’ve considered two approaches.

1. Simply, per iteration, declare a new multivariate normal distribution using the distributions package. i.e. declaring a new variable per iteration.
2. Computing the cholesky decomposition of my covariance matrix, and then sample from a standard normal of appropiate dimensions. Regardless, this seems to have much a higher cost than it is worht.

I would really value any insights.

---

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [April 17, 2025, 5:07am UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/3 "2025-04-17T05:07:19Z")

</div>

Hi and welcome to the Julia Discourse!

Wouldn’t 1. compute something like a Cholesky decomposition down the hood anyway? I’d guess 2. can be marginally quicker, but it should be easy to benchmark this (edit: no, see below).

If you really want to cut off a bit more, you can probably copy what Distributions.jl does down the hood and avoid allocating intermediate things, but I’m not sure if the speed-up would be worth the time.

---

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [April 17, 2025, 6:55am UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/5 "2025-04-17T06:55:00Z")

</div>

Looks like using the package is more efficient:

```julia
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

```julia
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.

```

---

<div class="post-metadata">

**Author:** ![Daniel\_Johansson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_johansson/32/216620_2.png) [@Daniel\_Johansson](https://discourse.julialang.org/u/Daniel_Johansson)\
**Post date:** [May 1, 2025, 10:06pm UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/6 "2025-05-01T22:06:33Z")

</div>

Hi thanks for the warm welcome!

It makes total sense. Any idea of why would calling the package is faster? I’m still a begginer in understanding why things are as fast as they are different programming languages, but I would to learn everything I can.

Thanks for taking the time btw!

---

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [May 2, 2025, 1:05pm UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/7 "2025-05-02T13:05:32Z")

</div>

I’ve quickly gone through the code. If I’m not missing anything, the following explains the difference: if you go to the source for Distributions.jl, you find that they use the function `unwhiten!` from PDMats.jl (see [Distributions.jl/src/multivariate/mvnormal.jl at efff906e2e6aad180d0be6dcaa9c98d3c398510d · JuliaStats/Distributions.jl · GitHub](https://github.com/JuliaStats/Distributions.jl/blob/efff906e2e6aad180d0be6dcaa9c98d3c398510d/src/multivariate/mvnormal.jl#L272C5-L272C13)) to transform a standard MVN to the desired variance-covariance structure. If you then go to PDMats.jl, I believe that it uses `unwhiten!` as written here ([PDMats.jl/src/pdmat.jl at 7b61f73aa52d53af980ac5b47a93f93bb9fe253e · JuliaStats/PDMats.jl · GitHub](https://github.com/JuliaStats/PDMats.jl/blob/7b61f73aa52d53af980ac5b47a93f93bb9fe253e/src/pdmat.jl#L134)):

```julia
function unwhiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
    @check_argdims axes(r) == axes(x)
    @check_argdims a.dim == size(x, 1)
    if r === x
        return lmul!(chol_lower(cholesky(a)), r)
    else
        return mul!(r, chol_lower(cholesky(a)), x)
    end
end

```

which suggests that it uses an in-place multiplication based on the cholesky-decomposition. But perhaps someone more versed in Distributions.jl can weigh in on this?

---

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [May 2, 2025, 1:55pm UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/8 "2025-05-02T13:55:37Z")

</div>

So, if you’re looking for an even faster version; simply copy the code from the `_rand!` method and make it use your PD matrix immediately, rather than reading it from a MVN object. You could even strip down `unwhiten!` a bit more for your use-case

---

<div class="post-metadata">

**Author:** ![Daniel\_Johansson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_johansson/32/216620_2.png) [@Daniel\_Johansson](https://discourse.julialang.org/u/Daniel_Johansson)\
**Post date:** [May 3, 2025, 4:11pm UTC](https://discourse.julialang.org/t/fastest-way-to-sample-from-mvn-changing-parameters/128145/9 "2025-05-03T16:11:23Z")

</div>

Nice! I didnt know things could get so deep for efficency.

Thanks for the help!