Worked Problems in PCA Asymptotic Inference

From a Covariance Matrix to an Inference Calculation

PCA produces estimated eigenvalues and directions, but an estimate alone does not describe its sampling uncertainty. In Gaussian models, asymptotic approximations turn a sample covariance matrix into standard errors and intervals. The calculations depend on both the estimated eigenvalues and, for directions, the gaps between them.

Let $X_1,\ldots,X_n$ be independent observations from $N_p(\mu,\Sigma)$, where $p$ is fixed. Write the sample covariance matrix as $S$, using the usual divisor $n-1$, and order its eigenvalues as $\hat\lambda_1\geq\cdots\geq\hat\lambda_p$. The corresponding population eigenvalues are $\lambda_1\geq\cdots\geq\lambda_p$. The results below are large-sample approximations; they are not exact finite-sample confidence statements.

Working assumptions: The observations are independent and identically distributed Gaussian vectors, the dimension $p$ is fixed as $n$ grows, and the population covariance matrix has the distinct eigenvalues needed for the particular eigenvalue or direction being studied. For inference on the leading direction, in particular, assume $\lambda_1>\lambda_2$. Write $v_j$ for a unit population eigenvector corresponding to $\lambda_j$.

Two reusable approximations organize the examples. For a simple eigenvalue, the asymptotic standard error is estimated by $\hat\lambda_j\sqrt{2/n}$. For the leading direction, uncertainty along $v_j$, for $j>1$, is governed by the product $\lambda_1\lambda_j$ divided by the squared eigenvalue gap. These approximations follow from first-order perturbation of the sample covariance matrix.

Two First-Order Results

Theorem: Under the working assumptions, for a simple population eigenvalue $\lambda_j$, $$ \sqrt n(\hat\lambda_j-\lambda_j)\xrightarrow{d}N(0,2\lambda_j^2). $$ If $\lambda_1>\lambda_2$, choose the sign of the estimated leading unit eigenvector $\hat v_1$ so that $\hat v_1^Tv_1\geq 0$. For each $j>1$, $$ \sqrt n\,v_j^T\hat v_1\xrightarrow{d} N\left(0,\frac{\lambda_1\lambda_j}{(\lambda_1-\lambda_j)^2}\right). $$

Proof. Work in the orthonormal population eigenbasis, where $\Sigma$ is diagonal with entries $\lambda_1,\ldots,\lambda_p$. Put $E=S-\Sigma$. For fixed $p$, the multivariate central limit theorem applies to the entries of $E$. In this basis, the Gaussian fourth-moment identity gives $$ \operatorname{Cov}(\sqrt n E_{ab},\sqrt n E_{cd}) =\Sigma_{ac}\Sigma_{bd}+\Sigma_{ad}\Sigma_{bc}. $$ Using $n-1$ rather than $n$ in $S$ does not change this limit. In particular, $\sqrt n E_{jj}$ converges to a mean-zero normal variable with variance $2\lambda_j^2$, and, for $j\ne 1$, $\sqrt n E_{j1}$ converges to a mean-zero normal variable with variance $\lambda_j\lambda_1$.

For a simple eigenvalue, first-order eigenvalue perturbation gives $\hat\lambda_j-\lambda_j=v_j^TEv_j+O_p(\lVert E\rVert^2)$. For the leading eigenvector, first-order perturbation gives $$ v_j^T\hat v_1 =\frac{v_j^TEv_1}{\lambda_1-\lambda_j} +O_p(\lVert E\rVert^2),\qquad j>1. $$ These expansions apply with probability tending to one because $\lVert E\rVert=O_p(n^{-1/2})$ and the relevant eigenvalue gaps are positive. Thus $\sqrt n\,\lVert E\rVert^2=o_p(1)$, so the remainder vanishes after multiplication by $\sqrt n$. Applying the stated central limit limits to the first-order terms proves both claims. $\square$

The results suggest plug-in standard errors: replace population eigenvalues by sample eigenvalues and use $n$ in the denominator. For a direction, the eigenvalue gap is squared in the denominator; a small gap therefore means greater uncertainty. The sign convention makes a local comparison between directions possible, but does not turn an eigenvector into an ordinary signed parameter.

Worked Example: Standard Error for the Leading Eigenvalue

Suppose $n=100$ observations give the sample covariance matrix $$ S=\begin{pmatrix}4&1\\1&2\end{pmatrix}. $$ Calculate a plug-in standard error and a 95% asymptotic confidence interval for the leading population eigenvalue.

The eigenvalues solve $\det(S-\lambda I)=0$. The trace is $6$, the determinant is $7$, and hence $$ \hat\lambda_{1,2} =\frac{6\pm\sqrt{6^2-4(7)}}{2} =3\pm\sqrt2. $$ Thus $\hat\lambda_1=3+\sqrt2\approx4.4142$.

The plug-in standard error is $$ \widehat{\operatorname{se}}(\hat\lambda_1) =\hat\lambda_1\sqrt{\frac{2}{n}} =4.4142\sqrt{\frac{2}{100}} \approx0.6243. $$ Using the standard normal 95% critical value $1.96$, the interval is $$ 4.4142\pm1.96(0.6243) =4.4142\pm1.2236 \approx(3.1906,\ 5.6378). $$ This is a large-sample interval for the population leading eigenvalue under the Gaussian model. Its precision depends on the approximation being adequate at $n=100$; the calculation does not make the interval exact.

Uncertainty in a Principal Direction

A principal direction is a unit vector, so its uncertainty is naturally described locally in directions perpendicular to it. After choosing the sign of $\hat v_1$ to align with $v_1$, the coefficient $v_j^T\hat v_1$ measures its small deviation toward the population direction $v_j$. The plug-in standard error for this coefficient is $$ \widehat{\operatorname{se}}_j =\frac{\sqrt{\hat\lambda_1\hat\lambda_j}} {\sqrt n\,(\hat\lambda_1-\hat\lambda_j)},\qquad j>1. $$ This approximation is most useful when the leading eigenvalue is separated from the rest. It should not be interpreted as a global, sign-sensitive confidence interval for every coordinate of the eigenvector.

Worked Example: Component-Wise Uncertainty for a Leading Direction

Suppose $n=200$ and the sample covariance matrix is $$ S=\begin{pmatrix}9&0&0\\0&4&0\\0&0&1\end{pmatrix}. $$ The leading sample direction is $\hat v_1=(1,0,0)^T$, and the sample eigenvalues are $9,4,1$. Calculate the uncertainty in the leading direction along the second and third population eigendirections, using plug-in values.

For the second direction, the estimated standard error of the tangent coefficient is $$ \widehat{\operatorname{se}}_2 =\frac{\sqrt{9\cdot4}}{\sqrt{200}(9-4)} =\frac{6}{5\sqrt{200}} \approx0.0849. $$ For the third direction, it is $$ \widehat{\operatorname{se}}_3 =\frac{\sqrt{9\cdot1}}{\sqrt{200}(9-1)} =\frac{3}{8\sqrt{200}} \approx0.0265. $$ A component-wise, local 95% normal approximation gives ranges of about $\pm1.96(0.0849)=\pm0.1664$ and $\pm1.96(0.0265)=\pm0.0520$, respectively. These are approximations for the projections onto the two tangent directions, not a joint confidence region. The second direction is less precisely determined: its gap from the leading eigenvalue is $5$, compared with a gap of $8$ for the third direction.

Explained-Variance Shares

The leading explained-variance share is the leading eigenvalue divided by the total variance. Define the population share $\rho_1=\lambda_1/T$, where $T=\sum_{j=1}^p\lambda_j=\operatorname{tr}(\Sigma)$, and estimate it by $\hat\rho_1=\hat\lambda_1/\sum_j\hat\lambda_j$. Since the sum of the sample eigenvalues is $\operatorname{tr}(S)$, this estimate is also $\hat\lambda_1/\operatorname{tr}(S)$.

Result: Under the Gaussian working assumptions, with a simple leading eigenvalue, $$ \sqrt n(\hat\rho_1-\rho_1)\xrightarrow{d}N(0,V_\rho), \qquad V_\rho=\frac{2\lambda_1^2}{T^4} \left[(T-\lambda_1)^2+\sum_{j=2}^p\lambda_j^2\right]. $$ A plug-in standard error is obtained by replacing each population eigenvalue and $T$ by its sample estimate.

Derivation. In the population eigenbasis, the first-order change in the share is $$ \hat\rho_1-\rho_1 =\frac{E_{11}}{T}-\frac{\lambda_1\operatorname{tr}(E)}{T^2} +o_p(n^{-1/2}) =\frac{(T-\lambda_1)E_{11}-\lambda_1\sum_{j=2}^pE_{jj}}{T^2} +o_p(n^{-1/2}). $$ For Gaussian observations, the diagonal terms have asymptotic variances $2\lambda_j^2/n$, and distinct diagonal terms have asymptotic covariance zero. Therefore the variance of the displayed first-order term, multiplied by $n$, is $$ \frac{2(T-\lambda_1)^2\lambda_1^2+ 2\lambda_1^2\sum_{j=2}^p\lambda_j^2}{T^4}=V_\rho. $$ The central limit theorem and the delta method give the asserted normal limit. A Wald interval uses $\hat\rho_1\pm1.96\widehat{\operatorname{se}}(\hat\rho_1)$; it is an asymptotic interval, and it may extend outside the possible range $[0,1]$ in small samples. $\square$

Worked Example: A Confidence Interval for the Leading Variance Share

Suppose $n=150$ and $$ S=\begin{pmatrix}6&0&0\\0&3&0\\0&0&1\end{pmatrix}. $$ The estimated leading share is $\hat\rho_1=6/(6+3+1)=0.6$. Calculate its plug-in standard error and a 95% asymptotic interval.

Here $\hat T=10$, $T-\lambda_1$ is estimated by $4$, and $\sum_{j=2}^p\hat\lambda_j^2=3^2+1^2=10$. The estimated asymptotic variance is $$ \widehat V_\rho =\frac{2(6)^2}{(10)^4}\left[4^2+10\right] =\frac{72\cdot26}{10000} =0.1872. $$ Because $V_\rho$ is the variance in the limit for $\sqrt n(\hat\rho_1-\rho_1)$, the estimated variance of $\hat\rho_1$ is $0.1872/150=0.001248$. Thus $$ \widehat{\operatorname{se}}(\hat\rho_1) =\sqrt{0.001248}\approx0.0353. $$ The 95% Wald interval is $$ 0.6\pm1.96(0.0353)\approx(0.5308,\ 0.6692). $$ The interval quantifies uncertainty in the fraction of total population variance represented by the leading component, not uncertainty in the eigenvalue alone.

Common Mistakes and Practical Interpretation

A tempting shortcut is to use the eigenvalue standard error to describe uncertainty in the corresponding eigenvector. These are different inferential targets. Eigenvalue uncertainty depends on the eigenvalue’s size, while direction uncertainty depends strongly on the gap to neighbouring eigenvalues. A large leading eigenvalue does not by itself guarantee a stable direction if the second eigenvalue is nearly as large.

Using the wrong denominator. The standard error for an eigenvalue is $\hat\lambda_j\sqrt{2/n}$, not $\hat\lambda_j\sqrt{2}$. The factor $n$ is the source of the sampling precision.
Ignoring the eigenvalue gap. For a direction, omitting $\hat\lambda_1-\hat\lambda_j$ can severely understate uncertainty. If the leading gap is very small, the local normal approximation may be unreliable.
Treating a vector’s sign as meaningful. The vectors $v$ and $-v$ represent the same principal direction. Align signs before comparing estimated and population directions; a sign flip alone is not sampling error.
Calling a Wald interval exact. The intervals here rely on large-sample normal approximations under a Gaussian model. Their stated coverage is approximate, not guaranteed for a particular finite sample.

The practical lesson is to identify the inferential target before calculating a standard error. For eigenvalues, use the eigenvalue approximation; for directions, inspect the relevant gaps; for explained-variance shares, account for uncertainty in both the leading eigenvalue and the total variance.

Check Your Understanding

Use the assumptions and calculations in this tutorial to answer the following.

For a simple eigenvalue estimate of $5$ from $n=100$ Gaussian observations, what plug-in asymptotic standard error follows from the eigenvalue result?
Why can two estimated leading directions have different uncertainty even if their leading eigenvalues are equal?
For $S=\operatorname{diag}(8,6,2)$, which gap governs the leading direction’s uncertainty along the second sample eigendirection?
In the explained-share calculation, why does the variance include uncertainty from eigenvalues other than the leading one?
What does sign alignment resolve when comparing eigenvectors, and what uncertainty does it not remove?