It seems this topic keeps coming up over and over again (see here and here), with a lot of efforts to resolve it but it gets stuck.
My code is reliant on using MvNormal for sampling. I could write my own code that checks if the covariance matrix I am using is positive definite or positive semi-definite. Then create my own special constructor for MvNormal in the degenerate case and my own rand method, but I feel very uneasy about overwriting a package.
Is there a cleaner way to do this?
Basically it means that your multivariate gaussian random vector is a linear transformation of a lower-dimensional standard gaussian random vector.
To find the right linear transformation, remove first the mean and then compute a (r,d) square root Matrix of your covariance. There are several algorithms to do that, and if the one used by MvNormal is not working in this particular case then pick another one by googling “square root of semi-positive deginite matrix in Julia”.
Then use X = mu + Gamma * MvNormal(r), where mu is the mean, gamma the matrix square root, and r the rank.