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.

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

1 Answer 1


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.