Estimation of covariance matrices

The likelihood function is:

${\displaystyle {\mathcal {L}}(\mu ,\Sigma )=(2\pi )^{-{\frac {np}{2}}}\,\prod _{i=1}^{n}\det(\Sigma )^{-{\frac {1}{2}}}\exp \left(-{\frac {1}{2}}(x_{i}-\mu )^{\mathrm {T} }\Sigma ^{-1}(x_{i}-\mu )\right)}$

It is fairly readily shown that the maximum-likelihood estimate of the mean vector μ is the "sample mean" vector:

${\displaystyle {\overline {x}}={\frac {x_{1}+\cdots +x_{n}}{n}}.}$

See the section on estimation in the article on the normal distribution for details; the process here is similar.

Since the estimate ${\displaystyle {\bar {x}}}$ does not depend on Σ, we can just substitute it for μ in the likelihood function, getting

${\displaystyle {\mathcal {L}}({\overline {x}},\Sigma )\propto \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}(x_{i}-{\overline {x}})^{\mathrm {T} }\Sigma ^{-1}(x_{i}-{\overline {x}})\right),}$

and then seek the value of Σ that maximizes the likelihood of the data (in practice it is easier to work with log ${\displaystyle {\mathcal {L}}}$).

Now we come to the first surprising step: regard the scalar ${\displaystyle (x_{i}-{\overline {x}})^{\mathrm {T} }\Sigma ^{-1}(x_{i}-{\overline {x}})}$ as the trace of a 1×1 matrix. This makes it possible to use the identity tr(AB) = tr(BA) whenever A and B are matrices so shaped that both products exist. We get

${\displaystyle {\begin{aligned}{\mathcal {L}}({\overline {x}},\Sigma )&\propto \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}\left(\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\left(x_{i}-{\overline {x}}\right)\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}\operatorname {tr} \left(\left(x_{i}-{\overline {x}}\right)\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(\sum _{i=1}^{n}\left(x_{i}-{\overline {x}}\right)\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(S\Sigma ^{-1}\right)\right)\end{aligned}}}$

where

${\displaystyle S=\sum _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} }\in \mathbf {R} ^{p\times p}.}$

${\displaystyle S}$ is sometimes called the scatter matrix, and is positive definite if there exists a subset of the data consisting of ${\displaystyle p}$ affinely independent observations (which we will assume).

It follows from the spectral theorem of linear algebra that a positive-definite symmetric matrix S has a unique positive-definite symmetric square root S^1/2. We can again use the "cyclic property" of the trace to write

${\displaystyle \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(S^{\frac {1}{2}}\Sigma ^{-1}S^{\frac {1}{2}}\right)\right).}$

Let B = S^1/2 Σ ⁻¹ S^1/2. Then the expression above becomes

${\displaystyle \det(S)^{-{\frac {n}{2}}}\det(B)^{\frac {n}{2}}\exp \left(-{1 \over 2}\operatorname {tr} (B)\right).}$

The positive-definite matrix B can be diagonalized, and then the problem of finding the value of B that maximizes

${\displaystyle \det(B)^{\frac {n}{2}}\exp \left(-{1 \over 2}\operatorname {tr} (B)\right)}$

Since the trace of a square matrix equals the sum of eigenvalues ("trace and eigenvalues"), the equation reduces to the problem of finding the eigenvalues λ₁, ..., λ_p that maximize

${\displaystyle \lambda _{i}^{\frac {n}{2}}\exp \left(-{\frac {\lambda _{i}}{2}}\right).}$

This is just a calculus problem and we get λ_i = n for all i. Thus, assume Q is the matrix of eigen vectors, then

${\displaystyle B=Q(nI_{p})Q^{-1}=nI_{p}}$

i.e., n times the p×p identity matrix.

Finally we get

${\displaystyle \Sigma =S^{\frac {1}{2}}B^{-1}S^{\frac {1}{2}}=S^{\frac {1}{2}}\left({\frac {1}{n}}I_{p}\right)S^{\frac {1}{2}}={\frac {S}{n}},}$

i.e., the p×p "sample covariance matrix"

${\displaystyle {S \over n}={1 \over n}\sum _{i=1}^{n}(X_{i}-{\overline {X}})(X_{i}-{\overline {X}})^{\mathrm {T} }}$

is the maximum-likelihood estimator of the "population covariance matrix" Σ. At this point we are using a capital X rather than a lower-case x because we are thinking of it "as an estimator rather than as an estimate", i.e., as something random whose probability distribution we could profit by knowing. The random matrix S can be shown to have a Wishart distribution with n − 1 degrees of freedom.[5] That is:

${\displaystyle \sum _{i=1}^{n}(X_{i}-{\overline {X}})(X_{i}-{\overline {X}})^{\mathrm {T} }\sim W_{p}(\Sigma ,n-1).}$

An alternative derivation of the maximum likelihood estimator can be performed via matrix calculus formulae (see also differential of a determinant and differential of the inverse matrix). It also verifies the aforementioned fact about the maximum likelihood estimate of the mean. Re-write the likelihood in the log form using the trace trick:

${\displaystyle \ln {\mathcal {L}}(\mu ,\Sigma )=\operatorname {const} -{n \over 2}\ln \det(\Sigma )-{1 \over 2}\operatorname {tr} \left[\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }\right].}$

The differential of this log-likelihood is

${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=-{\frac {n}{2}}\operatorname {tr} \left[\Sigma ^{-1}\left\{d\Sigma \right\}\right]-{1 \over 2}\operatorname {tr} \left[-\Sigma ^{-1}\{d\Sigma \}\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }-2\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )\{d\mu \}^{\mathrm {T} }\right].}$

It naturally breaks down into the part related to the estimation of the mean, and to the part related to the estimation of the variance. The first order condition for maximum, ${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=0}$, is satisfied when the terms multiplying ${\displaystyle d\mu }$ and ${\displaystyle d\Sigma }$ are identically zero. Assuming (the maximum likelihood estimate of) ${\displaystyle \Sigma }$ is non-singular, the first order condition for the estimate of the mean vector is

${\displaystyle \sum _{i=1}^{n}(x_{i}-\mu )=0,}$

which leads to the maximum likelihood estimator

${\displaystyle {\widehat {\mu }}={\bar {X}}={1 \over n}\sum _{i=1}^{n}X_{i}.}$

This lets us simplify

${\displaystyle \sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }=\sum _{i=1}^{n}(x_{i}-{\bar {x}})(x_{i}-{\bar {x}})^{\mathrm {T} }=S}$

as defined above. Then the terms involving ${\displaystyle d\Sigma }$ in ${\displaystyle d\ln L}$ can be combined as

${\displaystyle -{1 \over 2}\operatorname {tr} \left(\Sigma ^{-1}\left\{d\Sigma \right\}\left[nI_{p}-\Sigma ^{-1}S\right]\right).}$

The first order condition ${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=0}$ will hold when the term in the square bracket is (matrix-valued) zero. Pre-multiplying the latter by ${\displaystyle \Sigma }$ and dividing by ${\displaystyle n}$ gives

${\displaystyle {\widehat {\Sigma }}={1 \over n}S,}$

which of course coincides with the canonical derivation given earlier.

Dwyer[6] points out that decomposition into two terms such as appears above is "unnecessary" and derives the estimator in two lines of working. Note that it may be not trivial to show that such derived estimator is the unique global maximizer for likelihood function.

Given a sample of n independent observations x₁,..., x_n of a p-dimensional zero-mean Gaussian random variable X with covariance R, the maximum likelihood estimator of R is given by

${\displaystyle {\hat {\mathbf {R} }}={1 \over n}\sum _{i=1}^{n}x_{i}x_{i}^{\mathrm {T} }.}$

The parameter ${\displaystyle R}$ belongs to the set of positive-definite matrices, which is a Riemannian manifold, not a vector space, hence the usual vector-space notions of expectation, i.e. "${\displaystyle \mathrm {E} [{\hat {\mathbf {R} }}]}$", and estimator bias must be generalized to manifolds to make sense of the problem of covariance matrix estimation. This can be done by defining the expectation of a manifold-valued estimator ${\displaystyle {\hat {\mathbf {R} }}}$ with respect to the manifold-valued point ${\displaystyle R}$ as

${\displaystyle \mathrm {E} _{\mathbf {R} }[{\hat {\mathbf {R} }}]\ {\stackrel {\mathrm {def} }{=}}\ \exp _{\mathbf {R} }\mathrm {E} \left[\exp _{\mathbf {R} }^{-1}{\hat {\mathbf {R} }}\right]}$

where

${\displaystyle \exp _{\mathbf {R} }({\hat {\mathbf {R} }})=\mathbf {R} ^{\frac {1}{2}}\exp \left(\mathbf {R} ^{-{\frac {1}{2}}}{\hat {\mathbf {R} }}\mathbf {R} ^{-{\frac {1}{2}}}\right)\mathbf {R} ^{\frac {1}{2}}}$

${\displaystyle \exp _{\mathbf {R} }^{-1}({\hat {\mathbf {R} }})=\mathbf {R} ^{\frac {1}{2}}\left(\log \mathbf {R} ^{-{\frac {1}{2}}}{\hat {\mathbf {R} }}\mathbf {R} ^{-{\frac {1}{2}}}\right)\mathbf {R} ^{\frac {1}{2}}}$

are the exponential map and inverse exponential map, respectively, "exp" and "log" denote the ordinary matrix exponential and matrix logarithm, and E[·] is the ordinary expectation operator defined on a vector space, in this case the tangent space of the manifold.[1]

The intrinsic bias vector field of the SCM estimator ${\displaystyle {\hat {\mathbf {R} }}}$ is defined to be

${\displaystyle \mathbf {B} ({\hat {\mathbf {R} }})=\exp _{\mathbf {R} }^{-1}\mathrm {E} _{\mathbf {R} }\left[{\hat {\mathbf {R} }}\right]=\mathrm {E} \left[\exp _{\mathbf {R} }^{-1}{\hat {\mathbf {R} }}\right]}$

The intrinsic estimator bias is then given by ${\displaystyle \exp _{\mathbf {R} }\mathbf {B} ({\hat {\mathbf {R} }})}$.

For complex Gaussian random variables, this bias vector field can be shown[1] to equal

${\displaystyle \mathbf {B} ({\hat {\mathbf {R} }})=-\beta (p,n)\mathbf {R} }$

where

${\displaystyle \beta (p,n)={\frac {1}{p}}\left(p\log n+p-\psi (n-p+1)+(n-p+1)\psi (n-p+2)+\psi (n+1)-(n+1)\psi (n+2)\right)}$

and ψ(·) is the digamma function. The intrinsic bias of the sample covariance matrix equals

${\displaystyle \exp _{\mathbf {R} }\mathbf {B} ({\hat {\mathbf {R} }})=e^{-\beta (p,n)}\mathbf {R} }$

and the SCM is asymptotically unbiased as n → ∞.

Similarly, the intrinsic inefficiency of the sample covariance matrix depends upon the Riemannian curvature of the space of positive-definite matrices.