# Random draws of multivariate normal with positive semi-definite covariance matrix

I am trying to generate random draws from a multivariate normal distribution with positive semi-definite covariance matrices. The Distributions.jl `MvNormal` only handles positive definite covariance matrices and throws an exception otherwise. There was a discussion about adding support for semi-definite covariance matrices a few years ago here but it is still not supported. While `MvNormal` does much more than generating random numbers, I am only interested in the random number generation.

My current solution is to do a pivoted cholesky decomposition and multiply by independent standard normals:

``````mu = [0,0]
Sig = [1 1;1 1]
chol = cholesky(Sig,Val(true),check=false)
random_draw = mu + chol.L[invperm(chol.p),1:chol.rank]*randn(chol.rank)
``````

I check that the eigenvalues of `Sig` are non-negative before I generate the draws: `minimum(eigvals(Sig)) >= 0`. This has worked properly in a few problems I have tested.

I would appreciate any feedback on my approach. Is there an easier and/or more efficient way to achieve the same results? Even better, is there any package that already generates multivariate normals with positive semi-definite covariance matrices?

The problem is that you may get slightly negative pivots due to roundoff errors, in which case LAPACK may terminate the factorization prematurely.

Ideally there would be a way to specify some tolerance so that slightly negative pivots could be treated as zeros, but I don’t know of any API for this?

You could use `sqrt(Hermitian(A))` instead of Cholesky. This is slower than Cholesky, but it should be more robust to small roundoff errors for semidefinite matrices because we added a tolerance check (which can be customized via an `rtol` keyword argument) to treat small negative eigenvalues as zero: RFC: treat small negative λ as 0 for sqrt(::Hermitian) by stevengj · Pull Request #35057 · JuliaLang/julia · GitHub

I am using the `check=false` option to ensure that the cholesky decomposition does not stop with small negative pivots. But I need to check that the matrix is positive semi-definite before I can use this, otherwise I may get wrong results without any exceptions if I pass a matrix that is not positive semi-definite. I am using the `eigvals(Sig)` to check for positive semi-definiteness.

The `sqrt` is doing the eigendecomposition already. Is there any optional argument to make it throw an exception if the eigenvalues are negative beyond the tolerance? By default it returns a complex matrix.

I’m not sure this does what you think. Read the LAPACK docs in the linked thread and LAPACK technical report #161 section 7: my understanding is that the Cholesky factorization will stop early if it encounters a negative pivot. `check=false` just means that you get the partially factorized result with no error thrown.

Wait, I just noticed that you are only looking at columns `1:chol.rank` of the `L` factor — I think you are right, it is okay in this case:

``````julia> B = randn(3,7); A = B'B; F = cholesky(A,Val(true),check=false);

julia> L = F.L[invperm(F.p),1:F.rank]; L*L' ≈ A
true
``````

First, if this doesn’t guarantee that the Cholesky factorization won’t encounter a slightly negative pivot due to roundoff errors. Second, you could easily have a matrix that is theoretically semidefinite but in practice has slightly negative eigenvalues, again due to roundoff errors. For example:

``````julia> B = randn(3,5); eigvals(Hermitian(B'B))
5-element Vector{Float64}:
-1.1122313274558508e-16
-1.6929564042211117e-17
0.39066211748341695
2.830678429996612
8.062953073438964
``````