Maximum-Likelihood and Unbiased Covariance Estimates

Two Covariance Estimates, Two Denominators

In a multivariate regression, each case contributes a vector of responses, and the covariance matrix describes variation and association among the response variables after accounting for the mean model. Once the regression has been fitted, the residual sum-of-squares-and-cross-products matrix summarizes that remaining variation. There are two common ways to scale it: divide by the number of cases, or divide by the number of residual degrees of freedom. These choices answer different questions.

The first scaling arises from Gaussian maximum likelihood. The second corrects for the fact that fitting the regression mean uses some of the available information. To distinguish them precisely, consider a fixed-design multivariate linear model with \(n\) independent response vectors and \(q\) outcomes.

The projection formulation makes the definitions valid whenever the design has rank \(r\), regardless of how many columns are used to represent its column space. Under the model, \(M_XXB=0\), so \(\hat U=M_XU\). The matrix \(E\) is \(q\times q\); its diagonal entries are residual sums of squares for individual outcomes, and its off-diagonal entries are residual cross-products between outcomes.

Why the Gaussian MLE Uses \(n\)

For a positive-definite covariance matrix \(\Sigma\), the Gaussian log-likelihood, apart from constants that do not depend on the parameters, is

For each fixed \(\Sigma\), maximizing over the fitted mean is equivalent to minimizing the weighted residual term. The least-squares fitted matrix is \(P_XY\), and the corresponding residual SSCP is \(E\). After fitting the mean, the covariance-dependent part of the log-likelihood is therefore

Proof. Set \(A=E^{1/2}\Sigma^{-1}E^{1/2}\), which is positive definite. Because \(|\Sigma|=|E|/|A|\) and \(\operatorname{tr}(\Sigma^{-1}E)=\operatorname{tr}(A)\), maximizing the profiled log-likelihood over \(\Sigma\) is equivalent to maximizing

If the eigenvalues of \(A\) are \(\lambda_1,\ldots,\lambda_q>0\), this expression is the sum, over \(j\), of \(n\log(\lambda_j)/2-\lambda_j/2\). Each term has derivative \(n/(2\lambda_j)-1/2\), which is zero at \(\lambda_j=n\), and second derivative \(-n/(2\lambda_j^2)<0\). Thus each term is uniquely maximized at \(n\), so \(A=nI_q\). Rearranging gives \(\Sigma=E/n\). \(\square\)

This likelihood calculation explains the denominator: the Gaussian likelihood has \(n\) response vectors, and its covariance parameter is fitted using the residual SSCP from those \(n\) vectors. Maximum likelihood does not automatically correct for the estimation of the mean model. As the next result shows, that omission produces a downward bias.

The Bias Correction From Residual Degrees of Freedom

The projection \(M_X\) is symmetric and idempotent, and its rank and trace are both \(n-r=\nu\). These facts yield the expectation of \(E\) directly, without requiring a distributional approximation.

Proof. Since \(M_XXB=0\), \(E=U'M_XU\). Write \(u_i'\) for row \(i\) of \(U\), and \(m_{ij}\) for entry \((i,j)\) of \(M_X\). For response coordinates \(a\) and \(b\),

In the second equality, independence and zero means make \(E[U_{ia}U_{jb}]=0\) for \(i\ne j\); for \(i=j\), the expectation is \(\Sigma_{ab}\). Thus the matrix expectation is \(E[E]=\nu\Sigma\). Dividing by \(n\) gives \(E[E/n]=(\nu/n)\Sigma=(1-r/n)\Sigma\), so the bias is \(-r\Sigma/n\). Dividing by \(\nu\) gives \(E[E/\nu]=\Sigma\). \(\square\)

The factor \(\nu\) counts the dimensions of the observation space left after projecting out the fitted mean. It is not generally the number of columns in \(X\): redundant design columns do not increase its rank and do not use additional dimensions. The correction is based on \(\operatorname{rank}(X)\), not simply on a nominal column count.

Worked Examples

Interpreting the Denominators

Neither denominator is universally “the right one” without specifying the purpose. If the task is to report the Gaussian maximum-likelihood estimate, use \(n\). If an unbiased covariance estimate is required, use \(n-r\). The correction does not make every aspect of the estimate optimal, nor does unbiasedness mean that its realized value equals the true covariance. It means that, over repeated samples under the stated model and fixed design, its matrix expectation equals \(\Sigma\).

A frequent error is to argue that the likelihood uses \(n\) cases, so \(E/n\) must also be unbiased. The likelihood calculation establishes where the maximum occurs, not the expectation of that maximum across repeated samples. The expectation calculation instead shows that fitting the mean leaves only \(\nu\) residual dimensions, so \(E\) averages to \(\nu\Sigma\), not \(n\Sigma\). Another error is to subtract the number of outcomes \(q\) from \(n\). Here each case contributes a \(q\)-vector, but the regression mean removes \(r\) dimensions in the case space; the residual degrees of freedom are \(n-r\), not \(n-r-q\).

For a positive-definite Gaussian MLE, the residual SSCP must itself be positive definite. Under the continuous Gaussian model this occurs with probability one when \(\nu\ge q\). If \(\nu<q\), \(E\) must be singular, as in the third example, and the likelihood has no finite maximizer over positive-definite covariance matrices. The unbiased estimate \(E/\nu\) may still be defined, but it too will be singular in that situation.