matrix normal distribution

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

:$$

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 $\backslash 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 $\backslash mathbb\{R\}^\{n\backslash times\; p\}$, i.e.: the measure corresponding to integration with respect to $dx\_\{11\}\; dx\_\{21\}\backslash dots\; dx\_\{n1\}\; dx\_\{12\}\backslash dots\; dx\_\{n2\}\backslash dots\; dx\_\{np\}$.

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

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

if and only if

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

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

The equivalence between the above matrix normal and multivariate normal density functions can be shown using several properties of the trace and Kronecker product, as follows. We start with the argument of the exponent of the matrix normal PDF:

:$\backslash begin\{align\}$

&\;\;\;\;-\frac12\text{tr}\left[ \mathbf{V}^{-1} (\mathbf{X} - \mathbf{M})^{T} \mathbf{U}^{-1} (\mathbf{X} - \mathbf{M}) \right]\\

&= -\frac12\text{vec}\left(\mathbf{X} - \mathbf{M}\right)^T

\text{vec}\left(\mathbf{U}^{-1} (\mathbf{X} - \mathbf{M}) \mathbf{V}^{-1}\right) \\

&= -\frac12\text{vec}\left(\mathbf{X} - \mathbf{M}\right)^T

\left(\mathbf{V}^{-1}\otimes\mathbf{U}^{-1}\right)\text{vec}\left(\mathbf{X} - \mathbf{M}\right) \\

&= -\frac12\left[\text{vec}(\mathbf{X}) - \text{vec}(\mathbf{M})\right]^T

\left(\mathbf{V}\otimes\mathbf{U}\right)^{-1}\left[\text{vec}(\mathbf{X}) - \text{vec}(\mathbf{M})\right]

\end{align}

which is the argument of the exponent of the multivariate normal PDF with respect to Lebesgue measure in $\backslash mathbb\{R\}^\{n\; p\}$. The proof is completed by using the determinant property: $|\backslash mathbf\{V\}\backslash otimes\; \backslash mathbf\{U\}|\; =\; |\backslash mathbf\{V\}|^n\; |\backslash mathbf\{U\}|^p.$

If $\backslash mathbf\{X\}\; \backslash sim\; \backslash mathcal\{MN\}\_\{n\backslash times\; p\}(\backslash mathbf\{M\},\; \backslash mathbf\{U\},\; \backslash mathbf\{V\})$, then we have the following properties:{{cite book|author1=A K Gupta|author2=D K Nagar|title=Matrix Variate Distributions|url=https://books.google.com/books?id=PQOYnT7P1loC|access-date=23 May 2014|date=22 October 1999|publisher=CRC Press|isbn=978-1-58488-046-2|chapter=Chapter 2: MATRIX VARIATE NORMAL DISTRIBUTION}}{{cite journal|last=Ding|first=Shanshan|author2=R. Dennis Cook|title=Dimension folding PCA and PFC for matrix-valued predictors|journal=Statistica Sinica|date=2014|volume=24|issue=1|pages=463–492|jstor=26432553}}

The mean, or expected value is:

:$E[\backslash mathbf\{X\}]\; =\; \backslash mathbf\{M\}$

and we have the following second-order expectations:

:$E[(\backslash mathbf\{X\}\; -\; \backslash mathbf\{M\})(\backslash mathbf\{X\}\; -\; \backslash mathbf\{M\})^\{T\}]$

= \mathbf{U}\operatorname{tr}(\mathbf{V})

:$E[(\backslash mathbf\{X\}\; -\; \backslash mathbf\{M\})^\{T\}\; (\backslash mathbf\{X\}\; -\; \backslash mathbf\{M\})]$

= \mathbf{V}\operatorname{tr}(\mathbf{U})

where $\backslash operatorname\{tr\}$ denotes trace.

More generally, for appropriately dimensioned matrices A,B,C:

:$\backslash begin\{align\}$

E[\mathbf{X}\mathbf{A}\mathbf{X}^{T}]

&= \mathbf{U}\operatorname{tr}(\mathbf{A}^T\mathbf{V}) + \mathbf{MAM}^T\\

E[\mathbf{X}^T\mathbf{B}\mathbf{X}]

&= \mathbf{V}\operatorname{tr}(\mathbf{U}\mathbf{B}^T) + \mathbf{M}^T\mathbf{BM}\\

E[\mathbf{X}\mathbf{C}\mathbf{X}]

&= \mathbf{V}\mathbf{C}^T\mathbf{U} + \mathbf{MCM}

\end{align}

Transpose transform:

:$\backslash mathbf\{X\}^T\; \backslash sim\; \backslash mathcal\{MN\}\_\{p\backslash times\; n\}(\backslash mathbf\{M\}^T,\; \backslash mathbf\{V\},\; \backslash mathbf\{U\})$

Linear transform: let D (r-by-n), be of full rank r ≤ n and C (p-by-s), be of full rank s ≤ p, then:

:$\backslash mathbf\{DXC\}\backslash sim\; \backslash mathcal\{MN\}\_\{r\backslash times\; s\}(\backslash mathbf\{DMC\},\; \backslash mathbf\{DUD\}^T,\; \backslash mathbf\{C\}^T\backslash mathbf\{VC\})$

The product of two matrix normal distributions

:$\backslash mathcal\{MN\}(\backslash mathbf\{M\_1\},\; \backslash mathbf\{U\_1\},\; \backslash mathbf\{V\_1\})\backslash cdot\; \backslash mathcal\{MN\}(\backslash mathbf\{M\_2\},\; \backslash mathbf\{U\_2\},\; \backslash mathbf\{V\_2\})\; \backslash propto\; \backslash mathcal\{N\}(\backslash mu\_c,\; \backslash Sigma\_c)$

is proportional to a normal distribution with parameters:

:$\backslash Sigma\_c\; =\; (V\_1^\{-1\}\; \backslash otimes\; U\_1^\{-1\}\; +\; V\_2^\{-1\}\; \backslash otimes\; U\_2^\{-1\})^\{-1\},$

:$\backslash mu\_c\; =\; \backslash Sigma\_c\; \backslash big((V\_1^\{-1\}\; \backslash otimes\; U\_1^\{-1\})\; \backslash operatorname\{vec\}(M\_1)\; +\; (V\_2^\{-1\}\; \backslash otimes\; U\_2^\{-1\})\backslash operatorname\{vec\}(M\_2)\backslash big).$

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

:$\backslash mathbf\{Y\}\_i\; \backslash sim\; \backslash mathcal\{N\}\_p(\{\backslash boldsymbol\; \backslash mu\},\; \{\backslash boldsymbol\; \backslash Sigma\})\; \backslash text\{\; with\; \}\; i\; \backslash in\; \backslash \{1,\backslash ldots,n\backslash \}$.

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

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

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

Given k matrices, each of size n × p, denoted $\backslash mathbf\{X\}\_1,\; \backslash mathbf\{X\}\_2,\; \backslash ldots,\; \backslash 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:

:$$

\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

:$$

\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:

:$$

\mathbf{U} = \frac{1}{kp} \sum_{i=1}^k(\mathbf{X}_i-\mathbf{M})\mathbf{V}^{-1}(\mathbf{X}_i-\mathbf{M})^T

and

:$$

\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 {{cite arXiv| last1=Glanz|first1=Hunter |last2=Carvalho|first2=Luis |title=An Expectation-Maximization Algorithm for the Matrix Normal Distribution |year=2013 |class=stat.ME |eprint=1309.6609}} and references therein. The covariance parameters are non-identifiable in the sense that for any scale factor, s>0, we have:

:$$

\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 $\backslash mathbf\{X\}$ be an n by p matrix of np independent samples from the standard normal distribution, so that

:$$

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

Then let

:$$

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

so that

:$$

\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

{{reflist}}

• {{cite journal

|last=Dawid |first=A.P. |author-link=Philip Dawid

|year=1981

|title=Some matrix-variate distribution theory: Notational considerations and a Bayesian application

|journal=Biometrika

|volume=68 |issue=1 |pages=265–274

|doi=10.1093/biomet/68.1.265 |mr=614963 | jstor = 2335827

}}

• {{cite journal

|last=Dutilleul |first=P

|year=1999

|title=The MLE algorithm for the matrix normal distribution

|journal=Journal of Statistical Computation and Simulation

|volume=64 |issue=2 |pages=105–123

|doi=10.1080/00949659908811970

}}

• {{Citation

|last=Arnold |first=S.F.

|title=The theory of linear models and multivariate analysis

|publisher=John Wiley & Sons

|place=New York

|year=1981

|isbn=0471050652

}}

{{ProbDistributions|multivariate}}

Category:Random matrices

Category:Continuous distributions

Category:Multivariate continuous distributions