Deriving Gaussian Maximum-Likelihood Estimates

From the Gaussian Model to Maximum Likelihood

In multivariate regression, each case contributes a vector of responses, and the responses within a case may be correlated. The Gaussian likelihood accounts for that correlation through a common covariance matrix. Deriving its maximum clarifies two important facts: the coefficient estimate is the ordinary least-squares estimate, while the covariance estimate uses the residual sum-of-squares-and-cross-products matrix divided by the number of cases.

Let \(Y\) be an \(n\times p\) response matrix, \(X\) an \(n\times q\) design matrix of full column rank, \(B\) a \(q\times p\) coefficient matrix, and \(\Sigma\) a \(p\times p\) positive-definite covariance matrix. Assume the rows of \(Y\) are independent and that row \(i\) has a \(p\)-variate normal distribution with mean given by row \(i\) of \(XB\) and covariance \(\Sigma\). Equivalently, the errors have a matrix-normal distribution with independent rows and column covariance \(\Sigma\).

Assume throughout that \(X\) has rank \(q\). Define the ordinary least-squares coefficient matrix and its residual sum-of-squares-and-cross-products matrix by

The likelihood analysis will show that \(\hat B\) is also the maximum-likelihood estimate. A finite, positive-definite covariance estimate exists when \(S\) is positive definite. This requires \(n-q\ge p\), since the rank of \(\hat R\) cannot exceed \(n-q\); that dimension condition alone is not sufficient, so positive definiteness of the observed \(S\) must also be checked.

Why the Coefficient Estimate Is Least Squares

Fix any positive-definite \(\Sigma\). Maximizing the likelihood over \(B\) is equivalent to minimizing \(\operatorname{tr}\{\Sigma^{-1}(Y-XB)'(Y-XB)\}\). To see where that criterion is minimized, write \(B=\hat B+D\), where \(D\) is any \(q\times p\) matrix. The least-squares residual satisfies \(X'\hat R=0\), so

The cross-products vanish because \(\hat R'X=0\) and \(X'\hat R=0\). Thus the change in the trace criterion, relative to \(B=\hat B\), is

If \(D\ne 0\), then \(X'X\) is positive definite and \(\Sigma^{-1}\) is positive definite, making this trace strictly positive. For example, it equals the squared sum of the entries of \((X'X)^{1/2}D\Sigma^{-1/2}\), which is nonzero when \(D\ne0\). Therefore \(\hat B\) is the unique likelihood maximizer over \(B\) for each fixed \(\Sigma\). Importantly, this argument does not require uncorrelated responses: it holds for every positive-definite \(\Sigma\).

The same result can be read directly from the likelihood: after substituting \(\hat B\), the residual matrix in the likelihood is \(\hat R\), so its quadratic term becomes \(\operatorname{tr}(\Sigma^{-1}S)\). The remaining task is to maximize over \(\Sigma\).

Maximizing the Likelihood over the Covariance

Substitute \(B=\hat B\) into the log-likelihood. Terms that do not involve \(\Sigma\) can be ignored for maximization, leaving

Suppose now that \(S\) is positive definite. Define \(C=S^{-1/2}\Sigma S^{-1/2}\), which is positive definite, so \(\Sigma=S^{1/2}CS^{1/2}\). This change of parameterization gives

The second identity follows from \(\Sigma^{-1}=S^{-1/2}C^{-1}S^{-1/2}\) and the cyclic property of the trace. Thus, apart from the constant \(n\log|S|\), maximizing the likelihood is equivalent to minimizing the objective \(n\log|C|+\operatorname{tr}(C^{-1})\).

To prove this, let \(c_1,\ldots,c_p\) be the positive eigenvalues of \(C\). The determinant is the product of these eigenvalues, and the trace of \(C^{-1}\) is the sum of their reciprocals. The objective is therefore

For a single positive eigenvalue, define \(g(c)=n\log c+1/c\). Its derivative is

This derivative is negative for \(0<c<1/n\) and positive for \(c>1/n\). Also, \(g(c)\) tends to infinity as \(c\) decreases to zero or increases without bound. It follows that \(g\) has its unique global minimum at \(c=1/n\). Each summand is minimized there, so the full objective has its unique minimum when all eigenvalues of \(C\) equal \(1/n\). A symmetric matrix with all eigenvalues equal to \(1/n\) is \(C=I_p/n\). Transforming back gives \(\Sigma=S^{1/2}(I_p/n)S^{1/2}=S/n\). This proves the result. \(\square\)

When Does a Finite Covariance Estimate Exist?

The positive-definiteness condition on \(S\) is not a technical convenience: it determines whether the likelihood has a finite maximizer in the positive-definite covariance parameter space. If \(S\) is singular, choose an orthonormal eigenbasis in which \(S\) has at least one zero eigenvalue. For illustration, let the eigenvalues be \(s_1,\ldots,s_r>0\), followed by \(p-r\) zeros. In this basis, consider covariance matrices with eigenvalues \(s_j/n\) in the positive-eigenvalue directions and \(\varepsilon\) in each zero-eigenvalue direction, where \(\varepsilon>0\). Then \(\operatorname{tr}(\Sigma^{-1}S)=rn\), which stays fixed, while \(|\Sigma|\) tends to zero as \(\varepsilon\) tends to zero. Since the likelihood includes the factor \(|\Sigma|^{-n/2}\), the likelihood tends to infinity. Thus no finite positive-definite covariance MLE exists when \(S\) is singular.

Common Pitfall and Practical Interpretation

A common mistake is to use the residual degrees of freedom in the covariance likelihood estimate and report \(S/(n-q)\). That divisor is associated with an unbiased covariance estimate under the model, not with the likelihood maximization performed here. The likelihood contains \(n\log|\Sigma|\), because the determinant contribution appears once for each of the \(n\) independent rows. Its optimization gives \(S/n\). The distinction between the maximum-likelihood and unbiased covariance estimates is important in applications; the likelihood derivation alone does not make the two estimates interchangeable.

Another tempting error is to omit the factor \(n\) while minimizing the transformed objective. The actual objective is \(n\log|C|+\operatorname{tr}(C^{-1})\), not \(\log|C|+\operatorname{tr}(C^{-1})\). Omitting \(n\) changes the eigenvalue minimizer and contradicts the likelihood’s determinant factor. Keeping the likelihood objective intact yields \(c_j=1/n\) for every eigenvalue and hence \(\hat\Sigma_{\mathrm{ML}}=S/n\).