Deriving The Gaussian Precision-Matrix Likelihood

From a Gaussian model to a precision-matrix likelihood

In a Gaussian graphical model, the precision matrix is the parameter whose entries describe conditional associations: its diagonal entries govern conditional variances, while its off-diagonal entries encode conditional dependence. Estimation therefore requires a likelihood written directly in terms of this matrix. The resulting objective has two parts: a log-determinant term that reflects the volume of the Gaussian distribution, and a trace term that measures fit to the observed covariance.

Let $X_1,\ldots,X_n$ be independent observations in $\mathbb{R}^p$ from $N_p(\mu,\Sigma)$, where $\Sigma$ is positive definite. Define the precision matrix $\Omega=\Sigma^{-1}$, which is also positive definite. For observations $x_1,\ldots,x_n$, write

\ell(\mu,\Omega) = -\frac{np}{2}\log(2\pi) +\frac{n}{2}\log\det(\Omega) -\frac{1}{2}\sum_{i=1}^n (x_i-\mu)^\mathsf{T}\Omega(x_i-\mu).

This is the ordinary Gaussian log-likelihood after substituting $\Sigma=\Omega^{-1}$: because $\det(\Sigma)^{-1/2}=\det(\Omega)^{1/2}$, the covariance determinant contributes $+\tfrac{n}{2}\log\det(\Omega)$. The parameter space is the set of symmetric positive-definite $p\times p$ matrices. The positive-definiteness restriction is essential; it ensures that the matrix defines a valid nonsingular Gaussian distribution.

How the sample covariance enters

If the mean $\mu$ is known, define the centered second-moment matrix

S_\mu=\frac{1}{n}\sum_{i=1}^n(x_i-\mu)(x_i-\mu)^\mathsf{T}.

If $\mu$ is unknown, its maximum-likelihood estimate is the sample mean $\bar{x}=n^{-1}\sum_{i=1}^n x_i$. In that case define the maximum-likelihood sample covariance

S=\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x})^\mathsf{T}.

The denominator here is $n$, not $n-1$, because this matrix appears in the likelihood maximized over the mean. The familiar unbiased sample covariance uses a different normalization and should not be substituted without adjusting the likelihood.

For any vector $y$, the quadratic form $y^\mathsf{T}\Omega y$ equals $\operatorname{tr}(\Omega yy^\mathsf{T})$. Thus summing these forms turns the data-dependent part of the likelihood into a trace. With a known mean, the result is

\ell(\mu,\Omega) = -\frac{np}{2}\log(2\pi) +\frac{n}{2} \left\{\log\det(\Omega)-\operatorname{tr}(S_\mu\Omega)\right\}.

With an unknown mean, profiling out $\mu$ gives the same expression with $S$ in place of $S_\mu$:

\ell_{\mathrm{prof}}(\Omega) = -\frac{np}{2}\log(2\pi) +\frac{n}{2} \left\{\log\det(\Omega)-\operatorname{tr}(S\Omega)\right\}.

Definition: The Gaussian precision-matrix log-likelihood, up to a constant that does not depend on $\Omega$, is $ \frac{n}{2}\{\log\det(\Omega)-\operatorname{tr}(S\Omega)\}$. Here $S$ is the centered second-moment matrix when the mean is estimated, or $S_\mu$ when the mean is known. The optimization is over symmetric positive-definite $\Omega$.

Worked Example: Computing the Unconstrained Estimate

Suppose the mean is estimated and the observed sample covariance, using denominator $n$, is $S=\begin{pmatrix}4&1\\1&2\end{pmatrix}$. Find the unconstrained maximum-likelihood estimate of $\Omega$, and evaluate the data-dependent part of the log-likelihood there.

First, $S$ is positive definite: its first leading principal minor is $4>0$, and its determinant is $4(2)-1^2=7>0$. Therefore it is invertible, with

S^{-1} = \frac{1}{7} \begin{pmatrix} 2&-1\\ -1&4 \end{pmatrix}.

The unconstrained estimate is $\widehat{\Omega}=S^{-1}$. Its determinant is $1/7$, so $\log\det(\widehat{\Omega})=-\log 7$. Also,

\operatorname{tr}(S\widehat{\Omega}) = \operatorname{tr}(SS^{-1}) = \operatorname{tr}(I_2) = 2.

Hence the log-likelihood after removing the constant $-np\log(2\pi)/2$ is $ \frac{n}{2}(-\log 7-2)$. The calculation checks the scale as well as the estimate: the trace at the optimum is the dimension $p=2$.

The matrix score

The score is obtained by differentiating the log-likelihood with respect to the symmetric matrix $\Omega$. For a symmetric perturbation $H$, the differential of the log-determinant is $d\log\det(\Omega)[H]=\operatorname{tr}(\Omega^{-1}H)$. This identity follows by writing $\det(\Omega+tH)=\det(\Omega)\det(I+t\Omega^{-1}H)$ and using $\det(I+tA)=1+t\operatorname{tr}(A)+o(t)$ as $t$ tends to zero. The differential of $\operatorname{tr}(S\Omega)$ in direction $H$ is $\operatorname{tr}(SH)$.

Result: For the profiled Gaussian log-likelihood, the differential in any symmetric direction $H$ is $$ d\ell_{\mathrm{prof}}(\Omega)[H] = \frac{n}{2}\operatorname{tr}\{(\Omega^{-1}-S)H\}. $$ Under the Frobenius inner product $\langle A,B\rangle=\operatorname{tr}(A^\mathsf{T}B)$, the symmetric matrix score is therefore $ \frac{n}{2}(\Omega^{-1}-S)$.

This is a matrix gradient, not a collection of derivatives that ignores the symmetry constraint. In particular, if $\Omega_{jk}=\Omega_{kj}$ is treated as one free off-diagonal parameter for $j\ne k$, changing it changes both matrix entries. The coordinate derivative is then $n(\Omega^{-1}-S)_{jk}$, twice the corresponding off-diagonal entry of the matrix gradient. A diagonal coordinate derivative is $ \frac{n}{2}(\Omega^{-1}-S)_{jj}$.

Worked Example: Evaluating the Score at a Candidate Matrix

Let $n=10$, $S=\begin{pmatrix}2&1/2\\1/2&1\end{pmatrix}$, and consider the candidate $\Omega=\begin{pmatrix}1&1/5\\1/5&2\end{pmatrix}$. Its determinant is $2-1/25=49/25>0$, and its first diagonal entry is positive, so it is positive definite. Calculate the score.

The inverse of the candidate is

\Omega^{-1} = \frac{1}{49/25} \begin{pmatrix} 2&-1/5\\ -1/5&1 \end{pmatrix} = \frac{1}{49} \begin{pmatrix} 50&-5\\ -5&25 \end{pmatrix}.

Subtracting $S$ gives

\Omega^{-1}-S = \begin{pmatrix} -48/49&-59/98\\ -59/98&-24/49 \end{pmatrix}.

Thus the symmetric matrix score is $5\begin{pmatrix}-48/49&-59/98\\-59/98&-24/49\end{pmatrix}$. For the single off-diagonal coordinate $\Omega_{12}=\Omega_{21}$, its ordinary coordinate derivative is $10(-59/98)=-295/49$, not $5(-59/98)$. The factor of two comes from varying both symmetric entries.

Why the unconstrained maximum is the inverse covariance

When $S$ is positive definite, setting the score to zero gives $\Omega^{-1}=S$, suggesting $\widehat{\Omega}=S^{-1}$. The stationary equation alone does not establish a maximum, however. The following argument verifies that this candidate is the unique global maximizer.

Theorem: If $S$ is positive definite, the function $g(\Omega)=\log\det(\Omega)-\operatorname{tr}(S\Omega)$, over symmetric positive-definite $\Omega$, has the unique maximizer $\widehat{\Omega}=S^{-1}$.

Transform the parameter.
Let $A=S^{1/2}\Omega S^{1/2}$. Since $S$ and $\Omega$ are positive definite, $A$ is positive definite. Conversely, every positive-definite $A$ corresponds to $\Omega=S^{-1/2}AS^{-1/2}$.

Rewrite both terms.
Determinant multiplication gives $\log\det(\Omega)=\log\det(A)-\log\det(S)$. Cyclic invariance of trace gives $\operatorname{tr}(S\Omega)=\operatorname{tr}(A)$.

Optimize over eigenvalues.
If $\lambda_1,\ldots,\lambda_p>0$ are the eigenvalues of $A$, then $g(\Omega)=-\log\det(S)+\sum_{j=1}^p(\log\lambda_j-\lambda_j)$. For every $\lambda>0$, $\log\lambda-\lambda\le -1$, with equality only at $\lambda=1$: its derivative is $1/\lambda-1$, positive below one and negative above one.

Identify equality.
The sum reaches its maximum only when every eigenvalue equals one. A positive-definite symmetric matrix with all eigenvalues equal to one is $A=I_p$, which corresponds to $\Omega=S^{-1}$. This proves the unique global maximum. $\square$

The determinant term rewards larger precision volume through $\log\det(\Omega)$, while the trace term penalizes precision matrices that assign a large quadratic cost to the observed covariance. At the optimum the two terms balance: $S\Omega=I_p$, and the fitted trace is $p$. This balance is the likelihood version of matching the fitted covariance to the observed second moments.

Worked Example: A Singular Sample Covariance

Suppose $p=2$, $n=2$, and the observations are $x_1=(0,0)^\mathsf{T}$ and $x_2=(1,0)^\mathsf{T}$, with the mean estimated. Then $\bar{x}=(1/2,0)^\mathsf{T}$, and $S=\begin{pmatrix}1/4&0\\0&0\end{pmatrix}$. Consider the valid positive-definite matrices $\Omega_t=\operatorname{diag}(4,t)$, where $t>0$.

For every $t>0$, $\log\det(\Omega_t)=\log(4t)$, while $\operatorname{tr}(S\Omega_t)=1$. Consequently, the part of the log-likelihood depending on $t$ is

\frac{n}{2}\{\log(4t)-1\} = \log(4t)-1.

This tends to $+\infty$ as $t$ tends to infinity. Thus the likelihood has no finite maximum over positive-definite precision matrices in this example. The zero eigenvalue of $S$ leaves a direction in which the trace penalty does not grow, while the log-determinant can grow without bound.

Common mistakes and interpretation

Reversing the determinant sign. In precision form the term is $+\log\det(\Omega)$, not $-\log\det(\Omega)$. The covariance-form term is $-\log\det(\Sigma)$, and substituting $\Sigma=\Omega^{-1}$ changes its sign.
Using the wrong covariance normalization. The likelihood expression uses the centered sum divided by $n$, including when the mean is estimated. Using the unbiased denominator $n-1$ changes the trace term.
Ignoring symmetry in an off-diagonal derivative. A single free off-diagonal coordinate changes two entries of $\Omega$. The coordinate derivative is twice the off-diagonal entry of the symmetric matrix gradient.
Assuming the score equation always gives a finite estimate. If $S$ is singular, $S^{-1}$ does not exist, and a stationary equation does not provide an interior solution. The likelihood can be unbounded, as the example shows.

These calculations matter because the likelihood separates the empirical information in the data, summarized by $S$, from the parameter, $\Omega$. This form makes the unconstrained estimation problem transparent and supplies the score used to assess candidate precision matrices. In high dimensions, the existence issue is not merely technical: a singular empirical covariance can prevent ordinary likelihood maximization from producing a finite estimate.

Check Your Understanding

Use the definitions and results above to answer these questions.

Why does the precision-form likelihood contain $+\log\det(\Omega)$, rather than a negative log-determinant?
For an unknown mean, why does the likelihood use a covariance matrix with denominator $n$?
Write the differential of the profiled log-likelihood in a symmetric direction $H$, and state its matrix score.
When $S$ is positive definite, what transformation reduces the global maximization to separate scalar inequalities?
How can the likelihood increase without bound when the sample covariance has a zero eigenvalue?