Estimation of covariance matrices

In multivariate statistics, the importance of the Wishart distribution stems in part from the fact that it is the probability distribution of the maximum likelihood estimator of the covariance matrix of a multivariate normal distribution. Although no one is surprised that the estimator of the population covariance matrix is simply the sample covariance matrix, the mathematical derivation is perhaps not widely known and is surprisingly subtle and elegant.

Contents

1 The multivariate normal distribution

2 Maximum-likelihood estimation

2.1 First steps
2.2 The trace of a 1 × 1 matrix
2.3 Using the spectral theorem
2.4 Concluding steps

The multivariate normal distribution

A random vector X ∈ R^p×1 (a p×1 "column vector") has a multivariate normal distribution with a nonsingular covariance matrix Σ precisely if Σ ∈ R^{p × p} is a positive-definite matrix and the probability density function of X is

$f(x)=[\mathrm{constant}]\cdot \det(\Sigma)^{-1/2} \exp\left({-1 \over 2} (x-\mu)^T \Sigma^{-1} (x-\mu)\right)$

where μ ∈ R^p×1 is the expected value. The matrix Σ is the higher-dimensional analog of what in one dimension would be the variance.

Maximum-likelihood estimation

Suppose now that X₁, ..., X_n are independent and identically distributed with the distribution above. Based on the observed values x₁, ..., x_n of this sample, we wish to estimate Σ (we adhere to the convention of writing random variables as capital letters and data as lower-case letters).

First steps

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

$\overline{x}=(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 of μ does not depend on Σ, we can just substitute it for μ in the likelihood function

$L(\mu,\Sigma)=[\mathrm{constant}]\cdot \prod_{i=1}^n \det(\Sigma)^{-1/2} \exp\left({-1 \over 2} (x_i-\mu)^T \Sigma^{-1} (x_i-\mu)\right)$

$\propto \det(\Sigma)^{-n/2} \exp\left({-1 \over 2} \sum_{i=1}^n (x_i-\mu)^T \Sigma^{-1} (x_i-\mu) \right)$

and then seek the value of Σ that maximizes this.

We have

$L(\overline{x},\Sigma) \propto \det(\Sigma)^{-n/2} \exp\left({-1 \over 2} \sum_{i=1}^n (x_i-\overline{x})^T \Sigma^{-1} (x_i-\overline{x})\right).$

The trace of a 1 × 1 matrix

Now we come to the first surprising step.

Regard the scalar $(x_i-\overline{x})^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

$\det(\Sigma)^{-n/2} \exp\left({-1 \over 2} \sum_{i=1}^n \operatorname{tr}((x_i-\mu)^T \Sigma^{-1} (x_i-\mu)) \right)$

$=\det(\Sigma)^{-n/2} \exp\left({-1 \over 2} \sum_{i=1}^n \operatorname{tr}((x_i-\mu) (x_i-\mu)^T \Sigma^{-1}) \right)$

(so now we are taking the trace of a p×p matrix!)

$=\det(\Sigma)^{-n/2} \exp\left({-1 \over 2} \operatorname{tr} \left( \sum_{i=1}^n (x_i-\mu) (x_i-\mu)^T \Sigma^{-1} \right) \right)$

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

where

$S=\sum_{i=1}^n (x_i-\overline{x}) (x_i-\overline{x})^T \in \mathbf{R}^{p\times p}.$

Using the spectral theorem

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

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

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

$\det(S)^{-n/2} \det(B)^{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

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

reduces to the problem of finding the values of the diagonal entries λ₁, ..., λ_p that maximize

$\lambda_i^{n/2} \exp(-\lambda_i/2).$

This is just a calculus problem and we get λ_i = n, so that B = n I_p, i.e., n times the p×p identity matrix.

Concluding steps

Finally we get

Σ = S 1 / 2 B - 1 S 1 / 2 = S 1 / 2 ((1 / n) I p) S 1 / 2 = S / n,

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

${S \over n} = {1 \over n}\sum_{i=1}^n (X_i-\overline{X})(X_i-\overline{X})^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. This random matrix can be shown to have a Wishart distribution with n − 1 degrees of freedom.