Simulating Matrix-Normal Observations

Why the Simulation Construction Matters

A matrix-normal model allows entries to be dependent in two structured ways: observations in different rows can share row covariance, and variables in different columns can share column covariance. To simulate from the model, one needs to impose both dependence patterns without changing the desired mean. Cholesky factors provide a direct construction from independent standard normal variables.

The key is to keep the multiplication order consistent with the chosen convention for row and column covariance. In this tutorial, the first covariance parameter describes rows and the second describes columns. The construction below makes that convention explicit and gives a way to check it algebraically.

Definition: Let \(M\) be an \(r\times c\) matrix, and let \(U\) and \(V\) be positive-definite matrices of dimensions \(r\times r\) and \(c\times c\), respectively. A random \(r\times c\) matrix \(X\) has a matrix-normal distribution with mean \(M\), row covariance \(U\), and column covariance \(V\), written \(X\sim MN_{r,c}(M,U,V)\), when \(\operatorname{vec}(X)\) is multivariate normal with mean \(\operatorname{vec}(M)\) and covariance \(V\otimes U\). Here \(\operatorname{vec}\) stacks the columns of a matrix.

Equivalently, the entrywise covariance convention is

E[X_{ij}]=M_{ij}, \qquad \operatorname{Cov}(X_{ij},X_{k\ell})=U_{ik}V_{j\ell}.

Thus \(U\) relates entries in different rows, while \(V\) relates entries in different columns. In particular, \(\operatorname{Var}(X_{ij})=U_{ii}V_{jj}\). These marginal variances are products of the corresponding diagonal entries; \(U\) and \(V\) are not, in general, the covariance matrices of individual rows or columns without a scale factor.

The Cholesky Construction

For positive-definite \(U\) and \(V\), let \(A\) and \(B\) be their lower-triangular Cholesky factors:

U=AA^{\mathsf T}, \qquad V=BB^{\mathsf T}.

Generate an \(r\times c\) matrix \(Z\) whose entries are mutually independent \(N(0,1)\) random variables. Then set \(X=M+AZB^{\mathsf T}\). The sizes are \(A:r\times r\), \(Z:r\times c\), and \(B^{\mathsf T}:c\times c\), so the resulting matrix has the required \(r\times c\) shape.

Theorem: If \(Z\) has mutually independent standard normal entries and \(U=AA^{\mathsf T}\), \(V=BB^{\mathsf T}\), then \(X=M+AZB^{\mathsf T}\) has distribution \(MN_{r,c}(M,U,V)\).

Proof. Column-wise vectorization satisfies \(\operatorname{vec}(AZB^{\mathsf T})=(B\otimes A)\operatorname{vec}(Z)\). The entries of \(Z\) are independent standard normal variables, so \(\operatorname{vec}(Z)\) is multivariate normal with mean zero and covariance \(I_{rc}\). A linear transformation of a multivariate normal vector is multivariate normal. Therefore \(\operatorname{vec}(X)\) is normal with mean \(\operatorname{vec}(M)\), and its covariance is

(B\otimes A)I_{rc}(B\otimes A)^{\mathsf T} =(BB^{\mathsf T})\otimes(AA^{\mathsf T}) =V\otimes U.

This is the defining mean and covariance of \(MN_{r,c}(M,U,V)\), which proves the result. \(\square\)

Covariance check: For any two indices \(i,j\) and \(k,\ell\), the construction gives \(\operatorname{Cov}(X_{ij},X_{k\ell})=U_{ik}V_{j\ell}\). This entrywise identity provides a direct check of which factor controls rows and which controls columns.

Worked Simulations

Worked Example: A Two-by-Two Draw With Correlated Rows and Columns

Take

U= \begin{pmatrix}4&1\\1&1\end{pmatrix}, \qquad V= \begin{pmatrix}1&1/2\\1/2&1\end{pmatrix}, \qquad M= \begin{pmatrix}10&20\\-1&3\end{pmatrix}.

Both covariance matrices are positive definite. For example, their leading principal minors are \(4,3\) for \(U\), and \(1,3/4\) for \(V\), all positive. Their Cholesky factors are

A= \begin{pmatrix}2&0\\1/2&\sqrt{3}/2\end{pmatrix}, \qquad B= \begin{pmatrix}1&0\\1/2&\sqrt{3}/2\end{pmatrix}.

Multiplication verifies \(AA^{\mathsf T}=U\) and \(BB^{\mathsf T}=V\). For one simulated draw, suppose the independent standard normal matrix is

Z=\begin{pmatrix}1&-1\\2&0\end{pmatrix}.

First multiply on the left:

AZ= \begin{pmatrix}2&-2\\1/2+\sqrt{3}&-1/2\end{pmatrix}.

Since \(B^{\mathsf T}=\begin{pmatrix}1&1/2\\0&\sqrt{3}/2\end{pmatrix}\), multiplying on the right gives

AZB^{\mathsf T} = \begin{pmatrix}2&1-\sqrt{3}\\1/2+\sqrt{3}&1/4+\sqrt{3}/4\end{pmatrix}.

For instance, the lower-right entry is \((1/2+\sqrt{3})(1/2)+(-1/2)(\sqrt{3}/2)=1/4+\sqrt{3}/4\). Adding \(M\), the simulated observation is

X= \begin{pmatrix} 12&21-\sqrt{3}\\ -1/2+\sqrt{3}&13/4+\sqrt{3}/4 \end{pmatrix} \approx \begin{pmatrix}12&19.268\\1.232&3.683\end{pmatrix}.

The simulated values are one random draw, not covariance estimates. The target covariance can nevertheless be checked exactly from \(U\) and \(V\). For example, \(\operatorname{Var}(X_{11})=U_{11}V_{11}=4\), \(\operatorname{Cov}(X_{11},X_{21})=U_{12}V_{11}=1\), and \(\operatorname{Cov}(X_{11},X_{12})=U_{11}V_{12}=2\). In column-stacked order \((X_{11},X_{21},X_{12},X_{22})\), the full covariance is

V\otimes U= \begin{pmatrix} 4&1&2&1/2\\ 1&1&1/2&1/2\\ 2&1/2&4&1\\ 1/2&1/2&1&1 \end{pmatrix}.

This calculation checks both covariance directions: changing row indices uses entries of \(U\), while changing column indices uses entries of \(V\).

Worked Example: Diagonal Covariances and Entrywise Scaling

Let \(r=2\), \(c=3\), \(M=0\), \(U=\operatorname{diag}(1,4)\), and \(V=\operatorname{diag}(1,4,9)\). The Cholesky factors are \(A=\operatorname{diag}(1,2)\) and \(B=\operatorname{diag}(1,2,3)\). Suppose the generated standard normal matrix is

Z=\begin{pmatrix}1&-1&2\\1/2&0&-1\end{pmatrix}.

Because both factors are diagonal, left multiplication scales the second row by \(2\), while right multiplication scales the second and third columns by \(2\) and \(3\). Thus

AZ= \begin{pmatrix}1&-1&2\\1&0&-2\end{pmatrix}, \qquad X=AZB^{\mathsf T} =\begin{pmatrix}1&-2&6\\1&0&-6\end{pmatrix}.

The target variance of the third entry in the first row is \(U_{11}V_{33}=1\cdot9=9\). The covariance between the first and second entries in that row is \(U_{11}V_{12}=1\cdot0=0\). The covariance between \(X_{11}\) and \(X_{21}\) is \(U_{12}V_{11}=0\). In fact, all distinct entries have zero covariance here. Because the entries are jointly normal, zero covariance between distinct entries means they are independent; their variances are \(1,4,9\) in row one and \(4,16,36\) in row two. The displayed numerical matrix is one draw, while these variances and independence statements describe the distribution of repeated draws.

Worked Example: Checking the Row and Column Scaling for a Fixed Entry

Suppose \(U=\begin{pmatrix}2&1\\1&3\end{pmatrix}\) and \(V=\begin{pmatrix}5&-1\\-1&2\end{pmatrix}\). These are positive definite: their determinants are \(5\) and \(9\), respectively, and their diagonal entries are positive. Consider any \(2\times2\) simulated matrix \(X=M+AZB^{\mathsf T}\) constructed from their Cholesky factors.

The variance of \(X_{12}\) is \(U_{11}V_{22}=2\cdot2=4\). Its covariance with \(X_{22}\), which is in the same column but a different row, is \(U_{12}V_{22}=1\cdot2=2\). Its covariance with \(X_{11}\), which is in the same row but a different column, is \(U_{11}V_{21}=2(-1)=-2\). The negative value is valid: the column covariance \(V_{21}\) is negative, so entries in different columns can be negatively associated. The covariance with \(X_{21}\), which differs in both indices, is \(U_{12}V_{21}=1(-1)=-1\). These four calculations verify that the same construction handles row dependence, column dependence, and their combination.

A Practical Simulation Procedure

For each draw, use the following sequence. Reuse the Cholesky factors across draws when \(M,U,V\) are fixed; generate a fresh independent standard normal matrix each time.

Check dimensions and definiteness.
Confirm that \(U\) is \(r\times r\), \(V\) is \(c\times c\), and both are positive definite so their Cholesky factors exist.

Factor both covariance matrices.
Compute lower-triangular factors \(A\) and \(B\) satisfying \(U=AA^{\mathsf T}\) and \(V=BB^{\mathsf T}\).

Generate independent noise.
Draw \(rc\) independent \(N(0,1)\) values and arrange them as an \(r\times c\) matrix \(Z\).

Transform and shift.
Calculate \(X=M+AZB^{\mathsf T}\), preserving this multiplication order.

Why the Covariance Check Matters

In a simulation study, the construction determines whether the simulated data match the assumed model. When many independent matrices are generated, empirical means and covariances can be compared with their targets. For example, the empirical covariance between the \((i,j)\) and \((k,\ell)\) entries should approach \(U_{ik}V_{j\ell}\) as the number of simulated matrices grows, under the usual law of large numbers conditions. A finite simulation will fluctuate around the target; agreement is expected in the limit, not exactly in every finite run.

The Kronecker covariance \(V\otimes U\) is also a useful implementation check. Compute it in column-stacked order and confirm that its entries match \(U_{ik}V_{j\ell}\). This catches errors that may not be apparent from inspecting a single simulated matrix.

Common Mistake: Reversing a Factor or Its Transpose

A frequent error is to use \(AZB\) instead of \(AZB^{\mathsf T}\), or to put the column factor on the left. In general, \(AZB\) does not produce column covariance \(V=BB^{\mathsf T}\): its covariance uses \(B^{\mathsf T}B\) on the column side. Those matrices need not agree for a nonsymmetric Cholesky factor. The correction is to keep the factor and transpose paired as \(B\) and \(B^{\mathsf T}\) in the product \(AZB^{\mathsf T}\).

Another tempting but incorrect check is to claim that the rows of \(X\) have covariance \(U\) and the columns have covariance \(V\) without qualification. The precise entrywise statement is \(\operatorname{Cov}(X_{ij},X_{k\ell})=U_{ik}V_{j\ell}\). For example, the covariance between two entries in column \(j\) is \(V_{jj}U_{ik}\); the column scale \(V_{jj}\) is part of the result.

Key takeaway: To simulate \(MN_{r,c}(M,U,V)\), generate independent standard normal entries in \(Z\) and set \(X=M+AZB^{\mathsf T}\), where \(AA^{\mathsf T}=U\) and \(BB^{\mathsf T}=V\). The construction has vectorized covariance \(V\otimes U\), or, entry by entry, covariance \(U_{ik}V_{j\ell}\).

Check Your Understanding

Use the stated row-first, column-second covariance convention.

What are the dimensions of \(A\), \(Z\), and \(B\) when \(X\) is \(r\times c\)?
For \(X\sim MN_{r,c}(M,U,V)\), express \(\operatorname{Cov}(X_{ij},X_{k\ell})\) and \(\operatorname{Var}(X_{ij})\).
Why is \(B^{\mathsf T}\), rather than \(B\), used on the right in the Cholesky construction?
If \(U_{12}=0\), what is \(\operatorname{Cov}(X_{1j},X_{2\ell})\) for arbitrary \(j,\ell\)?
Does one simulated matrix have to exhibit the target covariance? Explain the difference between a draw and a covariance across repeated draws.