Expectation of a function of the sample covariance matrix

Let $$X \sim N(\mu, \Sigma)$$ and let $$\hat\Sigma$$ be the sample covariance matrix calculated from an iid sample $$X_1, \dots, X_N$$. Also, let $$\phi_i$$ and $$\phi_j$$ be the orthonormal eigenvectors of $$\Sigma$$ associated with eigenvalues $$\lambda_i$$ and $$\lambda_j$$ respectively, where $$i \neq j$$.

A Quadratic Discriminant Function Based on Bias Rectification of Eigenvalues (Sakai et al. 2000) equation 12 states this result, without proof:

$$\mathbb E\left[\left(\phi_i^\top\hat\Sigma\phi_j\right)^2\right] = \frac{\lambda_i\lambda_j}{N-1}.$$

How should the proof be approached? Any explanations of the approach or sources are appreciated.

• @User1865345 thanks for adding the link to the paper and fixing the LaTeX! It's the little things that count... Oct 18, 2023 at 19:31

This is a quite challenging and interesting problem as it calls for many classical results of multivariate Gaussian distribution and technical matrix operations. To begin with, note that since $$\hat{\Sigma}$$ is invariant under translation, we can assume $$\mu = 0$$ in the calculation.

Suppose $$X$$ is a $$m$$-variate Gaussian vector. Introduce the following notations: \begin{align*} Z = \begin{bmatrix} X_1^T \\ \vdots \\ X_N^T \end{bmatrix} \in \mathbb{R}^{N \times m}, \; e = \begin{bmatrix} 1 \\ \vdots \\ 1 \end{bmatrix} \in \mathbb{R}^{N \times 1}, \; \bar{X} = N^{-1}\sum_{k = 1}^NX_k \in \mathbb{R}^{m \times 1}. \end{align*} By definition, \begin{align*} (N - 1)\hat{\Sigma} = \sum_{k = 1}^N(X_k - \bar{X})(X_k - \bar{X})^T = Z^T(I_{(N)} - N^{-1}ee^T)Z := Z^TPZ \in \mathbb{R}^{m \times m}. \end{align*} Therefore \begin{align*} \phi_i^T\hat{\Sigma}\phi_j = \frac{1}{N - 1}(Z\phi_i)^TP(Z\phi_j). \tag{1}\label{1} \end{align*}

My strategy is trying to leverage the expectation and variance of Gaussian quadratic forms. To this end, we need to rewrite $$\eqref{1}$$ in terms of Kronecker products and vectorization of $$Z^T$$ so that $$\phi_i, \phi_j$$ can be absorbed into the central non-random matrix. In detail, denote $$\operatorname{vec}(Z^T) \in \mathbb{R}^{N^2 \times 1}$$ by $$\xi$$, then $$\xi \sim N_{N^2}(0, I_{(N)} \otimes \Sigma)$$. It can be verified that \begin{align*} G &:= (N - 1)\phi_i^T\hat{\Sigma}\phi_j \\ &= \phi_i^TZ^TP(\phi_j^TZ^T)^T \\ &= (\operatorname{vec}(\phi_i^TZ^T))^TP\operatorname{vec}(\phi_j^TZ^T) \\ &= ((I_{(N)} \otimes \phi_i^T)\operatorname{vec}(Z^T))^TP (I_{(N)} \otimes \phi_j^T)\operatorname{vec}(Z^T) \\ &= \xi^T (I_{(N)} \otimes \phi_i)P(I_{(N)} \otimes \phi_j^T)\xi \\ &= \xi^T (P \otimes \phi_i)(I_{(N)} \otimes \phi_j^T)\xi \\ &= \xi^T (P\otimes (\phi_i\phi_j^T))\xi. \tag{2}\label{2} \end{align*} In the derivation above, the third equality used the fact that $$v = \operatorname{vec}(v)^T$$ for a row vector $$v$$. The fourth equality used the identity $$\operatorname{vec}(ABC) = (C^T \otimes A)\operatorname{vec}(B)$$. The fifth equality used $$(A \otimes B)^T = A^T \otimes B^T$$ and the definition $$\xi = \operatorname{vec}(Z^T)$$. The sixth and seventh equality used the identity $$(A \otimes B)(C \otimes D) = (AC)\otimes(BD)$$ of Kronecker product.

It then follows by the expectation of Gaussian quadratic form that \begin{align*} E[G] &= \operatorname{tr}((P\otimes (\phi_i\phi_j^T))(I_{(N)}\otimes\Sigma)) \\ &= \operatorname{tr}(P \otimes (\phi_i\phi_j^T\Sigma)) \\ &= \operatorname{tr}(P \otimes (\phi_i\lambda_j\phi_j^T)) \\ &= \operatorname{tr}(P)\operatorname{tr}(\lambda_j\phi_i\phi_j^T) \\ &= (N - 1)\lambda_j\phi_j^T\phi_i = 0. \tag{3}\label{3} \end{align*} In the derivation above, the second equality used the condition $$\Sigma\phi_j = \lambda_j\phi_j$$. The fourth equality used the property $$\operatorname{tr}(A \otimes B) = \operatorname{tr}(A)\operatorname{tr}(B)$$ of Kronecker product and trace. The fifth equality used the property $$\operatorname{tr}(AB) = \operatorname{tr}(BA)$$ of trace and the fact $$\operatorname{tr}(P) = N - 1$$. In the last equality, we used the orthogonality between $$\phi_i$$ and $$\phi_j$$. The same properties will be used repeatedly in the variance calculation below, where I will skip explaining them during the application.

To compute the variance of $$G$$, note that since the matrix $$\Lambda := P\otimes (\phi_i\phi_j^T)$$ is non-symmetric because $$i \neq j$$ and the variance formula only applies to the symmetric matrix, we need to symmetrize $$\Lambda$$ first (see this link) then apply the variance formula as follows:
\begin{align*} & \operatorname{Var}(G) = \operatorname{Var}\left(\xi^T\frac{\Lambda + \Lambda^T}{2}\xi\right) \\ =& 2\operatorname{tr}\left(\frac{1}{4}(\Lambda + \Lambda^T)(I_{(N)} \otimes \Sigma) (\Lambda + \Lambda^T)(I_{(N)} \otimes \Sigma)\right)\\ =& \frac{1}{2}\operatorname{tr}\left((\Lambda + \Lambda^T)(I_{(N)} \otimes \Sigma) (\Lambda + \Lambda^T)(I_{(N)} \otimes \Sigma)\right)\\ =& \frac{1}{2}\operatorname{tr}\left(\left(P \otimes (\phi_i\phi_j^T\Sigma) + P \otimes (\phi_j\phi_i^T\Sigma)\right)^2\right) \\ =& \frac{1}{2}\operatorname{tr}\left(\left(\lambda_j P \otimes (\phi_i\phi_j^T) + \lambda_i P \otimes (\phi_j\phi_i^T)\right)^2\right) \\ =& \frac{1}{2}\lambda_i\lambda_j\operatorname{tr}\left(P \otimes (\phi_i\phi_i^T) + P \otimes (\phi_j\phi_j^T) \right) \\ =& \frac{1}{2}\lambda_i\lambda_j\left(\operatorname{tr}\left(P \otimes (\phi_i\phi_i^T)\right) + \operatorname{tr}\left(P \otimes (\phi_j\phi_j^T) \right)\right) \\ =& \frac{1}{2}\lambda_i\lambda_j\left(\operatorname{tr}(P)\operatorname{tr}(\phi_i\phi_i^T) + \operatorname{tr}(P)\operatorname{tr}(\phi_j\phi_j^T)\right) \\ =& \frac{1}{2}\lambda_i\lambda_j\left(\operatorname{tr}(P)\operatorname{tr}(\phi_i^T\phi_i) + \operatorname{tr}(P)\operatorname{tr}(\phi_j^T\phi_j)\right) \\ =& \frac{1}{2}\lambda_i\lambda_j\left((N - 1) \times 1 + (N - 1) \times 1\right) \\ =& \lambda_i\lambda_j(N - 1). \tag{4}\label{4} \end{align*}

In the above derivation, the fifth equality used the condition $$\Sigma\phi_i = \lambda_i\phi_i, \Sigma\phi_j = \lambda_j\phi_j$$. The sixth equality used the fact $$\phi_i^T\phi_j = \phi_j^T\phi_i = 0$$. The tenth equality used the facts $$\operatorname{tr}(P) = N - 1$$ and $$\phi_i^T\phi_i = \phi_j^T\phi_j = 1$$.

Combining $$\eqref{2}, \eqref{3}, \eqref{4}$$ then yields \begin{align*} E[(\phi_i\hat{\Sigma}\phi_j)^2] = \frac{1}{(N - 1)^2}E[G^2] = \frac{1}{(N - 1)^2}(\operatorname{Var}(G) + (E[G])^2) = \frac{\lambda_i\lambda_j}{N - 1}, \end{align*} as desired.