Matrix-Normal Errors and Separable Covariance

Why Model Errors as a Matrix?

In a multivariate regression, the error term is naturally an \(n\times p\) matrix: rows correspond to cases and columns to responses. The General Multivariate Linear Model writes the response matrix as \(Y=XB+U\), where \(U\) is the error matrix. A covariance model for \(U\) must describe dependence both across rows and across responses. The matrix-normal model does this by assigning one covariance factor to each direction.

This structure is useful when dependence among cases has a common form across responses, and dependence among responses has a common form across cases. It is also a substantive restriction: not every covariance matrix for the \(np\) entries of \(U\) can be expressed using just these two factors.

Definition: Let \(U\) be an \(n\times p\) random matrix. It has a matrix-normal distribution with mean \(M\), row covariance \(\Omega\), and response covariance \(\Sigma\), written \(U\sim MN_{n,p}(M,\Omega,\Sigma)\), if \(\operatorname{vec}(U)\sim N_{np}(\operatorname{vec}(M),\Sigma\otimes\Omega)\). Here \(\operatorname{vec}(U)\) stacks the columns of \(U\) in order; \(\Omega\) is \(n\times n\), and \(\Sigma\) is \(p\times p\). For the usual nonsingular density, \(\Omega\) and \(\Sigma\) are symmetric positive definite.

The order of the covariance factors depends on the vectorization convention. With column stacking, the response covariance appears first in \(\Sigma\otimes\Omega\). In the regression error model, \(M=0\), so the convention is \(U\sim MN_{n,p}(0,\Omega,\Sigma)\). The entrywise covariance formula makes the roles of the factors explicit.

Theorem: If \(U\sim MN_{n,p}(M,\Omega,\Sigma)\), then for all row indices \(i,k\) and response indices \(j,l\), \(\operatorname{Cov}(U_{ij},U_{kl})=\Omega_{ik}\Sigma_{jl}\). In particular, \(\operatorname{Var}(U_{ij})=\Omega_{ii}\Sigma_{jj}\).

Derivation. Under column stacking, the entry \(U_{ij}\) is coordinate \(i+(j-1)n\) of \(\operatorname{vec}(U)\). The covariance matrix of that vector is \(\Sigma\otimes\Omega\). By the definition of the Kronecker product, its block for response columns \(j,l\) is \(\Sigma_{jl}\Omega\). The entry of that block in row \(i\), column \(k\) is \(\Sigma_{jl}\Omega_{ik}\), which is exactly the covariance of \(U_{ij}\) and \(U_{kl}\). Setting \(k=i\) and \(l=j\) gives the stated variance. \(\square\)

Worked Example: Covariances Between Error Entries

Suppose \(n=p=2\), with \(\Omega=\begin{pmatrix}1&0.4\\0.4&2\end{pmatrix}\) and \(\Sigma=\begin{pmatrix}4&1\\1&9\end{pmatrix}\). Both matrices are positive definite: their diagonal entries are positive, and their determinants are \(1.84\) and \(35\), respectively.

The variance of \(U_{12}\), the entry in row 1 and response column 2, is \(\Omega_{11}\Sigma_{22}=1(9)=9\). The variance of \(U_{21}\) is \(\Omega_{22}\Sigma_{11}=2(4)=8\). Their covariance is \(\Omega_{12}\Sigma_{21}=0.4(1)=0.4\), so their correlation is \(0.4/\sqrt{9\cdot8}\approx0.047\).

The formula also separates two kinds of association. Within row 1, the covariance between \(U_{11}\) and \(U_{12}\) is \(\Omega_{11}\Sigma_{12}=1\). Within response 1, the covariance between \(U_{11}\) and \(U_{21}\) is \(\Omega_{12}\Sigma_{11}=1.6\). These are different quantities because one involves response covariance and the other involves row covariance.

The Density and the Kronecker Covariance

When both covariance factors are positive definite, the matrix-normal distribution has a density with respect to Lebesgue measure on the \(np\) entries of \(U\). For a general mean matrix \(M\), it is

f(U)=\frac{\exp\!\left[-\frac12\operatorname{tr}\!\left\{\Sigma^{-1}(U-M)^{\mathsf T}\Omega^{-1}(U-M)\right\}\right]} {(2\pi)^{np/2}|\Omega|^{p/2}|\Sigma|^{n/2}}.

To see why this is the multivariate normal density in matrix form, put \(u=\operatorname{vec}(U-M)\). The quadratic form in the vector density is \(u^{\mathsf T}(\Sigma^{-1}\otimes\Omega^{-1})u\). The vectorization identity gives \(u^{\mathsf T}(\Sigma^{-1}\otimes\Omega^{-1})u =\operatorname{tr}\{\Sigma^{-1}(U-M)^{\mathsf T}\Omega^{-1}(U-M)\}\). Also, \(|\Sigma\otimes\Omega|=|\Sigma|^n|\Omega|^p\), so the multivariate normal normalizing factor becomes the denominator displayed above. Thus the density and the column-stacked definition agree.

Worked Example: Constructing the Vector Covariance

Let \(n=p=2\), and take \(\Omega=\begin{pmatrix}1&0.5\\0.5&2\end{pmatrix}\) and \(\Sigma=\begin{pmatrix}3&1\\1&2\end{pmatrix}\). Both are positive definite: their determinants are \(1.75\) and \(5\), and their diagonal entries are positive. For column stacking, \(\operatorname{vec}(U)=(U_{11},U_{21},U_{12},U_{22})^{\mathsf T}\).

The covariance is \(\Sigma\otimes\Omega\). Writing out its four \(2\times2\) blocks gives

\Sigma\otimes\Omega = \begin{pmatrix} 3\Omega&1\Omega\\ 1\Omega&2\Omega \end{pmatrix} = \begin{pmatrix} 3&1.5&1&0.5\\ 1.5&6&0.5&2\\ 1&0.5&2&1\\ 0.5&2&1&4 \end{pmatrix}.

For example, the diagonal entry for \(U_{21}\) is \(6=\Omega_{22}\Sigma_{11}\). The entry in row 2, column 3 of this covariance matrix is \(0.5=\Omega_{21}\Sigma_{12}\), the covariance of \(U_{21}\) and \(U_{12}\). The block form shows that every response-pair block is a scalar multiple of the same row covariance matrix \(\Omega\).

What Separability Requires

A general covariance matrix for \(\operatorname{vec}(U)\) is any symmetric positive-definite \(np\times np\) matrix. A separable covariance is a much more restricted one: it must have the form \(\Sigma\otimes\Omega\), with positive-definite factors of sizes \(p\times p\) and \(n\times n\). Equivalently, partition the covariance into \(p\times p\) blocks of size \(n\times n\), one block for each pair of response columns. For separability, every block \(C_{jl}\) must equal \(\Sigma_{jl}\Omega\). Thus all these blocks must be scalar multiples of one common matrix.

Result: A symmetric positive-definite covariance matrix \(C\) for column-stacked \(U\) is separable if and only if it can be partitioned into response-pair blocks \(C_{jl}\) for which there are positive-definite matrices \(\Omega\) and \(\Sigma\) satisfying \(C_{jl}=\Sigma_{jl}\Omega\) for every \(j,l\).

Justification. If \(C=\Sigma\otimes\Omega\), the definition of the Kronecker product says that its \((j,l)\) response-pair block is \(\Sigma_{jl}\Omega\). Conversely, if every block has this form, assembling the blocks in response-column order produces exactly \(\Sigma\otimes\Omega\). Positive definiteness of the factors ensures this is a nonsingular covariance matrix. The requirement that the same \(\Omega\) works for every block is the substantive restriction. \(\square\)

Worked Example: A Valid Covariance That Is Not Separable

For \(n=p=2\), consider a covariance matrix for the column-stacked vector \((U_{11},U_{21},U_{12},U_{22})^{\mathsf T}\): \(\operatorname{diag}(1,1,1,2)\). It is positive definite because all four diagonal entries are positive.

If it were \(\Sigma\otimes\Omega\), its diagonal entries, in order, would be \(\Omega_{11}\Sigma_{11}\), \(\Omega_{22}\Sigma_{11}\), \(\Omega_{11}\Sigma_{22}\), and \(\Omega_{22}\Sigma_{22}\). Consequently, the product of the first and fourth entries would have to equal the product of the second and third: \[ (\Omega_{11}\Sigma_{11})(\Omega_{22}\Sigma_{22}) =(\Omega_{22}\Sigma_{11})(\Omega_{11}\Sigma_{22}). \] But for the proposed diagonal covariance, these products are \(1\cdot2=2\) and \(1\cdot1=1\). The required equality fails, so this positive-definite covariance is not separable. Positive definiteness alone does not imply a matrix-normal covariance structure.

Scale Ambiguity and Interpretation

The two covariance factors are not separately determined by their Kronecker product unless a scale convention is imposed. For any \(c>0\), \[ \left(\frac{\Sigma}{c}\right)\otimes(c\Omega) =\Sigma\otimes\Omega. \] Thus the matrix-normal distribution of \(U\) is unchanged by this simultaneous rescaling. This is a parameterization issue, not a change in the overall covariance of the errors. If separate numerical values for the factors must be reported, a scale constraint, such as fixing the trace of one factor, can remove the ambiguity.

The model’s interpretation is that row and response effects combine multiplicatively: the covariance between two entries is the row association \(\Omega_{ik}\) times the response association \(\Sigma_{jl}\). In particular, an identity row factor gives zero covariance between distinct rows, while a nonidentity row factor permits dependence across rows. The separable model assumes that the pattern of response covariance is shared across row pairs, up to the corresponding row-covariance multiplier. It does not allow every pair of entries to have its own unconstrained covariance.

Common Mistake

A frequent error is to write \(\Omega\otimes\Sigma\) for the covariance of \(\operatorname{vec}(U)\) while using column stacking. That reverses the dimensions and generally gives the wrong covariance: the \((j,l)\) response block must be \(\Sigma_{jl}\Omega\), so the full matrix is \(\Sigma\otimes\Omega\). Another mistake is to treat separability as a mere relabeling of any covariance matrix. The counterexample above shows that a positive-definite covariance can fail the required block structure. Check both the vectorization order and the common-factor restriction before interpreting a proposed model.

Key takeaway: With column-stacked errors, \(U\sim MN_{n,p}(0,\Omega,\Sigma)\) means \(\operatorname{Cov}\{\operatorname{vec}(U)\}=\Sigma\otimes\Omega\) and \(\operatorname{Cov}(U_{ij},U_{kl})=\Omega_{ik}\Sigma_{jl}\). This compact covariance model is separable, and therefore more restrictive than an arbitrary covariance for all \(np\) entries.

Check Your Understanding

Answer using the column-stacking convention stated above.

For \(U\sim MN_{n,p}(0,\Omega,\Sigma)\), write the covariance between \(U_{ij}\) and \(U_{kl}\).
Why does the covariance of \(\operatorname{vec}(U)\) have the order \(\Sigma\otimes\Omega\), rather than \(\Omega\otimes\Sigma\)?
What common block structure must hold for a covariance matrix of \(\operatorname{vec}(U)\) to be separable?
Explain why multiplying \(\Omega\) by a positive constant and dividing \(\Sigma\) by that same constant leaves the distribution unchanged.
Can every symmetric positive-definite covariance matrix for \(\operatorname{vec}(U)\) be represented as \(\Sigma\otimes\Omega\)? Explain.