When the Gaussian Covariance MLE Exists

Why the Existence Question Matters

For a multivariate Gaussian model, the unrestricted maximum-likelihood estimate of the covariance matrix is based on the scatter of the observations around their sample mean. But a covariance estimate must be positive definite to belong to the usual Gaussian model. If the observations occupy a lower-dimensional subspace, the scatter matrix is singular, and the likelihood can increase without bound as the fitted covariance collapses in directions with no observed variation.

The key issue is therefore not just how to calculate the covariance estimate, but whether the data provide variation in every direction. We establish the exact rank condition for a finite maximum and then show when that condition holds with probability one.

Definition: For observations $X_1,\ldots,X_n$ in $\mathbb{R}^p$, define the sample mean and sample scatter matrix by $$ \bar X=\frac{1}{n}\sum_{i=1}^n X_i, \qquad S=\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)'. $$ The scatter matrix is positive definite if $a'Sa>0$ for every nonzero $a\in\mathbb{R}^p$. Equivalently, it has rank $p$. In the unrestricted Gaussian model with unknown mean and positive-definite covariance, the covariance MLE, when it exists, is $S/n$.

The divisor here is $n$, as required for maximum likelihood. The issue in this tutorial is whether $S/n$ lies inside the parameter space of positive-definite covariance matrices, not whether another covariance estimator might be defined when $S$ is singular.

Worked Examples: Checking the Scatter Rank

Worked Example: Three Noncollinear Observations in Two Dimensions

Take $p=2$ and observations $X_1=(0,0)'$, $X_2=(2,0)'$, and $X_3=(0,2)'$. Their mean is $\bar X=(2/3,2/3)'$. The scatter matrix is $$ S= \begin{pmatrix} 8/3 & -4/3\\ -4/3 & 8/3 \end{pmatrix}. $$ Its determinant is $(8/3)^2-(-4/3)^2=48/9=16/3>0$, so it is positive definite. The covariance MLE is $$ \hat\Sigma=\frac{S}{3} = \begin{pmatrix} 8/9 & -4/9\\ -4/9 & 8/9 \end{pmatrix}. $$ Its determinant is $48/81=16/27>0$, confirming that this estimate is in the positive-definite parameter space.

Worked Example: Collinear Observations and an Unbounded Likelihood

Now take $p=2$ and $X_1=(0,0)'$, $X_2=(1,1)'$, and $X_3=(2,2)'$. The mean is $(1,1)'$, and direct calculation gives $$ S= \begin{pmatrix} 2&2\\ 2&2 \end{pmatrix}. $$ This matrix has rank one. The unit vectors $v=(1,1)'/\sqrt{2}$ and $u=(1,-1)'/\sqrt{2}$ are eigenvectors of $S$ with eigenvalues $4$ and $0$, respectively. For $\varepsilon>0$, choose a covariance matrix with eigenvalue $4/3$ along $v$ and eigenvalue $\varepsilon$ along $u$. This covariance is positive definite for every $\varepsilon>0$.

At the sample mean, the log-likelihood, apart from its constant term, is $-\frac{n}{2}\log|\Sigma|-\frac{1}{2}\operatorname{tr}(\Sigma^{-1}S)$. For this choice, $|\Sigma|=4\varepsilon/3$ and $\operatorname{tr}(\Sigma^{-1}S)=4/(4/3)=3$. With $n=3$, the log-likelihood is therefore $$ -\frac{3}{2}\log(4\varepsilon/3)-\frac{3}{2}+\text{constant}. $$ As $\varepsilon$ decreases to zero, this tends to $+\infty$. No positive-definite covariance matrix maximizes the likelihood.

Worked Example: The Boundary Case of Four Points in Three Dimensions

Let $p=3$, $n=4$, and observe $0,e_1,e_2,e_3$, where $e_1,e_2,e_3$ are the coordinate vectors. The mean is $(1/4,1/4,1/4)'$. The sum of the outer products $X_iX_i'$ is the identity matrix, and $4\bar X\bar X'=\frac14\mathbf{1}\mathbf{1}'$, where $\mathbf{1}=(1,1,1)'$. Hence $$ S=I-\frac14\mathbf{1}\mathbf{1}'. $$ Along $\mathbf{1}$, this matrix has eigenvalue $1-3/4=1/4$; in the two-dimensional orthogonal complement, it has eigenvalue $1$. Thus $S$ has rank three, and $S/4$ is positive definite. This example shows that $n=p+1$ can provide full rank; the almost-sure condition will be $n>p$.

The Exact Existence Condition

For independent observations from $N_p(\mu,\Sigma)$, with $\mu\in\mathbb{R}^p$ and $\Sigma$ positive definite, the log-likelihood is $$ \ell(\mu,\Sigma) =-\frac{np}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma| -\frac12\sum_{i=1}^n (X_i-\mu)'\Sigma^{-1}(X_i-\mu). $$ The decomposition $$ \sum_{i=1}^n (X_i-\mu)(X_i-\mu)' =S+n(\bar X-\mu)(\bar X-\mu)' $$ shows that, for any fixed $\Sigma$, the likelihood is uniquely maximized over $\mu$ at $\bar X$. The remaining question is whether maximizing over positive-definite $\Sigma$ produces a valid covariance matrix.

Theorem: For any observed $X_1,\ldots,X_n$, the unrestricted Gaussian likelihood over $\mu\in\mathbb{R}^p$ and positive-definite $\Sigma$ has a finite maximum if and only if $S$ has rank $p$. When it does, the unique maximizers are $\hat\mu=\bar X$ and $\hat\Sigma=S/n$. If $S$ is singular, the likelihood is unbounded above.

Proof. First suppose $S$ is positive definite and set $A=S/n$. At $\mu=\bar X$, compare any positive-definite $\Sigma$ with $A$. Let $c_1,\ldots,c_p$ be the positive eigenvalues of $A^{-1/2}\Sigma A^{-1/2}$. Since $S=nA$, the difference in log-likelihoods is $$ \ell(\bar X,\Sigma)-\ell(\bar X,A) =-\frac{n}{2}\sum_{j=1}^p\left(\log c_j+\frac{1}{c_j}-1\right). $$ For $c>0$, $g(c)=\log c+1/c-1$ has derivative $g'(c)=(c-1)/c^2$. It decreases up to $1$, increases after $1$, and $g(1)=0$; hence $g(c)\geq 0$, with equality only at $c=1$. The displayed difference is therefore nonpositive, with equality only if all $c_j=1$, which implies $\Sigma=A$. For a general $\mu$, the likelihood has the additional log-likelihood difference $$ -\frac{n}{2}(\mu-\bar X)'\Sigma^{-1}(\mu-\bar X)\leq 0, $$ with equality only at $\mu=\bar X$. Thus $(\bar X,S/n)$ is the unique global maximizer.

Now suppose $S$ has rank $r<p$. Choose an orthonormal eigenbasis for $S$, with positive eigenvalues $\lambda_1,\ldots,\lambda_r$ and $p-r$ zero eigenvalues. Set $\mu=\bar X$, and choose $\Sigma_\varepsilon$ diagonal in this basis, with eigenvalues $\lambda_j/n$ for $1\leq j\leq r$ and $\varepsilon$ for each zero-eigenvalue direction. For every $\varepsilon>0$, this is an allowed covariance matrix. Its determinant is $$ |\Sigma_\varepsilon|=\left(\prod_{j=1}^r\frac{\lambda_j}{n}\right)\varepsilon^{p-r}, $$ and $\operatorname{tr}(\Sigma_\varepsilon^{-1}S)=nr$, which does not depend on $\varepsilon$. Consequently, the term $-\frac{n}{2}\log|\Sigma_\varepsilon|$ tends to $+\infty$ as $\varepsilon$ tends to zero, while the quadratic term remains fixed. The likelihood is unbounded above, so no finite maximum is attained in the positive-definite parameter space. This proves both directions. $\square$

When Does Full Rank Occur Almost Surely?

The scatter rank has a simple upper bound: centering the observations imposes the relation $\sum_i(X_i-\bar X)=0$, so $\operatorname{rank}(S)\leq\min(p,n-1)$. In particular, $n\leq p$ makes full rank impossible. Under the nonsingular Gaussian model, $n>p$ is also sufficient with probability one.

Theorem: Suppose $X_1,\ldots,X_n$ are independent $N_p(\mu,\Sigma)$ observations with positive-definite $\Sigma$. Then $$ P\bigl(\operatorname{rank}(S)=p\bigr)= \begin{cases} 0,&n\leq p,\\ 1,&n>p. \end{cases} $$ Thus the unrestricted covariance MLE exists with probability one exactly when $n>p$.

Proof. Let $D$ be the $p\times n$ matrix whose columns are $X_1,\ldots,X_n$. Choose an orthogonal $n\times n$ matrix $H$ whose first column is $\mathbf{1}/\sqrt n$. Then $$ S=D\left(I-\frac{\mathbf{1}\mathbf{1}'}{n}\right)D' =(DH)_{[:,2:n]}(DH)_{[:,2:n]}'. $$ The notation $(DH)_{[:,2:n]}$ denotes the matrix formed from columns 2 through $n$ of $DH$. By orthogonal invariance of the joint Gaussian distribution, these $n-1$ columns are independent $N_p(0,\Sigma)$ vectors. If $n-1<p$, their matrix has rank at most $n-1<p$, so $S$ cannot have rank $p$.

If $n-1\geq p$, consider the first $p$ of those $n-1$ columns. Their $p\times p$ matrix has a joint density with respect to Lebesgue measure, since $\Sigma$ is positive definite. Its determinant is a nonzero polynomial in its entries: it is nonzero, for example, at the identity matrix. The zero set of a nonzero polynomial has Lebesgue measure zero, so this determinant is nonzero with probability one. The columns therefore have rank $p$ with probability one, as does the full residual matrix, and $S$ has rank $p$. This establishes the stated probabilities and, by the first theorem, the almost-sure existence condition. $\square$

Common Pitfall and Practical Interpretation

Wrong reasoning: “There are at least as many observations as dimensions, so the covariance MLE exists.” When $n=p$, centering leaves at most $n-1=p-1$ independent residual directions. The scatter matrix must be singular, regardless of the observed values.

Correction: With an unknown mean, the necessary and sufficient almost-sure condition under a nonsingular Gaussian model is $n>p$, not $n\geq p$. The extra observation beyond the dimension accounts for the loss of one rank from centering. More generally, even when $n>p$, a particular data set with a lower-dimensional arrangement has singular scatter and no finite unrestricted covariance MLE.

This condition matters when fitting a full covariance matrix: if the sample size is too small relative to dimension, the problem is not merely that the estimate may be unstable. The likelihood itself has no finite maximizer over positive-definite covariance matrices. Restricted or regularized covariance models can address that difficulty, but they change the parameter space and are distinct estimation problems.

Check Your Understanding

Use the rank condition and likelihood argument to answer these questions.

Why can the centered scatter matrix have rank at most $n-1$, even when $p$ is much larger than one?
For a data set with $p=4$ and $n=4$, can the unrestricted covariance MLE exist? Explain.
If $S$ is singular, what happens to the Gaussian likelihood when covariance eigenvalues shrink in directions in the null space of $S$?
Under independent $N_p(\mu,\Sigma)$ sampling with positive-definite $\Sigma$, what is the probability that $S$ has full rank when $n=p+1$?
Why does the almost-sure rank result require a positive-definite population covariance matrix?