# Equivalent method for maximum likelihood estimation of covariance parameters

My goal is to estimate the parameters of a covariance matrix $\Omega$, by maximizing the following log-likelihood function:

$$\log L(\vec\tau, \rho, \sigma \mid W, X) = -m\ln(\left | \Omega \right |) - \operatorname{Tr}(X^T\Omega(\vec\tau, \rho, \sigma)^{-1}X)$$

with $\Omega = (1-\rho)\operatorname{diag}(\vec\tau^2) + \rho\vec\tau\vec\tau^T + \sigma^2WW^T$

where $X$ and $W$ are known matrices in $\mathbb{R}^{m \times n}$ and $\mathbb{R}^{n \times k}$, respectively, and thus $\Omega$ is in $\mathbb{R}^{n\times n}$. To clarify my notation: by $\text{diag}(\vec\tau^2)$ I mean the matrix which along its diagonal has the elements of the vector $\tau$ squared, and off-diagonal entries equal to 0 (apologies if this is not the most conventional notation).

The problem is that the maximization of this likelihood has to be done numerically (at least, I have been unable to derive a closed-form expression for any of the parameters, but I would be very happy to be proven wrong), and each iteration of the optimization algorithm requires the inversion of $\Omega$. So I'm wondering whether there isn't a faster way of doing this optimization. Specifically, I observe that if I didn't assume any structure for $\Omega$ (i.e. if each of its entries were unconstrained rather than being a function of some set of parameters), the maximum likelihood estimate (MLE) of $\Omega$ would be equal to:

$$\hat{\Omega}_{ML} = S = \frac{1}{m}XX^T$$

So my question is: would it be equivalent, instead of maximizing the full log-likelihood of my parameters, to minimize the squared error between $\Omega(\vec\tau, \rho, \sigma)$ and $S$? That is, to compute:

$$\hat{\theta} = \arg\min_\theta \left \| S-\Omega(\theta)\right \|_F^2$$

With $\theta = \{\vec\tau, \rho, \sigma\}$

And if so, is there an easy way to prove this?

This would be very helpful to me as it would probably speed up my data analysis by several orders of magnitude.

• I'll probably change my mind about this a few minutes after starting to think about it, but I might begin by asking whether the methods on this page could be adapted to the present situation: en.wikipedia.org/wiki/… Dec 13, 2013 at 1:20
• You wrote "$X$ and $W$ are known matrices". The expression $L(\bullet\mid W,X)$ makes me think they're your data, i.e. realizations of random vectors, where you know a parametrized family of probability distributions, from one of which they were drawn. Is $m$ your sample size? For $X$, at least? Dec 13, 2013 at 1:36
• Well, strictly speaking I have a matrix $B$ in $\mathbb{R}^{m \times n}$ (where $m$ is indeed the sample size) which is my actual dataset. $W$ is a parameter that I estimate independently, and $X = B - WC$, where $C$ is an $m \times k$ matrix which forms a basis as a function of a variable $s$ (which in this (training) stage of the analysis is known, but later I will try to estimate given some test data and the parameter estimates). I was worried that all of that background might only make things more confusing so I left it out - do you think it matters for my question? Dec 13, 2013 at 14:04
• I'm not sure it matters, but I wondered to what extent this resembles things I've seen before. Dec 13, 2013 at 17:51