# How to sample a singular matrix-variate normal distribution?

What is known

The wikipedia article for the matrix-variate normal distribution has a section on how to efficiently sample from it without using the Kronecker product (which is computationally prohibitive in this case). To summarize, given the Cholesky decomposition of two positive definite covariance matrices $C = AA^{T}$ and $D = B^{T}B$ and a random sample from the standard normal $X$ (with the corresponding dimensions), a sample from the matrix-variate normal distribution is efficiently obtained via $Y = AXB$.

The problem

In the case of singular covariance matrices, due to linearly related samples or features in $Y$, it is not possible to perform the Cholesky decomposition. In a multivariate normal setting one would simply resolve this by:

taking a standard normal distribution, rescaling and rotating it, and finally embedding it isometrically into an affine subspace of a higher dimensional space. Algebraically, this is simply done by means of a singular value decomposition (see here).

However, in the case of a matrix-variate normal it is not clear how to proceed.

Simply exchanging the Cholesky decompositions for an SVD style square root operation does not seem to yield correct results (see here) $-$ possibly because of the difference between $AA^{T}$ and $B^{T}B$ in requiring the lower and upper triangular matrix of the Cholesky decomposition respectively and therefore some inherent difference to the SVD approach.

Potential solution

It's not possible to have $Y$ and $Y^{T}$ yield the standard deviations that were initially specified and from which (together with the respective correlation matrices) the covariance matrices were initially generated. Since, the expected covariances of $Y$ are different from $C$ and $D$ due to the additional trace($C$) and trace($D$) factors (see here).

Page 69 of this book may be helpful and this post.