Generating samples from singular Gaussian distribution

Question

Let random vector $x = (x_1,...,x_n)$ follow multivariate normal distribution with mean $m$ and covariance matrix $S$. If $S$ is symmetric and positive definite (which is the usual case) then one can generate random samples from $x$ by first sampling indepently $r_1,...,r_n$ from standard normal and then using formula $m + Lr$, where $L$ is the Cholesky lower factor so that $S=LL^T$ and $r = (r_1,...,r_n)^T$.

What about if one wants samples from singular Gaussian i.e. $S$ is still symmetric but not more positive definite (only positive semi-definite). We can assume also that the variances (diagonal elements of $S$) are strictly positive. Then some elements of $x$ must have linear relationship and the distribution actually lies in lower dimensional space with dimension $<n$, right?

It is obvious that if e.g. $n=2, m = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, S = \begin{bmatrix} 1 & 1 \\ 1 & 1\end{bmatrix}$ then one can generate $x_1 \sim N(0,1)$ and set $x_2=x_1$ since they are fully correlated. However, is there any good methods for generating samples for general case $n>2$? I guess one needs first to be able identify the lower dimensional subspace, then move to that space where one will have valid covariance matrix, then sample from it and finally deduce the values for the linearly dependent variables from this lower-dimensional sample. But what is the best way to that in practice? Can someone point me to books or articles that deal with the topic; I could not find one.

Answer 1

The singular Gaussian distribution is the push-forward of a nonsingular distribution in a lower-dimensional space. Geometrically, you can take a standard Normal distribution, rescale it, rotate it, and embed it isometrically into an affine subspace of a higher dimensional space. Algebraically, this is done by means of a Singular Value Decomposition (SVD) or its equivalent.

---

Let $\Sigma$ be the covariance matrix and $\mu$ the mean in $\mathbb{R}^n$. Because $\Sigma$ is non-negative definite and symmetric, the SVD will take the form

$$\Sigma = U \Lambda^2 U^\prime$$

for an orthogonal matrix $U\in O(n)$ and a diagonal matrix $\Lambda$. $\Lambda$ will have $m$ nonzero entries, $0\le m \le n$.

Let $X$ have a standard Normal distribution in $\mathbb{R}^m$: that is, each of its $m$ components is a standard Normal distribution with zero mean and unit variance. Abusing notation a little, extend the components of $X$ with $n-m$ zeros to make it an $n$-vector. Then $U\Lambda X$ is in $\mathbb{R}^n$ and we may compute

$$\text{Cov}(U\Lambda X) = U \Lambda\text{Cov}(X) \Lambda^\prime U^\prime = U \Lambda^2 U^\prime = \Sigma.$$

Consequently

$$Y = \mu + U\Lambda X$$

has the intended Gaussian distribution in $\mathbb{R}^n$.

It is of interest that this works when $n=m$: that is to say, this is a (standard) method to generate multivariate Normal vectors, in any dimension, for any given mean $\mu$ and covariance $\Sigma$ by using a univariate generator of standard Normal values.

---

As an example, here are two views of a thousand simulated points for which $n=3$ and $m=2$:

The second view, from edge-on, demonstrates the singularity of the distribution. The R code that produced these figures follows the preceding mathematical exposition.

#
# Specify a Normal distribution.
#
mu <- c(5, 5, 5)
Sigma <- matrix(c(1, 2, 1,
                  2, 3, 1,
                  1, 1, 0), 3)
#
# Analyze the covariance.
#
n <- dim(Sigma)[1]
s <- svd((Sigma + t(Sigma))/2) # Guarantee symmetry
s$d <- abs(zapsmall(s$d))
m <- sum(s$d > 0)
#$
# Generate a standard Normal `x` in R^m.
#
n.sample <- 1e3 # Number of points to generate
x <- matrix(rnorm(m*n.sample), nrow=m)
#
# Embed `x` in R^n and apply the square root of Sigma obtained from its SVD.
#
x <- rbind(x, matrix(0, nrow=n-m, ncol=n.sample))
y <- s$u %*% diag(sqrt(s$d)) %*% x + mu
#
# Plot the results (presuming n==3).
#
library(rgl)
plot3d(t(y), type="s", size=1, aspect=TRUE, 
       xlab="Y1", ylab="Y2", zlab="Y3", box=FALSE,
       col="Orange")