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

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.

marked as duplicate by whuber♦Oct 6 '16 at 13:24

• When you say, "this method breaks down." Is the basic problem you can't do a Cholesky decomposition of a singular matrix? – Matthew Gunn Oct 5 '16 at 19:13
• Why not remove duplicate variables before Cholesky (or similar matrix square root)? Those with 100% correlation you can simply copy from their twins afterwards – Aksakal Oct 5 '16 at 19:53
• @MatthewGunn I guess that seems to be it. Then the solution is basically similar to what was posted for the multivariate Gaussian. – baf84b4c Oct 6 '16 at 8:50
• @Aksakal How would you do this specifically in a simple example? – baf84b4c Oct 6 '16 at 13:09
• A comment I wrote in the duplicate thread points out that matrix square roots are not unique. Therefore, your question is not any different from the question there: it is immaterial that the SVD might produce a different square root than the Cholesky decomposition. Moreover, your question explicitly concerns semidefinite matrices: they are positive but of reduced rank, exactly as in the duplicate. The duplicate completely and thoroughly answers your questions. Moreover, some of the answers you have already received here duplicate answers there! – whuber Oct 6 '16 at 15:13

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}$.

• Actually, the lower triangular obtained in this way ($US^{1/2}$) is different from the one produced by the Cholesky decomposition. See: math.stackexchange.com/a/307208. Not sure if this will do for sampling from the matrix-variate normal in the way given here: en.wikipedia.org/wiki/…. – baf84b4c Oct 6 '16 at 12:53
• @baf84b4c The SVD doesn't produce a lower triangular matrix. It is not a Cholesky decomposition. – Matthew Gunn Oct 6 '16 at 15:45
• Ah, yes. So, theoretically, if I had a positive definite matrix I could do both Cholesky decomposition and the SVD approach (which would be overkill). This means I would get a different kind of sampling by doing $LX$, where $L$ is the lower triangular obtained from the Cholesky decomposition and $US^{1/2}X$ according to your approach. And both should be correct? – baf84b4c Oct 6 '16 at 19:55
• Let $\mathbf{z}$ be distributed standard multivariate normal. $\mathrm{Var}(A'\mathbf{Z}) =A'\mathrm{Var}(\mathbf{Z}) A= A'A$. So ANY matrix $A$ such that $A'A=\Sigma$ will make $\mathbf{x} = A\mathbf{z}$ have covariance matrix $\Sigma$ that you want. And $A$ is not unique!. Perhaps do out the two-dimensional case by hand if you want to see that. If $x_1 = a z_1 + bz_2$ and $x_2 = c z_1 + d z_2$ you have four variables: $a$, $b$, $c$ and $d$. But there are only three constraints (i.e. (i) variance of $x_1$ (ii) variance of $x_2$ and (iii) covariance of $x_1$ and $x_2$) – Matthew Gunn Oct 6 '16 at 20:06

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
• For the mentioned case of binary block structured covariance matrices, there seems to be an even simpler solution that doesn't involve any Cholesky decomposition, since it would just be a diagonal matrix after removing the duplicate pairs? – baf84b4c Oct 6 '16 at 14:06
• If it's all 1's and 0's, then yes, it's even simpler. You don't need any Cholesky at all. Just get the IIDs and duplicate as necessary. – Aksakal Oct 6 '16 at 14:36
• Thanks. I see now that I can't simply reduce the rank by removing duplicates because the duplicates are actually in the correlation matrix and not the covariance matrix. For the later case an SVD based approach might be necessary. – baf84b4c Oct 6 '16 at 15:00