Is there a quicker way to obtain normalized multivariate random normal draws?

General Usage Performance
,
11

For example, let’s compare your method:

function doit(mu, V, n)
    X = rand(MvNormal(mu, Hermitian(V)), n)
    foreach(normalize!, eachcol(X))
    return X
end

to a version that expects static arrays, and returns an array of 3-component SVectors instead of a matrix:

using StaticArrays
function doit(mu::SVector, V::SMatrix, n)
    L = cholesky(V).U'
    X = [normalize(mu + L * @SVector(randn(3))) for _=1:n]
    return X
end

On my machine, I find that

mu_s = SA[2., 2., 2.]
V_s = SA[0.74 -0.08 0.34; -0.08 1.57 -0.5; 0.34 -0.5 1.16]
@btime doit($mu_s, $V_s, $n);

is almost twice as fast as @btime doit($mu, $V, $n); and if you pre-allocate X then you gain another 10%.

(And, as I said above, it may be even better if you merge this loop with any additional processing you want to do on these random vectors, rather than the traditional “vectorized” style of doing many individual passes over the data with built-in “vectorized” operations, which is only necessary in languages where loops are slow.)

Note that this implementation is essentially pulling out the internals of the MvNormal algorithm: if V has the Cholesky factorization V = L L^*, and r = randn(3), then x = \mu + L r has mean \mu and covariance V.

1 Like