# Covarince matrix has to be be rank-deficient to maximize multivariate gaussian likelihood?

I am facing a similar situation as that addressed in this question, and the accepted answer there has helped me a lot, but I need to resolve a doubt.

The accepted answer draws on the excellent resource "Matrix Cookbook" to show that

$$\frac{\partial \mathbf{L}}{\partial \mathbf{\Sigma}}= -1/2\left(\mathbf{\Sigma^{-1}-\Sigma^{-1}(y-\mu)(y-\mu)'\Sigma^{-1}}\right)$$

where $\mathbf{L}$ is the log likelihood of the gaussian vector $\mathbf{y}$ with covariance matrix $\mathbf{\Sigma}$ and mean $\mathbf{\mu}$

If I'm not mistaken, to solve for the $\mathbf{\Sigma}$ that maximizes $\mathbf{L}$, one would then set $\frac{\partial \mathbf{L}}{\partial \mathbf{\Sigma}}=0$ and get

$$\mathbf{\Sigma^{-1}=\Sigma^{-1}(y-\mu)(y-\mu)'\Sigma^{-1}}$$

Pre- and post-multiplying by $\mathbf{\Sigma}$ gives

$$\mathbf{\Sigma=(y-\mu)(y-\mu)'}$$

Which brings me to my doubt. The RHS of this is an outer product with determinant 0 and thus not invertible. But the covariance matrix $\mathbf{\Sigma}$ on the LHS must be invertible. Likewise, the intermediate equation

$$\mathbf{I=\Sigma^{-1}(y-\mu)(y-\mu)'}$$

seems to be a contradiction since $\mathbf{[(y-\mu)(y-\mu)']^{-1}}$ does not exist.

Anyways, what this seems to say to me is that the condition for maximum $\mathbf{L}$ requires that $\mathbf{\Sigma}$ be rank deficient with determinant 0, in which case it could not really be called a covariance matrix, and $\mathbf{L}$ would be undefined at its maximum.

And yet the answerer says that he/she uses these formulae all the time for ML parameter estimation, so I guess I am missing something. Please help.

The derivative you describe is with respect to the density (and hence the likelihood) of a single observation. If instead you suppose that you have $n$ observations independently drawn from the same distribution, the full likelihood is $$L(\theta|y_1, \dots, y_n)=\prod_{i=1}^n p(y_i|\theta),$$ by which $$\log L(\theta|y_1, \dots, y_n)=\sum_{i=1}^n \log p(y_i|\theta).$$
The derivative in your post is $$\frac{\partial \log p(y_i|\mu, \Sigma)}{\partial \Sigma}$$ which is summed over $i=1, \dots, n$ yielding $$\frac{\partial \log L(\mu, \Sigma|y_1, \dots, y_n)}{\partial \Sigma}=\sum_{i=1}^n \frac{\partial \log p(y_i|\mu, \Sigma)}{\partial \Sigma}=-\frac{n}{2}\Sigma^{-1}+\frac{1}{2}\Sigma^{-1}\left[\sum_{i=1}^n(y_i-\mu)(y_i-\mu)'\right]\Sigma^{-1}.$$ Solving the same way you did results in $$\Sigma=\frac{1}{n}\sum_{i=1}^n(y_i-\mu)(y_i-\mu)'.$$
If you now let $Z=(z_1, \dots, z_n)'$ (a $n\times k$ matrix) where $z_i=y_i-\mu$, you may equivalently express this as $$\Sigma=\frac{1}{n}Z'Z.$$ Note here that the rank of $Z'Z$ is at most $k$ as it's a $k\times k$ matrix. Further, the rank of $Z$ is at most $\min(k, n)$ (i.e. the smallest of the column and row rank). In fact, with the rows of $Z$ being independent draws from a normal distribution the rank of $Z$ will be exactly $\min(n, k)$ with probability one. Thus, for the example in your original post, you have $n=1$ so that $rank(Z)=1$ and hence $rank(Z'Z)=1$ and the estimate of the covariance matrix is singular. For non-singularity, you will need $n\geq k$.