PCA Inference and Asymptotics · Tutorial 510 of 975

Asymptotic Covariance of Principal Directions

Derive how sampling fluctuations in a covariance matrix translate into joint uncertainty in its estimated principal directions.

Advanced 11 min read

What You'll Learn

State conditions for a joint asymptotic normal limit of sample eigenvectors
Derive the first-order perturbation of a signed eigenvector
Express cross-covariances through covariance fluctuations and eigenvalue gaps
Obtain the simpler directional covariance formula for Gaussian observations
Calculate how a non-Gaussian fourth moment changes directional uncertainty
Recognize the sign and orthogonality constraints in joint eigenvector limits

Why Principal Directions Have Sampling Uncertainty

A sample principal direction is an estimate of a population direction, and its precision depends on more than the sample size. Even a small perturbation of a covariance matrix can substantially rotate an eigenvector when its eigenvalue is close to a neighboring eigenvalue. The asymptotic covariance makes this relationship precise: fluctuations in the covariance matrix are projected onto the directions perpendicular to each eigenvector and scaled by eigenvalue gaps.

We consider a fixed dimension $p$ and a population covariance matrix with distinct eigenvalues. The eigenvectors are defined only up to sign, so each sample eigenvector will be signed to point in the same direction as its population counterpart. This convention gives a local, signed limit; it does not remove the underlying sign ambiguity of principal components.

Definition: Let $X_1,\ldots,X_n$ be independent and identically distributed $p$-dimensional random vectors with mean $\mu$, covariance matrix $\Sigma$, and $E\|X_i-\mu\|^4<\infty$. Write the spectral decomposition as $$ \Sigma=\sum_{j=1}^p \lambda_j u_j u_j', \qquad \lambda_1>\lambda_2>\cdots>\lambda_p, $$ where the $u_j$ are orthonormal unit vectors. Let $\hat\lambda_j,\hat u_j$ be the correspondingly ordered eigenvalues and unit eigenvectors of the sample covariance matrix $S_n$, with the sign of $\hat u_j$ chosen so that $\hat u_j'u_j\geq 0$. The principal-direction error is $\hat u_j-u_j$. Its asymptotic covariance is the covariance matrix of the Gaussian limit of $\sqrt n(\hat u_j-u_j)$, when that limit exists.

For specificity, take $S_n=n^{-1}\sum_{i=1}^n(X_i-\bar X)(X_i-\bar X)'$. Using the divisor $n-1$ instead gives the same first-order limit. Let $Y_i=X_i-\mu$, and express its coordinates in the population eigenbasis as $Y_{ia}=u_a'Y_i$. The matrix fluctuation has a multivariate central limit theorem when the fourth moments are finite. Because the sample mean correction is negligible at this scale, estimating $\mu$ does not change the limit. We use this covariance-matrix limit to derive the joint direction limit.

First-Order Change in an Eigenvector

Let $E_n=S_n-\Sigma$. Under the stated conditions, $E_n=O_p(n^{-1/2})$. For a simple eigenvalue, the first-order change in $u_j$ has no component along $u_j$: normalization fixes that component to first order. Its component along each other $u_k$ is the covariance perturbation in that direction divided by the corresponding eigenvalue gap.

Theorem: Suppose $\lambda_j$ is distinct from every other eigenvalue and $\|E_n\|=O_p(n^{-1/2})$. With the sign convention above, $$ \hat u_j-u_j = \sum_{k\ne j}u_k \frac{u_k'E_nu_j}{\lambda_j-\lambda_k} +O_p(\|E_n\|^2). $$

To see the expansion, write $\hat u_j=a_j u_j+\sum_{k\ne j}b_{kj}u_k$, where $a_j=\hat u_j'u_j\geq0$. The sample eigenvalue equation, projected onto $u_k$, gives $$ (\lambda_k-\hat\lambda_j)b_{kj}+u_k'E_n\hat u_j=0, \qquad b_{kj}=\frac{u_k'E_n\hat u_j}{\hat\lambda_j-\lambda_k}. $$ The eigenvalue perturbation bound gives $\hat\lambda_j-\lambda_j=O_p(\|E_n\|)$; the fixed positive gaps then keep every denominator bounded away from zero with probability tending to one. Thus $b_{kj}=O_p(\|E_n\|)$, so $\hat u_j-u_j=O_p(\|E_n\|)$. In the numerator, replacing $\hat u_j$ by $u_j$ therefore introduces an $O_p(\|E_n\|^2)$ error. Replacing $\hat\lambda_j-\lambda_k$ by $\lambda_j-\lambda_k$ also introduces an $O_p(\|E_n\|^2)$ error, since the numerator is $O_p(\|E_n\|)$ and the denominator changes by $O_p(\|E_n\|)$. Finally, unit length gives $a_j=(1-\sum_{k\ne j}b_{kj}^2)^{1/2}=1+O_p(\|E_n\|^2)$. Combining these coordinate expansions proves the result. $\square$

The Joint Asymptotic Covariance

Define a centered Gaussian symmetric matrix $Z$ by its entries in the population eigenbasis. For indices $a,b,c,d$, let $$ \Gamma_{ab,cd} = \operatorname{Cov}(Y_{ia}Y_{ib},Y_{ic}Y_{id}) = E[Y_{ia}Y_{ib}Y_{ic}Y_{id}] -\lambda_a\lambda_c\,\mathbf{1}_{\{a=b\}}\mathbf{1}_{\{c=d\}}, $$ where $\mathbf{1}_{A}$ is 1 when $A$ holds and 0 otherwise. The covariance-matrix central limit theorem says that the entries of $\sqrt n E_n$, expressed in this basis, converge jointly to the entries of $Z$, with $\operatorname{Cov}(Z_{ab},Z_{cd})=\Gamma_{ab,cd}$. Define the Gaussian limit for each direction by $$ V_j=\sum_{k\ne j}u_k\frac{Z_{kj}}{\lambda_j-\lambda_k}. $$

Theorem: Under the iid, finite-fourth-moment and distinct-eigenvalue conditions above, the signed eigenvectors have the joint limit $$ \sqrt n \begin{pmatrix} \hat u_1-u_1\\ \vdots\\ \hat u_p-u_p \end{pmatrix} \xrightarrow{d} \begin{pmatrix} V_1\\ \vdots\\ V_p \end{pmatrix}. $$ In particular, for any $j,r$, $$ \operatorname{Cov}(V_j,V_r) = \sum_{k\ne j}\sum_{l\ne r} \frac{\Gamma_{kj,lr}} {(\lambda_j-\lambda_k)(\lambda_r-\lambda_l)} u_k u_l'. $$

For the limit, apply the joint central limit theorem to the finitely many entries of $n^{-1/2}\sum_i(Y_iY_i'-\Sigma)$. The sample-centering term is $\bar Y\bar Y'$, and $\sqrt n\|\bar Y\bar Y'\|=O_p(n^{-1/2})=o_p(1)$, since $\sqrt n\bar Y=O_p(1)$. The theorem’s first-order expansion has remainder $O_p(\|E_n\|^2)$, so multiplying that remainder by $\sqrt n$ gives $O_p(n^{-1/2})=o_p(1)$. Applying the continuous linear transformation in the expansion to the matrix central limit theorem gives the stated joint limit. Taking covariances of its finite sums gives the displayed covariance formula. $\square$

The gaps appear as denominators in the covariance: a small gap amplifies the effect of the same covariance fluctuation. For Gaussian observations the formula simplifies further. In the eigenbasis, Gaussian coordinates are independent with variances $\lambda_a$, and $$ \Gamma_{ab,cd} = \lambda_a\lambda_b \left(\mathbf{1}_{\{a=c\}}\mathbf{1}_{\{b=d\}} +\mathbf{1}_{\{a=d\}}\mathbf{1}_{\{b=c\}}\right) \quad (a\ne b,\ c\ne d). $$ Consequently, the marginal covariance of direction $j$ is $$ \operatorname{Cov}(V_j) = \sum_{k\ne j} \frac{\lambda_j\lambda_k}{(\lambda_j-\lambda_k)^2}u_ku_k'. $$ For two different directions, the Gaussian cross-covariance is zero unless their coordinates correspond to the same unordered pair of eigenvectors. In that case, orthogonality makes the corresponding signed coordinates perfectly negatively related. In general, the covariance formula depends on fourth moments, not only on $\Sigma$.

Worked Example: A Two-Dimensional Gaussian Covariance

Suppose $X_i\sim N(\mu,\Sigma)$, where $\Sigma=\operatorname{diag}(4,1)$. Then $\lambda_1=4$, $\lambda_2=1$, and the population directions are $u_1=(1,0)'$ and $u_2=(0,1)'$. For the leading direction, the only perpendicular coordinate is $$ u_2'V_1=\frac{Z_{21}}{4-1}. $$ Since the observations are Gaussian, $\operatorname{Var}(Z_{21})=\lambda_2\lambda_1=4$. Therefore $$ \operatorname{Var}(u_2'V_1)=\frac{4}{3^2}=\frac49. $$ Thus the leading sample direction’s perpendicular error has limiting variance $4/9$ after multiplication by $\sqrt n$; its approximate finite-sample variance is $4/(9n)$. The second direction has the coordinate $u_1'V_2=Z_{12}/(1-4)=-Z_{21}/3$, so it has the same marginal variance $4/9$. The minus sign is required by the opposite gap.

Worked Example: All Directional Variances in Three Dimensions

Let the population covariance be Gaussian with eigenvalues $9,4,1$ and eigenvectors $u_1,u_2,u_3$. For $V_1$, the two perpendicular coordinates have variances $$ \operatorname{Var}(u_2'V_1)=\frac{9\cdot4}{(9-4)^2}=\frac{36}{25}, \qquad \operatorname{Var}(u_3'V_1)=\frac{9\cdot1}{(9-1)^2}=\frac9{64}. $$ Hence $\operatorname{Cov}(V_1)=\frac{36}{25}u_2u_2'+\frac9{64}u_3u_3'$. For the second direction, $$ \operatorname{Var}(u_1'V_2)=\frac{4\cdot9}{(4-9)^2}=\frac{36}{25}, \qquad \operatorname{Var}(u_3'V_2)=\frac{4\cdot1}{(4-1)^2}=\frac49. $$ For the third direction, $$ \operatorname{Var}(u_1'V_3)=\frac{1\cdot9}{(1-9)^2}=\frac9{64}, \qquad \operatorname{Var}(u_2'V_3)=\frac{1\cdot4}{(1-4)^2}=\frac49. $$ These are the six marginal tangent-coordinate variances. The joint limit also has cross-covariances: $u_2'V_1=Z_{21}/5$ and $u_1'V_2=Z_{12}/(-5)$, so their covariance is $-36/25$. The other paired coordinates have covariances $-9/64$ and $-4/9$, respectively; coordinates from different unordered pairs have covariance zero. Thus separate marginal error summaries do not describe the whole joint distribution.

Worked Example: A Fourth Moment Changes the Answer

Keep the covariance $\operatorname{diag}(4,1)$, but use a non-Gaussian construction. Let $A$ and $B$ be independent Rademacher variables, each taking $1$ and $-1$ with equal probability. Independently, let $R$ be 0 or $\sqrt2$, each with probability $1/2$, and set $Y_1=2AR$, $Y_2=BR$. Both means are zero. Since $E[R^2]=1$ and $E[AB]=0$, the variances are $E[Y_1^2]=4$, $E[Y_2^2]=1$, and the covariance is zero. Also $E[R^4]=2$, so $$ \operatorname{Var}(Y_1Y_2) =E[Y_1^2Y_2^2] =4E[R^4] =8. $$ The leading direction’s asymptotic perpendicular variance is therefore $8/(4-1)^2=8/9$, rather than the Gaussian value $4/9$. The covariance eigenvalues and their gap are unchanged; the difference comes from the fourth moment of the observations.

Interpretation and a Common Pitfall

The eigenvalue gap is a geometric measure of how strongly the covariance matrix distinguishes one direction from its neighbors. When the gap is large, a given off-diagonal fluctuation $u_k'E_nu_j$ produces a smaller rotation. When the gap is small, the same fluctuation can produce much greater directional uncertainty. Distinct eigenvalues are essential for this individual-eigenvector limit; when eigenvalues coincide, the individual directions are not uniquely determined, and the relevant limit concerns eigenspaces instead.

Pitfall: using the Gaussian formula for every distribution. Knowing $\Sigma$ gives the eigenvalues and gaps, but does not generally determine $\Gamma$. The correction is to use the fourth-moment covariance formula unless Gaussian sampling is justified.
Pitfall: ignoring the sign convention. A raw sample eigenvector can point in either direction, so an ordinary signed difference from $u_j$ may not be small. The local limit uses sign alignment; sign-invariant inference is a separate issue.
Pitfall: treating directions as jointly independent. For Gaussian data, direction errors in different eigenvector planes can be uncorrelated, but coordinates for a shared plane are tied together by symmetry and orthogonality. The joint covariance formula captures this dependence.

Key takeaway: Under distinct population eigenvalues, each first-order eigenvector error is a sum of covariance fluctuations in perpendicular directions, divided by the corresponding eigenvalue gaps. Its joint asymptotic covariance depends on both those gaps and the fourth moments of the observations.

Check Your Understanding

Use the definitions and results above to work through these questions.

Why does the first-order error in $u_j$ have no component along $u_j$?
For Gaussian data, what happens to the asymptotic variance of a direction coordinate as its eigenvalue gap decreases?
In the three-dimensional example, why is the covariance between $u_2'V_1$ and $u_1'V_2$ negative?
Which moment information, in addition to the population covariance, is needed for the general joint covariance formula?
Why must sample eigenvectors be sign-aligned before using the stated signed limit?

Back to Tutorials