Linear Covariance Components and Likelihood Equations

Covariance Components as Likelihood Directions

A covariance model can be structured without assigning a separate parameter to every covariance entry. Instead, it can be built from known symmetric matrices, each describing a particular pattern of variation. The unknown coefficients then control how much each pattern contributes. Likelihood estimation asks whether the observed data support more or less variation in each of those component directions.

We consider independent observations $Y_1,\ldots,Y_n$, each a $p$-dimensional vector with known mean $\mu$ and Gaussian distribution. Keeping the mean known isolates the covariance score; the same calculation applies to residual vectors when the mean has been specified.

Definition: A linear covariance-component model has covariance matrix $$ \Sigma(\theta)=\sum_{k=1}^{q}\theta_k A_k, $$ where the $A_k$ are known symmetric $p\times p$ component matrices and $\theta=(\theta_1,\ldots,\theta_q)^\mathsf{T}$ is an unknown parameter vector. The admissible parameter set consists of values for which $\Sigma(\theta)$ is positive definite. In variance-component models the $A_k$ are commonly positive semidefinite and the $\theta_k$ are nonnegative, but the linear representation itself does not require that particular parameterization.

Write $r_i=Y_i-\mu$, and suppose $Y_i\sim N_p(\mu,\Sigma(\theta))$ independently. Define the residual sum-of-products matrix $S=\sum_{i=1}^{n}r_i r_i^\mathsf{T}$. The likelihood is a function of $\theta$ through $\Sigma(\theta)$. A likelihood score equation sets a partial derivative of the log-likelihood to zero. Such equations are necessary at an interior local maximum, but they do not by themselves guarantee that a solution exists or is a global maximum.

Deriving the Gaussian Score

For a positive-definite covariance matrix, the density of one observation is proportional to $|\Sigma|^{-1/2}\exp\{-r_i^\mathsf{T}\Sigma^{-1}r_i/2\}$. Independence makes the log-likelihood a sum over observations. Using $\operatorname{tr}(Brr^\mathsf{T})=r^\mathsf{T}Br$, it can be written compactly as follows.

\ell(\theta) =-\frac{np}{2}\log(2\pi) -\frac{n}{2}\log|\Sigma(\theta)| -\frac{1}{2}\operatorname{tr}\!\left(\Sigma(\theta)^{-1}S\right).

To differentiate this expression, consider a perturbation $\Sigma+hA_k$, with $h$ small enough that the perturbed matrix remains positive definite. The two required directional derivatives are

Matrix derivative identities: For symmetric positive-definite $\Sigma$ and a fixed matrix $A$, $$ \left.\frac{d}{dh}\log|\Sigma+hA|\right|_{h=0} =\operatorname{tr}(\Sigma^{-1}A), \qquad \left.\frac{d}{dh}(\Sigma+hA)^{-1}\right|_{h=0} =-\Sigma^{-1}A\Sigma^{-1}. $$

For the first identity, factor $\Sigma+hA=\Sigma(I+h\Sigma^{-1}A)$. The determinant expansion $|I+hC|=1+h\operatorname{tr}(C)+O(h^2)$ gives the derivative of the log-determinant. For the second, differentiate $(\Sigma+hA)(\Sigma+hA)^{-1}=I$ at $h=0$. The resulting equation is $A\Sigma^{-1}+\Sigma\,\frac{d}{dh}(\Sigma+hA)^{-1}|_{0}=0$; multiplication on the left by $\Sigma^{-1}$ gives the stated formula.

Theorem (Gaussian covariance score equations): Suppose $\Sigma(\theta)=\sum_{k=1}^{q}\theta_k A_k$ is positive definite on an open parameter set, the $A_k$ are fixed symmetric matrices, and the observations are independent $N_p(\mu,\Sigma(\theta))$ with known $\mu$. Then $$ \frac{\partial\ell}{\partial\theta_k} =\frac{1}{2}\operatorname{tr}\!\left[ \left\{\Sigma^{-1}S\Sigma^{-1}-n\Sigma^{-1}\right\}A_k \right]. $$ Every interior stationary point therefore satisfies, for $k=1,\ldots,q$, $$ \operatorname{tr}\!\left[ \left\{\Sigma^{-1}S\Sigma^{-1}-n\Sigma^{-1}\right\}A_k \right]=0. $$

Proof. Since $\partial\Sigma/\partial\theta_k=A_k$, the log-determinant term contributes $-\frac{n}{2}\operatorname{tr}(\Sigma^{-1}A_k)$. By the inverse derivative identity, the derivative of $\operatorname{tr}(\Sigma^{-1}S)$ is $-\operatorname{tr}(\Sigma^{-1}A_k\Sigma^{-1}S)$. Its contribution to the log-likelihood is therefore $+\frac12\operatorname{tr}(\Sigma^{-1}A_k\Sigma^{-1}S)$. The cyclic property of trace rewrites this as $+\frac12\operatorname{tr}(\Sigma^{-1}S\Sigma^{-1}A_k)$. Adding the terms gives the score formula. At an interior stationary point each partial derivative is zero, giving the displayed equations. $\square$

The common matrix in braces measures a discrepancy between the observed residual sum of products and the covariance model, after weighting by $\Sigma^{-1}$. Each score tests that discrepancy in one particular component direction $A_k$. Thus the equations require trace orthogonality to the modeled directions; they do not generally require every entry of $S$ to equal the corresponding entry of $n\Sigma$.

Worked Examples

Worked Example: Independent Coordinates With Separate Variances

Suppose $p=2$, $\mu=0$, and the two coordinates have separate variances but zero covariance: $$ \Sigma(\theta)=\theta_1 \begin{pmatrix}1&0\\0&0\end{pmatrix} +\theta_2 \begin{pmatrix}0&0\\0&1\end{pmatrix} = \begin{pmatrix}\theta_1&0\\0&\theta_2\end{pmatrix}, \qquad \theta_1,\theta_2>0. $$ For $n=2$, take observations $Y_1=(1,2)^\mathsf{T}$ and $Y_2=(3,0)^\mathsf{T}$. The diagonal entries of $S$ are $1^2+3^2=10$ and $2^2+0^2=4$, so $S_{11}=10$ and $S_{22}=4$.

For the first component, the score equation reduces to $S_{11}/\theta_1^2-n/\theta_1=0$. Multiplying by $\theta_1^2$ gives $S_{11}=n\theta_1$, so $\hat\theta_1=10/2=5$. Similarly, $S_{22}=n\theta_2$, giving $\hat\theta_2=4/2=2$. Both estimates are positive, and substitution gives zero for each score: $$ \frac{1}{2}\left(\frac{10}{5^2}-\frac{2}{5}\right)=0, \qquad \frac{1}{2}\left(\frac{4}{2^2}-\frac{2}{2}\right)=0. $$

In this example the component matrices select separate diagonal entries, so the trace equations simplify to one variance equation per coordinate. That simple separation is not typical when components induce correlations.

Worked Example: Estimating Random-Intercept and Residual Variances

Consider repeated measurements at two occasions, with covariance $$ \Sigma(a,b)=a \begin{pmatrix}1&1\\1&1\end{pmatrix} +b \begin{pmatrix}1&0\\0&1\end{pmatrix}, \qquad a>0,\quad b>0. $$ Here $a$ is the random-intercept variance and $b$ is the residual variance. Suppose $\mu=0$, $n=2$, and the observations are $Y_1=(3,1)^\mathsf{T}$ and $Y_2=(1,3)^\mathsf{T}$.

Use the orthonormal directions $e_+=(1,1)^\mathsf{T}/\sqrt{2}$ and $e_-=(1,-1)^\mathsf{T}/\sqrt{2}$. In these directions the covariance eigenvalues are $\lambda_+=b+2a$ and $\lambda_-=b$. The transformed observations are $e_+^\mathsf{T}Y_1=2\sqrt{2}$, $e_+^\mathsf{T}Y_2=2\sqrt{2}$, $e_-^\mathsf{T}Y_1=\sqrt{2}$, and $e_-^\mathsf{T}Y_2=-\sqrt{2}$. Thus the sums of squared transformed observations are $Q_+=8+8=16$ and $Q_-=2+2=4$.

In the transformed coordinates the log-likelihood, up to a constant, is $$ \ell(a,b)=-\frac{n}{2}\log\lambda_+-\frac{Q_+}{2\lambda_+} -\frac{n}{2}\log\lambda_- -\frac{Q_-}{2\lambda_-}. $$ Differentiating with respect to an eigenvalue gives $-n/(2\lambda)+Q/(2\lambda^2)$. Because $\partial\lambda_+/\partial a=2$, $\partial\lambda_-/\partial a=0$, and both eigenvalues have derivative 1 with respect to $b$, the two score equations imply $$ \lambda_+=Q_+/n=16/2=8, \qquad \lambda_-=Q_-/n=4/2=2. $$ Recover the variance components from $b=\lambda_-$ and $a=(\lambda_+-\lambda_-)/2$: $$ \hat b=2,\qquad \hat a=(8-2)/2=3. $$ Both are positive, and they produce eigenvalues $2+2(3)=8$ and $2$, as required. The covariance estimate is therefore $\begin{pmatrix}5&3\\3&5\end{pmatrix}$.

Worked Example: Evaluating Scores in Correlated Component Directions

Let $n=1$, $p=2$, $\mu=0$, and use the components $A_1=I_2$ and $A_2=\begin{pmatrix}1&1\\1&1\end{pmatrix}$. At $\theta_1=\theta_2=1$, the covariance is $\Sigma=\begin{pmatrix}2&1\\1&2\end{pmatrix}$. For the observation $Y_1=(1,0)^\mathsf{T}$, $S=\begin{pmatrix}1&0\\0&0\end{pmatrix}$, while $$ \Sigma^{-1}=\frac{1}{3}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}, \qquad \Sigma^{-1}S\Sigma^{-1} =\begin{pmatrix}4/9&-2/9\\-2/9&1/9\end{pmatrix}. $$

The score discrepancy matrix is $$ B=\Sigma^{-1}S\Sigma^{-1}-\Sigma^{-1} =\begin{pmatrix}-2/9&1/9\\1/9&-5/9\end{pmatrix}. $$ For $A_1=I_2$, $\operatorname{tr}(BA_1)=\operatorname{tr}(B)=-7/9$, so $\partial\ell/\partial\theta_1=-7/18$. For $A_2$, $\operatorname{tr}(BA_2)$ is the sum of the four entries of $B$, which is $(-2+1+1-5)/9=-5/9$; hence $\partial\ell/\partial\theta_2=-5/18$. Neither score is zero, so this parameter value is not a stationary point. Each derivative reflects the whole component pattern, including its off-diagonal entries.

What the Equations Establish—and a Common Pitfall

The equations are useful because they turn a covariance likelihood into one scalar condition per unknown component coefficient. Their meaning can also be checked at the population level. Under the stated model, $E[S]=n\Sigma$, since each residual has mean zero and covariance $\Sigma$. Therefore $$ E[\Sigma^{-1}S\Sigma^{-1}-n\Sigma^{-1}] =\Sigma^{-1}(n\Sigma)\Sigma^{-1}-n\Sigma^{-1}=0. $$ Taking the trace against any fixed $A_k$ shows that each score has expectation zero when evaluated at the true covariance parameter. This is a property of the score under the model; it does not mean the observed score is zero at the true parameter for every finite sample.

Common mistake: matching covariance entries instead of solving component score equations. It is tempting to set $S/n=\Sigma$ and equate matrix entries. That is generally not the likelihood equation when the covariance is restricted to a structured family: $S/n$ may not belong to that family, and the likelihood weights discrepancies through $\Sigma^{-1}$. The correction is to compute $\Sigma^{-1}S\Sigma^{-1}-n\Sigma^{-1}$, then take its trace against each $A_k$. Only when the component structure makes those equations reduce to entrywise matching is that shortcut justified.

A second caution is that zero scores characterize interior stationary points, not automatically maxima. The parameter domain may include a boundary, the likelihood may have multiple stationary points, and a stationary point may be a saddle point. The equations derived here are the central calculus conditions; deciding which admissible candidate is the maximum requires checking the likelihood and the parameter domain.

Key takeaway: In a Gaussian linear covariance-component model, the score for component $A_k$ is one half the trace of the weighted covariance discrepancy against $A_k$. An interior likelihood solution makes each such componentwise trace equal zero.

Check Your Understanding

Use the score formula and the examples above to answer the following.

Why must the covariance matrix remain positive definite when differentiating the Gaussian log-likelihood?
For independent observations with known mean, what matrix summarizes the residual products in the covariance likelihood?
In the random-intercept example, why does the contrast direction $e_-$ have variance $b$ rather than $a+b$?
Why does a zero score in every component direction not, by itself, prove that a parameter value is a global maximum?
When is equating $S/n$ to the fitted covariance matrix a valid shortcut, and why does it fail for many structured covariance families?

Linear Covariance Components and Likelihood Equations

What You'll Learn