Matrix normal distribution

The probability density function for the random matrix X (n × p) that follows the matrix normal distribution ${\displaystyle {\mathcal {MN}}_{n,p}(\mathbf {M} ,\mathbf {U} ,\mathbf {V} )}$ has the form:

${\displaystyle p(\mathbf {X} \mid \mathbf {M} ,\mathbf {U} ,\mathbf {V} )={\frac {\exp \left(-{\frac {1}{2}}\,\mathrm {tr} \left[\mathbf {V} ^{-1}(\mathbf {X} -\mathbf {M} )^{T}\mathbf {U} ^{-1}(\mathbf {X} -\mathbf {M} )\right]\right)}{(2\pi )^{np/2}|\mathbf {V} |^{n/2}|\mathbf {U} |^{p/2}}}}$

where ${\displaystyle \mathrm {tr} }$ denotes trace and M is n × p, U is n × n and V is p × p, and the density is understood as the probability density function with respect to the standard Lebesgue measure in ${\displaystyle \mathbb {R} ^{n\times p}}$, i.e.: the measure corresponding to integration with respect to ${\displaystyle dx_{11}dx_{21}\dots dx_{n1}dx_{12}\dots dx_{n2}\dots dx_{np}}$.

The matrix normal is related to the multivariate normal distribution in the following way:

${\displaystyle \mathbf {X} \sim {\mathcal {MN}}_{n\times p}(\mathbf {M} ,\mathbf {U} ,\mathbf {V} ),}$

if and only if

${\displaystyle \mathrm {vec} (\mathbf {X} )\sim {\mathcal {N}}_{np}(\mathrm {vec} (\mathbf {M} ),\mathbf {V} \otimes \mathbf {U} )}$

where ${\displaystyle \otimes }$ denotes the Kronecker product and ${\displaystyle \mathrm {vec} (\mathbf {M} )}$ denotes the vectorization of ${\displaystyle \mathbf {M} }$.

If ${\displaystyle \mathbf {X} \sim {\mathcal {MN}}_{n\times p}(\mathbf {M} ,\mathbf {U} ,\mathbf {V} )}$, then we have the following properties:[1][2]

Let's imagine a sample of n independent p-dimensional random variables identically distributed according to a multivariate normal distribution:

${\displaystyle \mathbf {Y} _{i}\sim {\mathcal {N}}_{p}({\boldsymbol {\mu }},{\boldsymbol {\Sigma }}){\text{ with }}i\in \{1,\ldots ,n\}}$.

When defining the n × p matrix ${\displaystyle \mathbf {X} }$ for which the ith row is ${\displaystyle \mathbf {Y} _{i}}$, we obtain:

${\displaystyle \mathbf {X} \sim {\mathcal {MN}}_{n\times p}(\mathbf {M} ,\mathbf {U} ,\mathbf {V} )}$

where each row of ${\displaystyle \mathbf {M} }$ is equal to ${\displaystyle {\boldsymbol {\mu }}}$, that is ${\displaystyle \mathbf {M} =\mathbf {1} _{n}\times {\boldsymbol {\mu }}^{T}}$, ${\displaystyle \mathbf {U} }$ is the n × n identity matrix, that is the rows are independent, and ${\displaystyle \mathbf {V} ={\boldsymbol {\Sigma }}}$.

Given k matrices, each of size n × p, denoted ${\displaystyle \mathbf {X} _{1},\mathbf {X} _{2},\ldots ,\mathbf {X} _{k}}$, which we assume have been sampled i.i.d. from a matrix normal distribution, the maximum likelihood estimate of the parameters can be obtained by maximizing:

${\displaystyle \prod _{i=1}^{k}{\mathcal {MN}}_{n\times p}(\mathbf {X} _{i}\mid \mathbf {M} ,\mathbf {U} ,\mathbf {V} ).}$

The solution for the mean has a closed form, namely

${\displaystyle \mathbf {M} ={\frac {1}{k}}\sum _{i=1}^{k}\mathbf {X} _{i}}$

but the covariance parameters do not. However, these parameters can be iteratively maximized by zero-ing their gradients at:

${\displaystyle \mathbf {U} ={\frac {1}{kp}}\sum _{i=1}^{k}(\mathbf {X} _{i}-\mathbf {M} )\mathbf {V} ^{-1}(\mathbf {X} _{i}-\mathbf {M} )^{T}}$

and

${\displaystyle \mathbf {V} ={\frac {1}{kn}}\sum _{i=1}^{k}(\mathbf {X} _{i}-\mathbf {M} )^{T}\mathbf {U} ^{-1}(\mathbf {X} _{i}-\mathbf {M} ),}$

See for example [3] and references therein. The covariance parameters are non-identifiable in the sense that for any scale factor, s>0, we have:

${\displaystyle {\mathcal {MN}}_{n\times p}(\mathbf {X} \mid \mathbf {M} ,\mathbf {U} ,\mathbf {V} )={\mathcal {MN}}_{n\times p}(\mathbf {X} \mid \mathbf {M} ,s\mathbf {U} ,{\tfrac {1}{s}}\mathbf {V} ).}$

Sampling from the matrix normal distribution is a special case of the sampling procedure for the multivariate normal distribution. Let ${\displaystyle \mathbf {X} }$ be an n by p matrix of np independent samples from the standard normal distribution, so that

${\displaystyle \mathbf {X} \sim {\mathcal {MN}}_{n\times p}(\mathbf {0} ,\mathbf {I} ,\mathbf {I} ).}$

Then let

${\displaystyle \mathbf {Y} =\mathbf {M} +\mathbf {A} \mathbf {X} \mathbf {B} ,}$

so that

${\displaystyle \mathbf {Y} \sim {\mathcal {MN}}_{n\times p}(\mathbf {M} ,\mathbf {AA} ^{T},\mathbf {B} ^{T}\mathbf {B} ),}$

where A and B can be chosen by Cholesky decomposition or a similar matrix square root operation.

Dawid (1981) provides a discussion of the relation of the matrix-valued normal distribution to other distributions, including the Wishart distribution, inverse-Wishart distribution and matrix t-distribution, but uses different notation from that employed here.

• Multivariate normal distribution

• Dawid, A.P. (1981). "Some matrix-variate distribution theory: Notational considerations and a Bayesian application". Biometrika. 68 (1): 265–274. doi:10.1093/biomet/68.1.265. JSTOR 2335827. MR 0614963.
• Dutilleul, P (1999). "The MLE algorithm for the matrix normal distribution". Journal of Statistical Computation and Simulation. 64 (2): 105–123. doi:10.1080/00949659908811970.
• Arnold, S.F. (1981), The theory of linear models and multivariate analysis, New York: John Wiley & Sons, ISBN 0471050652