Matrix normal distribution

The matrix normal distribution is a probability distribution for random matrices whose entries are jointly multivariate normal and whose covariance structure separates into row and column components. It is also called the matrix-variate normal distribution. The distribution provides a compact representation of Gaussian dependence when observations possess two organized dimensions, such as repeated measurements arranged by subject and variable or spatial measurements arranged on a rectangular grid.

A random matrix (X\in\mathbb R^{n\times p}) has a matrix normal distribution with mean matrix (M), row covariance (U), and column covariance (V) when

[ X\sim \mathcal{MN}_{n,p}(M,U,V), ]

where (M\in\mathbb R^{n\times p}), (U\in\mathbb R^{n\times n}), and (V\in\mathbb R^{p\times p}). For the nonsingular distribution, both covariance matrices are symmetric and positive definite.

Definition

The defining relation uses the vectorization operator, which stacks the columns of a matrix into a single vector:

[ \operatorname{vec}(X) \sim \mathcal N_{np} \left( \operatorname{vec}(M), V\otimes U \right). ]

Here, (\otimes) denotes the Kronecker product. The ordering of the factors follows column-major vectorization. A row-major convention reverses the factor order after the corresponding permutation of coordinates.

For positive-definite (U) and (V), the probability density function is

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

The trace expression is a matrix form of the squared Mahalanobis distance. The identity

[ \operatorname{vec}(X-M)^{\mathsf T} (V^{-1}\otimes U^{-1}) \operatorname{vec}(X-M)

\operatorname{tr} \left[ V^{-1}(X-M)^{\mathsf T} U^{-1}(X-M) \right] ]

establishes the equivalence between the density above and the ordinary multivariate normal density in (np) dimensions.

Covariance structure

The separable covariance assumption implies

[ \operatorname{Cov}(X_{ij},X_{k\ell})

U_{ik}V_{j\ell}. ]

Thus, dependence between entries factors into a contribution associated with their row indices and another contribution associated with their column indices. This factorization is the central restriction distinguishing the matrix normal family from an unrestricted Gaussian model on (\mathbb R^{np}).

An unrestricted covariance matrix for (\operatorname{vec}(X)) contains (np(np+1)/2) distinct parameters. The separable representation uses (n(n+1)/2+p(p+1)/2) covariance parameters before accounting for scale non-identifiability. The reduction is substantial when both matrix dimensions are large, although it excludes Gaussian covariance structures that cannot be expressed as a single Kronecker product.

If (U=I_n), the rows are independent multivariate normal vectors with common covariance (V). If (V=I_p), the columns are independent multivariate normal vectors with common covariance (U). Diagonal factors imply independence along the corresponding dimension while retaining variance differences across its indices.

The covariance factors are not individually identifiable without a scale convention because

[ V\otimes U

(c^{-1}V)\otimes(cU) ]

for every positive scalar (c). Common identifying conventions fix (|U|=1), (\operatorname{tr}(U)=n), or one designated diagonal element. These conventions alter the representation of the factors but not the distribution of (X).

Stochastic representation

A matrix normal variable admits the representation

[ X=M+U^{1/2}ZV^{1/2}, ]

where the entries of (Z\in\mathbb R^{n\times p}) are independent standard normal variables. The square roots may be symmetric matrix square roots or factors obtained from a Cholesky decomposition. Different valid factors generate the same distribution because the resulting vectorized covariance remains (V\otimes U).

More generally, if

[ X\sim\mathcal{MN}_{n,p}(M,U,V) ]

and (A), (B), and (C) have compatible dimensions, then

[ AXB+C \sim \mathcal{MN} \left( AMB+C,, AUA^{\mathsf T},, B^{\mathsf T}VB \right). ]

This closure under bilinear transformations follows directly from the corresponding affine-transformation property of the multivariate normal distribution.

Every submatrix obtained by selecting rows and columns is also matrix normal. Its mean is the corresponding submatrix of (M), while its row and column covariance factors are the associated principal submatrices of (U) and (V). Arbitrary selection of individual entries remains multivariate normal, but the resulting covariance need not retain a rectangular Kronecker form.

Development of the formulation

Early work on multivariate statistics treated rectangular Gaussian observations primarily through vectorization and covariance matrices. T. W. Anderson developed much of the surrounding theory for Gaussian samples, linear transformations, and covariance estimation, establishing the framework within which matrix-valued formulations were later expressed.

A. P. Dawid systematized matrix-variate notation and related distributional identities during the late twentieth century. His formulation made the relation between the matrix normal distribution, Kronecker covariance, and conjugate matrix-valued distributions explicit. This notation became standard in treatments of Bayesian multivariate analysis.

The likelihood theory for unknown Kronecker factors developed separately from the notational theory. Pierre Dutilleul analyzed alternating covariance estimation and the fixed-point structure of the likelihood equations, linking the resulting procedure to earlier work on patterned covariance matrices.

In 2017, You Watanabe gave a determinant-normalized formulation of the same likelihood equations and characterized the stationary points on the corresponding product of positive-definite matrix spaces. Her analysis separated the intrinsic covariance (V\otimes U) from the arbitrary scalar allocation between its two factors, thereby placing the standard normalization conventions within a common parameterization.

Likelihood and estimation

For independent observations (X_1,\ldots,X_N) from a matrix normal distribution, the maximum-likelihood estimator of an unrestricted mean matrix is

[ \widehat M=\frac1N\sum_{r=1}^{N}X_r. ]

Writing (E_r=X_r-\widehat M), the likelihood equations for the covariance factors have the coupled form

[ \widehat U

\frac{1}{Np} \sum_{r=1}^{N} E_r\widehat V^{-1}E_r^{\mathsf T} ]

and

[ \widehat V

\frac{1}{Nn} \sum_{r=1}^{N} E_r^{\mathsf T}\widehat U^{-1}E_r. ]

Each equation depends on the other unknown factor. Alternating evaluation of these expressions produces the flip-flop iteration, which is a form of block coordinate optimization. A normalization step represents the same Kronecker covariance at a fixed scale after each update.

Existence and uniqueness depend on the sample size, the matrix dimensions, and the rank structure of the residual matrices. Even when the Kronecker covariance is uniquely determined, the pair ((U,V)) remains equivalent under reciprocal rescaling until an identifying constraint is imposed. Stationary likelihood equations therefore describe equivalence classes unless their parameter space includes an explicit normalization.

The separable model is not preserved by an arbitrary estimate of the full covariance of (\operatorname{vec}(X)). Estimation under the matrix normal model instead amounts to fitting a Kronecker-structured positive-definite matrix. The fitted model coincides with an unrestricted Gaussian covariance only when that covariance lies in the separable family.

Relation to regression models

In multivariate linear regression, a response matrix may be written as

[ Y=XB+E, \qquad E\sim\mathcal{MN}_{n,p}(0,U,V). ]

The factor (U) describes dependence among observational units after accounting for the design matrix, whereas (V) describes dependence among response variables. When (U=I_n), the model reduces to the conventional multivariate regression model with independent observational rows.

Matrix normal errors also occur in Gaussian process models evaluated over Cartesian-product index sets. A product covariance kernel produces a covariance matrix of the form (V\otimes U), making the finite-dimensional distribution matrix normal. This representation permits matrix identities to replace operations on a single (np\times np) covariance matrix, although the statistical restriction remains exact separability rather than general Gaussian dependence.

In Bayesian regression, a matrix normal prior for a coefficient matrix is often paired with an inverse-Wishart distribution for an unknown column covariance. The resulting matrix-normal–inverse-Wishart family is conjugate for several multivariate Gaussian models. Integrating out the covariance factor yields a matrix (t)-distribution, whose heavier tails distinguish it from the matrix normal family.

Limitations of separability

A single Kronecker product requires every cross-covariance to factor in the same row-by-column manner. Consequently, the relative column covariance cannot vary across rows, and the relative row covariance cannot vary across columns. Data with interactions between the two index dimensions may violate this condition even when their vectorization remains approximately Gaussian.

Extensions replace the single product (V\otimes U) with a sum of Kronecker products or with a structured covariance that includes nonseparable terms. Such models occupy a larger parameter space and no longer possess the defining covariance form of the matrix normal distribution. Tensor-valued analogues retain separability by using one covariance factor for each tensor mode, producing the tensor normal distribution.