Stein's loss for multivariate normal covariance estimator

Question

In James and Stein (1961) the authors consider the following loss function for an estimator $\hat{\Sigma}$ of the covariance matrix $\Sigma$ of a multivariate normal distribution:

$$L(\hat{\Sigma}) = tr[\hat{\Sigma}\Sigma^{-1}] - \log|\hat{\Sigma}\Sigma^{-1}| - p.$$

This loss function has since been referred to as `Stein's loss' in several papers on regularized covariance estimation. Is there any intuitive justification for this loss function?

I noticed that the loss function resembles the KL divergence between two multivariate normal distributions with the same means and with covariances $\hat{\Sigma}$ and $\Sigma$:

$$2KL(N(0,\hat{\Sigma}) || N(0,\Sigma)) = \log|\hat{\Sigma}^{-1}\Sigma| - C\int_x [x^T (\hat{\Sigma}^{-1} - \Sigma^{-1})x ]\exp\{-\frac{1}{2}x^T \hat{\Sigma}^{-1} x\} dx$$

where

$$C = (2\pi)^{-k/2}|\hat{\Sigma}^{-1}|^{1/2}$$

however I am not sure how to simplify the integral.

Answer 1

Wikipedia would seem to confirm your suspicion; Stein's loss is the KLD (up to a 1/2 multiplier):

$$D_{KL}(\mathcal{N}_0||\mathcal{N}_1)=\frac{1}{2}\bigg(\text{tr}(\Sigma_1^{-1}\Sigma_0)+(\mu_1-\mu_0)^T \Sigma_1^{-1}(\mu_1-\mu_0) - \ln\bigg(\frac{\text{det}\Sigma_0}{\text{det}\Sigma_1}\bigg)-k\bigg)$$

I'll confess to not having done the math, but I'm pretty sure that you can arrange terms so it's equivalent to the expectation of a quadratic form of Gaussian rv's.

Answer 2

I worked out the details omitted from JMS's answer.

Let $P_0 = \Sigma_0^{-1}$ and $P_1 = \Sigma_1^{-1}$

We have

$$KL(N_0||N_1) = \int N_0(x) \frac{1}{2}[\ln |P_0| - \ln|P_1| + (x-\mu_0)^T P_0 (x-\mu_0) - (x-\mu_1)^T P_1 (x-\mu_1)] dx$$ $$ = \frac{1}{2}\ln|P_0\Sigma_1| - \frac{1}{2} \mathbb{E}[(x-\mu_0)^T P_0 (x-\mu_0) - (x-\mu_1)^T P_1 (x-\mu_1)] $$ $$ = \frac{1}{2}\ln|P_0\Sigma_1| + \frac{1}{2}(\mu_0-\mu_1)^T P_1 (\mu_0-\mu_1)- \frac{1}{2} \mathbb{E}[x^T P_0 x - x^T P_1 x - \mu_0^T (P_0 - P_1) \mu_0] $$ $$ = \frac{1}{2}\ln|P_0\Sigma_1| + \frac{1}{2}(\mu_0-\mu_1)^T P_1 (\mu_0-\mu_1)- \frac{1}{2} \mathbb{E}[(x-\mu_0)^T (P_0-P_1) (x-\mu_0)] $$

To simplify notation, assume $\mu_0 = 0$. It remains to show that $$\mathbb{E}[x^T P_0x-x^T P_1 x] = k - tr(P_1\Sigma_0)$$

But $x^T P_0 x = x^T \Sigma_0 x$ has a Chi-squared distribution and therefore has expectation $k$. Meanwhile,

$$\mathbb{E}[x^T P_1 x] = \mathbb{E}[tr(x^T P_1 x)] = \mathbb{E}[tr(P_1 xx^T )]= tr(\mathbb{E}(P_1 xx^T)) = tr(P_1 \Sigma_0)$$

where the second equality can be obtained by treating $x$ as a square matrix padded by zeroes.