Sampling from matrix-variate normal distribution with singular covariances? [duplicate]

Question

The matrix-variate normal distribution can be sampled indirectly by utilizing the Cholesky decomposition of two positive definite covariance matrices. However, if one or both of the covariance matrices are positive semi-definite and not positive definite (for example a block structure due to several pairs of perfectly correlated features and samples) the Cholesky decomposition fails, e.g.

$\Sigma_{A} = \begin{bmatrix} 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 1 & 1 \end{bmatrix} \quad $or another example:$ \quad \Sigma_{B} = \begin{bmatrix} 4 & 14 & 0 & 0\\ 14 & 49 & 0 & 0\\ 0 & 0 & 25 & 20\\ 0 & 0 & 20 & 16 \end{bmatrix}$

Where $\Sigma_{B}$ is generated from $R$ (correlation matrix this time) = $\Sigma_{A} $

and $D$ (standard deviations) = $\begin{bmatrix} 2 & 0 & 0 & 0\\ 0 & 7 & 0 & 0\\ 0 & 0 & 5 & 0\\ 0 & 0 & 0 & 4 \end{bmatrix}$ via $RDR$.

Is it possible to adapt the SVD based sampling technique for the multivariate normal case that overcomes this difficulty to the matrix-variate case?

---

This question is different from this post in that it is not clear if the lower diagonal produced by the SVD based sampling technique will suffice, since it is potentially quite different from one produced by a Cholesky decomposition that might be performed in this case by removing duplicate features and/or samples from the covariance matrices, performing the decomposition, and putting them back in. Also, the mentioned post is not concerned with positive semi-definite matrices.

Answer 1

This sounds more like an issue with singular covariance matrices than with random matrices vs. random vectors. To handle the latter issue, do everything as a random vector, and then in the last step, reshape the vector into a matrix.

To handle the former problem:

If your desired covariance matrix is singular...

Let $\Sigma$ be a singular covariance matrix. Because it's singular, you can't do a Cholesky decomposition. But you can do a singular value decomposition.

[U, S, V] = svd(Sigma)

The singular value decomposition will construct matrices $U$, $S$, and $V$ such that $ \Sigma = U S V'$ and $S$ is diagonal. Furthermore, $U = V$ (because $\Sigma$ is symmetric). The number of positive singular values will be the rank of your covariance matrix.

You can then construct $n$ random vectors of length $k$ with.

X = randn(n, k) * sqrt(S) * U'

Let $\mathbf{z}$ be a standard multivariate normal vector. The basic idea is that:

\begin{align*} \mathrm{Var}\left(US^{\frac{1}{2}} \mathbf{z} \right) &= US^{\frac{1}{2}} \mathrm{Var}\left( \mathbf{x} \right)S^{\frac{1}{2}}U' \\ &= U S U'\\ &= \Sigma \end{align*}

Once you get a vector $US^{\frac{1}{2}}\mathbf{z}$ you can simply reshape it to the dimensions of your matrix. (Eg. a 6 by 1 vector could become a 3 by 2 matrix.)

The SVD on a symmetric matrix $C$ is a way to find another matrix $A$ such that $AA' = C$.

Optional step to be a more clever mathematician and more efficient coder

Find the singular values below some tolerance and remove them and their corresponding columns from $S$ and $U$ respectively. This way you can generate less than a $k$ dimensional random vector $\mathbf{z}$.

Answer 2

If the only problem is in failure of Cholesky decomposition, then I'd try reducing the rank of matrix then duplicating the removed variables. Here's MATLAB example:

CODE:

%% singular cov matrix
C = [1 0.5 0.5; 0.5 1 1; 0.5 1 1]
chol(C)

%% reduce rank by removing duplicates
c = C(1:2,1:2)
h = chol(c)

%% bring back the twins

H = [h h(:,2)]

%% test Cholesky
r = randn(500,2);
cr = r*H;

cov(cr)

OUTPUT

C =

    1.0000    0.5000    0.5000
    0.5000    1.0000    1.0000
    0.5000    1.0000    1.0000

Error using chol
Matrix must be positive definite.


c =

    1.0000    0.5000
    0.5000    1.0000


h =

    1.0000    0.5000
         0    0.8660


H =

    1.0000    0.5000    0.5000
         0    0.8660    0.8660


ans =

    1.0063    0.4705    0.4705
    0.4705    0.9914    0.9914
    0.4705    0.9914    0.9914

>> 

For the matrix in your example, H matrix will be very simple:

1 1 0 0
0 0 1 1