Stein's loss for multivariate normal covariance estimator 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}) \Vert N(0,\Sigma)) = \log|\hat{\Sigma}^{-1}\Sigma| - C\int_x [x^T (\hat{\Sigma}^{-1} - \Sigma^{-1})x ]\exp\left\{-\frac{1}{2}x^T \hat{\Sigma}^{-1} x\right\} dx$$
where
$$C = (2\pi)^{-k/2}|\hat{\Sigma}^{-1}|^{1/2}$$
however I am not sure how to simplify the integral.
 A: 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.
A: Wikipedia would seem to confirm your suspicion; Stein's loss is the KLD (up to a 1/2 multiplier):
$$D_{KL}(\mathcal{N}_0\Vert\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.
A: A simple way to compute the integral is to notice that it's equal to  $\mathbb{E}[\varepsilon^\top \Lambda \varepsilon]$ with $\varepsilon$ of mean $0$ and covariance $\hat \Sigma$, and $\Lambda = \hat \Sigma^{-1} - \Sigma^{-1}$.
Independently on the distribution of $\varepsilon$, this is equal to $\mathrm{tr} [\Lambda \,\mathrm{cov}[\varepsilon]]$
Indeed:
\begin{align*} 
  \mathbb{E}[\varepsilon^\top \Lambda \varepsilon] 
  &= \mathbb{E}[\mathrm{tr} \, \varepsilon^\top \Lambda \varepsilon] \\
  &= \mathbb{E}[\mathrm{tr} \,  \Lambda \varepsilon \varepsilon^\top] \\
  &= \mathrm{tr} \, \Lambda \mathrm{cov}[ \varepsilon] \\
  &= \mathrm{tr} \, [(\hat \Sigma^{-1} - \Sigma^{-1}) \hat \Sigma] \enspace,\\
\end{align*}
which gives the desired result.
