Why the matrix normal distribution is not estimated with this "method of moments"?

Question

According to the Wikipedia Page, a random matrix $\bf{X}\in \mathbb{R}^{p\times q}$ follows a matrix normal distribution $\cal{MN}(\bf M, \bf U, \bf V)$ means that $$ \text{vec}(\bf X) \sim \cal N ( \mathrm{vec} (\bf M), \bf V \otimes \bf U),$$ where $\bf M \in \mathbb{R}^{p\times q}$, $\bf U \in \mathbb{R}^{p\times p}$, and $\bf V \in \mathbb{R}^{q\times q}$.

Suppose we have observed matrix-valued data $\bf X_1, \dots, \bf X_n$, which are assumed to be i.i.d. realizations of $\cal{MN}(\bf M, \bf U, \bf V)$. Then, the MLE has the estimate $\hat{\bf M}$ as the usual sample mean, and $\bf U$, $\bf V$ given by the following iterative procedure: $$ \begin{aligned} \hat{\bf U} = \frac{1}{nq} \sum_{i=1}^n (\bf X_i - \bf M) \hat{\bf V}^{-1} (\bf X_i - \bf M)^\top, \\ \hat{\bf V} = \frac{1}{np} \sum_{i=1}^n (\bf X_i - \bf M)^\top \hat{\bf U}^{-1} (\bf X_i - \bf M). \end{aligned} $$

I understand the MLE framework, but I wonder why this method does not work: Notice from the same Wikipedia page that $$\begin{aligned} \bf E[(\bf X - \bf M) (\bf X - \bf M)^\top] = \bf U\, \mathrm{tr}(\bf V), \\ \bf E[(\bf X - \bf M)^\top (\bf X - \bf M)] = \bf V\, \mathrm{tr}(\bf U). \end{aligned}$$ As $\bf U$ and $\bf V$ are not identifiable up to a scaling factor, we pose the additional restriction $\mathrm{tr}(\bf V) = q$. Then, I think it is plausible to estimate by: $$\begin{aligned} \hat{\bf U} = \frac{1}{nq} \sum_{i=1}^n (\bf X_i - \bf M) (\bf X_i - \bf M)^\top, \\ \hat{\bf V} = \frac{1}{n\operatorname{tr}(\hat{\bf U})} \sum_{i=1}^n (\bf X_i - \bf M)^\top (\bf X_i - \bf M). \end{aligned}$$ However, the second approach is not aligned with the MLE. Where did it go wrong?

Answer 1

The second approach to estimating the parameters of a matrix normal distribution essentially attempts to bypass the complexities of maximum likelihood estimation (MLE) by relying solely on the raw second moments of the data and a normalization constraint on the trace of $\mathbf{V}$, where $\mathrm{tr}(\mathbf{V}) = q$. The matrix normal distribution, denoted as $\mathcal{MN}(\mathbf{M}, \mathbf{U}, \mathbf{V})$, has a density function dependent on $\mathbf{U}$ and $\mathbf{V}$ through their Kronecker product, $\mathbf{V} \otimes \mathbf{U}$. However, $\mathbf{U}$ and $\mathbf{V}$ are not identifiable due to a scaling ambiguity; replacing $\mathbf{U}$ with $a\mathbf{U}$ and $\mathbf{V}$ with $\frac{1}{a}\mathbf{V}$ leaves the distribution unchanged for any scalar $a>0$. To address this, the constraint $\mathrm{tr}(\mathbf{V}) = q$ is often imposed. One might be tempted to use the following moment conditions: $\mathbb{E}[(\mathbf{X}-\mathbf{M})(\mathbf{X}-\mathbf{M})^\top] = \mathbf{U}\,\mathrm{tr}(\mathbf{V})$ and $\mathbb{E}[(\mathbf{X}-\mathbf{M})^\top(\mathbf{X}-\mathbf{M})] = \mathbf{V}\,\mathrm{tr}(\mathbf{U})$. From these, one might derive estimates for $\mathbf{U}$ and $\mathbf{V}$ directly from the sample covariances, leading to $\hat{\mathbf{U}} = \frac{1}{nq}\sum_{i=1}^n (\mathbf{X}_i - \mathbf{M})(\mathbf{X}_i - \mathbf{M})^\top$ and $\hat{\mathbf{V}} = \frac{1}{n \operatorname{tr}(\hat{\mathbf{U}})}\sum_{i=1}^n (\mathbf{X}_i - \mathbf{M})^\top(\mathbf{X}_i - \mathbf{M})$. This intuitive approach, however, does not yield the MLE. The true MLE equations arise from maximizing the joint likelihood of the data, which involves the inverse covariance structure, $\mathbf{V}^{-1}\otimes \mathbf{U}^{-1}$. This results in a set of coupled nonlinear equations: $$\hat{\mathbf{U}} = \frac{1}{nq}\sum_{i=1}^n (\mathbf{X}_i - \mathbf{M}) \hat{\mathbf{V}}^{-1} (\mathbf{X}_i - \mathbf{M})^\top$$ and $$\hat{\mathbf{V}} = \frac{1}{np}\sum_{i=1}^n (\mathbf{X}_i - \mathbf{M})^\top \hat{\mathbf{U}}^{-1} (\mathbf{X}_i - \mathbf{M}).$$ These equations contain $\hat{\mathbf{U}}^{-1}$ and $\hat{\mathbf{V}}^{-1}$ within the summations, making them nonlinear and requiring an iterative solution. The simple second-moment approach bypasses these nonlinearities by simply plugging in sample moments and using the normalization constraint. It does not consider the effect of the inverse covariance matrices, $\mathbf{U}^{-1}$ and $\mathbf{V}^{-1}$, on the likelihood, and as a result, it cannot produce the MLE. The iterative procedure is necessary because the maximum likelihood conditions naturally involve these inverses, which cannot be estimated directly without iteration. Thus, only by solving these full MLE equations can we be sure that our estimates align with the true maximum likelihood solution.