Joint Limits for Distinct Sample Eigenvalues

From Covariance Fluctuations to Eigenvalue Fluctuations

A sample covariance matrix fluctuates around its population counterpart, and its eigenvalues fluctuate as a consequence. When the population eigenvalues are distinct, each one has a well-defined first-order response to a small matrix perturbation: its leading change is the perturbation projected onto that eigenvalue’s population eigenvector. Combining this fact with the covariance-matrix central limit theorem gives a joint limit for all the ordered sample eigenvalues.

The word joint matters. Even when the sample covariance entries have a multivariate normal limit, the eigenvalue limits need not be independent: their dependence is determined by the fourth moments of the observations. Under multivariate normal sampling, however, the limiting eigenvalue fluctuations are independent.

Definition: Let \(X_1,\ldots,X_n\) be independent and identically distributed vectors in \(\mathbb{R}^p\), with mean \(\mu\) and positive-definite covariance matrix \(\Sigma\). Write the ordered population eigenvalues as \(\lambda_1>\lambda_2>\cdots>\lambda_p>0\), with corresponding orthonormal eigenvectors \(u_1,\ldots,u_p\). Define \[ \bar X=\frac{1}{n}\sum_{i=1}^n X_i,\qquad S_n=\frac{1}{n}\sum_{i=1}^n(X_i-\bar X)(X_i-\bar X)'. \] Let \(\hat\lambda_1\geq\hat\lambda_2\geq\cdots\geq\hat\lambda_p\) be the ordered eigenvalues of \(S_n\). The distinct-eigenvalue asymptotic problem is to find the limiting distribution of the vector \(\sqrt n(\hat\lambda-\lambda)\), where \(\hat\lambda=(\hat\lambda_1,\ldots,\hat\lambda_p)'\) and \(\lambda=(\lambda_1,\ldots,\lambda_p)'\).

Throughout, assume \(E\|X_i-\mu\|^4<\infty\). The eigenvalues are strictly separated, and \(p\) is fixed as \(n\) increases. Using the divisor \(n-1\) instead of \(n\) for the sample covariance gives the same first-order limit: the two matrices differ by a factor \(n/(n-1)=1+O(n^{-1})\), which changes their eigenvalues by \(O_p(n^{-1})\), negligible after multiplication by \(\sqrt n\).

Why the First-Order Change Is a Projection

Consider a symmetric matrix \(A\) with distinct eigenvalues \(\lambda_1,\ldots,\lambda_p\), and let \(u_j\) be a unit eigenvector for \(\lambda_j\). If a small symmetric matrix \(E\) is added, the first-order change in the \(j\)th eigenvalue is \(u_j'Eu_j\). This perturbation formula is the key step: it extracts the diagonal fluctuation in the population eigenbasis, while off-diagonal fluctuations affect eigenvalues only at higher order.

Theorem: Let \(A\) be a real symmetric matrix with distinct ordered eigenvalues and let \(\delta_j=\min_{k\ne j}|\lambda_j-\lambda_k|>0\). For a symmetric perturbation \(E\) with sufficiently small operator norm, \[ \lambda_j(A+E)=\lambda_j(A)+u_j'Eu_j+O\!\left(\frac{\|E\|^2}{\delta_j}\right). \] Here the constant in the remainder may be chosen uniformly for \(E\) in a sufficiently small neighborhood of zero.

To justify the remainder, write \(\lambda=\lambda_j(A)\), and choose a unit eigenvector \(v\) of \(A+E\) for the eigenvalue \(\theta\) that remains near \(\lambda\). Decompose \(v=a u_j+w\), with \(u_j'w=0\), choosing its sign so \(a\geq0\). Weyl’s inequality gives \(|\theta-\lambda|\leq\|E\|\). For \(\|E\|<\delta_j/2\), projecting the eigenvector equation onto the orthogonal complement of \(u_j\) gives \[ \|(A-\theta I)w\|\geq \frac{\delta_j}{2}\|w\|, \qquad \|(A-\theta I)w\|\leq\|E v\|\leq\|E\|. \] Thus \(\|w\|\leq 2\|E\|/\delta_j\). For sufficiently small \(\|E\|\), this also ensures \(a=(1-\|w\|^2)^{1/2}\) is bounded away from zero. Projecting the eigenvector equation onto \(u_j\) yields \[ \theta=\lambda+u_j'Eu_j+\frac{u_j'Ew}{a}. \] The last term has absolute value at most \(\|E\|\|w\|/a=O(\|E\|^2/\delta_j)\), proving the expansion. \(\square\)

The Joint Limit Theorem

Let \(Z_i=U'(X_i-\mu)\), where \(U=(u_1,\ldots,u_p)\). In these coordinates, the population covariance is diagonal, and \(E[Z_{ij}^2]=\lambda_j\). Define the \(p\times p\) matrix \(\Omega\) by \[ \Omega_{jk}=\operatorname{Cov}(Z_{1j}^2,Z_{1k}^2) =E[Z_{1j}^2Z_{1k}^2]-\lambda_j\lambda_k,\qquad 1\leq j,k\leq p. \]

Theorem: Under the assumptions above, \[ \sqrt n \begin{pmatrix} \hat\lambda_1-\lambda_1\\ \vdots\\ \hat\lambda_p-\lambda_p \end{pmatrix} \xrightarrow{d}N_p(0,\Omega). \] The limiting multivariate normal distribution may be singular. Its covariance matrix depends on the joint fourth moments of the coordinates in the population eigenbasis.

Here is the full derivation. Put \(Y_i=X_i-\mu\). The centered sample covariance satisfies \[ S_n=\frac1n\sum_{i=1}^nY_iY_i'-\bar Y\bar Y', \qquad \bar Y=\frac1n\sum_{i=1}^nY_i. \] Since \(E\|\bar Y\|^2=O(n^{-1})\), we have \(\bar Y=O_p(n^{-1/2})\), and hence \(\|\bar Y\bar Y'\|=O_p(n^{-1})\). Therefore \[ \sqrt n\left\|S_n-\frac1n\sum_{i=1}^nY_iY_i'\right\|=o_p(1). \] The multivariate central limit theorem applied to the vector of diagonal entries in the \(U\)-basis gives \[ \sqrt n\big(u_j'(S_n-\Sigma)u_j\big)_{j=1}^p \xrightarrow{d}N_p(0,\Omega). \] The finite fourth-moment assumption guarantees finite variances for these entries.

Also, \(S_n-\Sigma=O_p(n^{-1/2})\) in operator norm because \(p\) is fixed and the covariance entries satisfy the multivariate central limit theorem. Apply the perturbation expansion to \(A=\Sigma\) and \(E=S_n-\Sigma\). For each \(j\), \[ \hat\lambda_j-\lambda_j =u_j'(S_n-\Sigma)u_j+O_p(\|S_n-\Sigma\|^2) =u_j'(S_n-\Sigma)u_j+O_p(n^{-1}). \] Multiplying the remainder by \(\sqrt n\) makes it \(o_p(1)\). The joint limit above therefore transfers to the ordered eigenvalues. Since the population eigenvalues are distinct, their gaps are positive; with probability tending to one, the sample eigenvalues remain in the corresponding order. This proves the theorem. \(\square\)

Worked Examples

Worked Example: The Gaussian Case

Suppose \(X_i\sim N_p(\mu,\Sigma)\), with distinct eigenvalues \(\lambda_1,\ldots,\lambda_p\). In the population eigenbasis, the coordinates \(Z_{ij}\) are independent normal variables with \(Z_{ij}\sim N(0,\lambda_j)\). For \(j=k\), a centered normal variable has fourth moment \(3\lambda_j^2\), so \[ \Omega_{jj}=E[Z_{1j}^4]-\lambda_j^2=3\lambda_j^2-\lambda_j^2=2\lambda_j^2. \] For \(j\ne k\), independence gives \(E[Z_{1j}^2Z_{1k}^2]=\lambda_j\lambda_k\), hence \(\Omega_{jk}=0\). Thus the joint limit is normal with diagonal covariance matrix, so its coordinates are independent.

For example, if the population eigenvalues are \(9,4,1\), then \[ \sqrt n \begin{pmatrix} \hat\lambda_1-9\\ \hat\lambda_2-4\\ \hat\lambda_3-1 \end{pmatrix} \xrightarrow{d} N_3\!\left(0, \begin{pmatrix} 162&0&0\\ 0&32&0\\ 0&0&2 \end{pmatrix}\right). \] At \(n=100\), the corresponding first-order standard errors are \(\sqrt{162/100}\approx1.273\), \(\sqrt{32/100}\approx0.566\), and \(\sqrt{2/100}\approx0.141\). These are asymptotic standard errors, not exact finite-sample standard deviations.

Worked Example: Fourth Moments Can Couple Eigenvalue Limits

Let \(R^2\) equal \(1/2\) or \(3/2\), each with probability \(1/2\), and let \(\varepsilon_1,\varepsilon_2\) be independent random signs, also independent of \(R\). Set \(Z_1=2R\varepsilon_1\) and \(Z_2=R\varepsilon_2\). The signs give mean zero and zero covariance between the coordinates. Since \(E[R^2]=1\), the covariance matrix is \(\operatorname{diag}(4,1)\), with distinct eigenvalues.

Also, \(E[R^4]=\{(1/2)^2+(3/2)^2\}/2=5/4\). The entries of the limiting covariance are therefore \[ \Omega_{11}=16E[R^4]-4^2=20-16=4,\qquad \Omega_{22}=E[R^4]-1^2=\frac14, \] and \[ \Omega_{12}=4E[R^4]-(4)(1)=5-4=1. \] Consequently, the joint limit has covariance matrix \[ \Omega=\begin{pmatrix}4&1\\1&1/4\end{pmatrix}. \] This matrix is singular: its determinant is \(4(1/4)-1=0\). Indeed, \(Z_1^2=4Z_2^2\) for every observation, so the first-order diagonal fluctuations are tied together. This example shows why zero population covariance between coordinates does not by itself imply independent eigenvalue limits.

Worked Example: First-Order Prediction Versus an Exact Perturbation

Take \(\Sigma=\operatorname{diag}(9,4)\) and the symmetric perturbation \[ E=\begin{pmatrix}0.1&0.3\\0.3&0\end{pmatrix}. \] The first-order formula predicts changes \(u_1'Eu_1=0.1\) and \(u_2'Eu_2=0\), because the population eigenvectors are the coordinate axes. The perturbed matrix is \[ \Sigma+E=\begin{pmatrix}9.1&0.3\\0.3&4\end{pmatrix}. \] Its trace is \(13.1\), and its determinant is \(9.1(4)-0.3^2=36.31\). Hence its eigenvalues are \[ \frac{13.1\pm\sqrt{13.1^2-4(36.31)}}{2} =\frac{13.1\pm\sqrt{26.37}}{2} \approx 9.1176,\ 3.9824. \] The exact changes are approximately \(0.1176\) and \(-0.0176\), not \(0.1\) and \(0\). The off-diagonal perturbation changes the eigenvalues at second order; it does not enter their first-order terms.

Interpretation and Common Pitfall

The theorem converts matrix-level sampling uncertainty into uncertainty for the spectrum. Its covariance formula also indicates what information is needed: the joint fourth moments of observations projected onto the population principal axes. In the Gaussian case these moments reduce to the familiar diagonal covariance \(2\lambda_j^2\). Outside that case, the off-diagonal entries of \(\Omega\) need not vanish.

Pitfall: treating all eigenvalue limits as independent. A diagonal population covariance does not imply independent squared coordinates, and uncorrelated coordinates need not have uncorrelated squares. The correction is to calculate \(\Omega_{jk}=E[Z_j^2Z_k^2]-\lambda_j\lambda_k\); independence follows in the Gaussian case, not merely from zero covariance.
Pitfall: using the first-order formula as an exact finite-sample identity. The projection \(u_j'Eu_j\) is the linear term. Off-diagonal perturbations generally contribute at second order, as the numerical example shows.
Pitfall: applying the formula at a repeated eigenvalue. The positive-gap condition is essential to this eigenvalue-by-eigenvalue expansion. When eigenvalues coincide, individual axes are not uniquely identified and the analysis must instead account for the whole cluster or eigenspace.

Thus, for distinct population eigenvalues, the joint eigenvalue limit is obtained by taking the diagonal entries of the covariance-matrix fluctuation in the population eigenbasis. The distinct gaps make this first-order description valid, while the fourth moments determine how the resulting eigenvalue fluctuations move together.

Check Your Understanding

Use the assumptions and results in this tutorial to answer the following questions.

Why does estimating the mean and centering the observations not change the first-order eigenvalue limit?
For two distinct population eigenvalues, write the limiting covariance between their eigenvalue fluctuations in terms of eigenbasis coordinates.
Under multivariate normal sampling, what is the asymptotic variance of the \(j\)th sample eigenvalue after multiplication by \(\sqrt n\)?
Why can two eigenvalue limits be dependent even when the corresponding population coordinates have zero covariance?
Which condition in the perturbation argument fails when two population eigenvalues are equal?