# How to sample from MvNormal without allocating?

**URL:** <https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646>\
**Category:** Statistics\
**Tags:** question\
**Created:** [March 27, 2023, 3:40am UTC](https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646 "2023-03-27T03:40:43Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![markmbaum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markmbaum/32/32745_2.png) [@markmbaum](https://discourse.julialang.org/u/markmbaum)\
**Post date:** [March 27, 2023, 3:40am UTC](https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646/1 "2023-03-27T03:40:44Z")

</div>

Even using the in-place method `rand!` to sample from a `MvNormal` distribution seems to allocate some memory and be orders of magnitude slower than the univariate case. I only have a two dimensional case, so I was hoping to find a way to optimize but can’t seem to find out how.

Suggestions?

---

<div class="post-metadata">

**Author:** ![oxinabox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oxinabox/32/206603_2.png) [@oxinabox](https://discourse.julialang.org/u/oxinabox)\
**Post date:** [March 27, 2023, 4:13am UTC](https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646/2 "2023-03-27T04:13:44Z")

</div>

There is no obvious reason to me that this should be allocating.  
Here is the code for reference if anyone wants to go digging

> <https://github.com/JuliaStats/Distributions.jl/blob/ec68da3a8d4a4776367f2d7ca5ec2d4666e29c78/src/multivariate/mvnormal.jl#L276-L290>

> <https://github.com/JuliaStats/PDMats.jl/blob/fff131e11e23403931a42f5bfb3384f0d2b114c9/src/generics.jl#L40C10-L43>

`cholesky` on a `PDMat` should be nonallocating as it does that upfront during construction

---

<div class="post-metadata">

**Author:** ![markmbaum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markmbaum/32/32745_2.png) [@markmbaum](https://discourse.julialang.org/u/markmbaum)\
**Post date:** [March 27, 2023, 4:25am UTC](https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646/3 "2023-03-27T04:25:08Z")

</div>

A simple example:

```julia
using BenchmarkTools, Distributions, Random
using Random: rand!

X = MvNormal([2 1; 1 3])
y = zeros(2)
rng = Xoshiro()
@btime rand!($rng, $X, $y);

```

which produces

```julia
  133.021 ns (2 allocations: 96 bytes)

```

It seems to matter that the distribution is a `ZeroMeanFullNormal` but I’m not sure why.

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [April 3, 2023, 3:10pm UTC](https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646/4 "2023-04-03T15:10:51Z")

</div>

> [@markmbaum](#):
>
> ```julia
> using BenchmarkTools, Distributions, Random
> using Random: rand!
> 
> X = MvNormal([2 1; 1 3])
> y = zeros(2)
> rng = Xoshiro()
> @btime rand!($rng, $X, $y);
> 
> ```

There seem to be two sources of allocations. The first is [`PDMats.chol_lower`](https://github.com/JuliaStats/PDMats.jl/blob/fff131e11e23403931a42f5bfb3384f0d2b114c9/src/chol.jl#L6-L8) called [here](https://github.com/JuliaStats/PDMats.jl/blob/fff131e11e23403931a42f5bfb3384f0d2b114c9/src/generics.jl#L42). Not certain why this allocates; maybe because the result is a typeunion? This allocation is 16 bytes.

The second is [this line](https://github.com/JuliaStats/Distributions.jl/blob/ec68da3a8d4a4776367f2d7ca5ec2d4666e29c78/src/multivariate/mvnormal.jl#L278), where the mean `μ` is added to the result via broadcast. The mean in this case is a `FillArrays.Zeros`, and it seems that `broadcasted` makes a copy:[https://github.com/JuliaArrays/FillArrays.jl/blob/c3b38add861d475aadc66a112e045d7e0db31372/src/fillbroadcast.jl#L205](https://github.com/JuliaArrays/FillArrays.jl/blob/c3b38add861d475aadc66a112e045d7e0db31372/src/fillbroadcast.jl#L205). This results in an 80 byte allocation. Seems like this could be improved in FillArrays.

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [April 11, 2023, 12:01pm UTC](https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646/5 "2023-04-11T12:01:07Z")

</div>

> [@sethaxen](#):
>
> The second is [this line](https://github.com/JuliaStats/Distributions.jl/blob/ec68da3a8d4a4776367f2d7ca5ec2d4666e29c78/src/multivariate/mvnormal.jl#L278), where the mean `μ` is added to the result via broadcast. The mean in this case is a `FillArrays.Zeros`, and it seems that `broadcasted` makes a copy:[FillArrays.jl/fillbroadcast.jl at c3b38add861d475aadc66a112e045d7e0db31372 · JuliaArrays/FillArrays.jl · GitHub](https://github.com/JuliaArrays/FillArrays.jl/blob/c3b38add861d475aadc66a112e045d7e0db31372/src/fillbroadcast.jl#L205). This results in an 80 byte allocation. Seems like this could be improved in FillArrays.

Seems when I posted this, @jishnub had already begun working on a PR to fix this: [don't materialize when broadcasting Zeros with Vector by jishnub · Pull Request #211 · JuliaArrays/FillArrays.jl · GitHub](https://github.com/JuliaArrays/FillArrays.jl/pull/211)

---

<div class="post-metadata">

**Author:** ![markmbaum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markmbaum/32/32745_2.png) [@markmbaum](https://discourse.julialang.org/u/markmbaum)\
**Post date:** [August 15, 2023, 2:29pm UTC](https://discourse.julialang.org/t/how-to-sample-from-mvnormal-without-allocating/96646/6 "2023-08-15T14:29:28Z")

</div>

Seems like this has improved. Running the same tiny example

```julia
using BenchmarkTools, Distributions
using Random: Xoshiro, rand!

X = MvNormal([2 1; 1 3])
y = zeros(2)
rng = Xoshiro()
@btime rand!($rng, $X, $y);

```

69.586 ns (1 allocation: 16 bytes)