Confidence Regions for Leading Eigenspaces

Why infer an eigenspace rather than its individual directions?

In principal component analysis, the leading \(k\)-dimensional eigenspace of a population covariance matrix \(\Sigma\) describes a \(k\)-dimensional subspace of directions. If the eigenvalues inside that group are equal, there is no unique choice of individual eigenvectors within the subspace. Even when they are distinct, eigenvector signs are arbitrary. The orthogonal projector onto the subspace avoids both issues: it is unchanged by sign changes or by a change of orthonormal basis within the subspace.

A confidence region for this projector must account for sampling variability and for the eigengap separating the leading group from the rest. The central point is to keep the coordinates used for the asymptotic limit fixed. Coordinates expressed in a sample eigenbasis can behave differently if that basis rotates randomly.

Write \(U=(U_1,U_2)\) for a fixed orthogonal eigenbasis of \(\Sigma\), with \(U_2\) spanning the complementary eigenspace. The basis is fixed at the population covariance for theoretical coordinates. Given a sample covariance \(S\), let \(\hat P\) be its projector onto the leading \(k\) eigenvectors. The matrix \(U_2^T\hat P U_1\) is a local coordinate for the estimated subspace relative to the population basis; it is not itself invariant to a change of basis in either group.

First-order projector coordinates and covariance

Assume \(p\) and \(k\) are fixed as the sample size \(n\) grows, the boundary eigengap is positive, and \(\sqrt n\,\operatorname{vec}(S-\Sigma)\xrightarrow{d}G\), where \(G\) is a mean-zero Gaussian vector with covariance matrix \(\Omega\). Here \(\operatorname{vec}\) stacks matrix columns. For example, this central limit theorem holds for the usual sample covariance under independent, identically distributed observations with finite fourth moments, with the appropriate covariance \(\Omega\). No distinctness of eigenvalues within either group is required.

Proof. Work first in the fixed eigenbasis \(U\), so that \(\Sigma\) is diagonal. Write \(\hat P\) in leading and complementary blocks as \(\begin{pmatrix}A&B^T\\B&D\end{pmatrix}\). The identity \(\hat P^2=\hat P\) implies \(BA+DB=B\). Since \(\hat P\) converges to \(P\), \(A\) converges to the identity and \(D\) to zero. The same identity, or the first-order perturbation equation for the invariant subspace, shows that the cross block is first order in \(S-\Sigma\); products of two perturbations are negligible after multiplication by \(\sqrt n\). The cross-block equation at first order is \[ B_{ai}=\frac{(S-\Sigma)_{k+a,i}}{\lambda_i-\lambda_{k+a}}+o_p(n^{-1/2}). \] Each denominator is nonzero because \(\lambda_i\geq\lambda_k>\lambda_{k+1}\geq\lambda_{k+a}\). This equation also follows by differentiating the eigenvector equation and projecting onto the complementary eigenspace. Applying the assumed multivariate central limit theorem to the finitely many entries gives the stated Gaussian limit and covariance. The formula is valid for repeated eigenvalues within either group: the basis is fixed, and the cross-group denominators remain separated from zero. \(\square\)

The result is a fixed-coordinate delta-method statement. It does not automatically imply a Gaussian limit for \( \hat U_2^T\hat P\hat U_1\), where \(\hat U_1,\hat U_2\) are estimated eigenvector bases. If eigenvalues are repeated within a group, those bases may rotate according to the same random perturbation that determines the cross-block entries. In general, a data-dependent rotation of a Gaussian limit need not remain Gaussian.

Wald regions for projectors

A coordinate ellipsoid based on \(V\) is natural for theory, but the population basis \(U\) is unknown. A way to state the region without choosing an estimated basis for its coordinates is to work in the ambient space of \(p\times p\) matrices. Let \(D_\Sigma\) denote the derivative of the map from a covariance matrix to its leading projector. For a symmetric perturbation \(E\), its action is \[ D_\Sigma[E]=\sum_{i=1}^k\sum_{j=k+1}^p \frac{u_j^T E u_i}{\lambda_i-\lambda_j} \left(u_j u_i^T+u_i u_j^T\right). \] This is the first-order change in the projector. Although the displayed expression uses an eigenbasis, the resulting linear map is unchanged by a different orthonormal basis within an exactly repeated eigenspace.

Let \(\mathcal W\) be the covariance matrix of \(\operatorname{vec}(D_\Sigma[G])\). It is generally singular as a \(p^2\times p^2\) matrix: first-order projector changes lie in a \(d=k(p-k)\)-dimensional tangent space. Suppose their covariance is positive definite on that tangent space, so \(\mathcal W\) has rank \(d\). Let \(\hat\Omega\) consistently estimate \(\Omega\), and construct \(\hat{\mathcal W}\) by substituting \(S\) and \(\hat\Omega\) into the derivative covariance formula. The derivative is continuous as long as the boundary eigengap stays positive; repeated eigenvalues within a group do not invalidate this ambient construction.

Justification. The projector derivative and the delta method give \(\sqrt n\,\operatorname{vec}(\hat P-P)\xrightarrow{d}N(0,\mathcal W)\). Consistency of \(S\) and \(\hat\Omega\), together with the nonzero eigengap, gives \(\hat{\mathcal W}\xrightarrow{p}\mathcal W\). On the \(d\)-dimensional range of \(\mathcal W\), the quadratic form using its Moore–Penrose inverse therefore converges in distribution to \(\chi^2_d\): after diagonalizing the covariance on that range, it is the sum of squares of \(d\) independent standard normal coordinates. Applying this convergence to \(Q=P\), whose residual is \(\hat P-P\), proves the coverage statement because the chi-square distribution has a continuous distribution function at its \(1-\alpha\) quantile. \(\square\)

The ambient formulation uses the full covariance of the projector perturbation and does not require reporting a random orientation for the estimated eigenvectors. In practice, \(\hat\Omega\) must match the sampling scheme and covariance estimator. For non-Gaussian observations, a covariance estimate based only on a Gaussian Wishart formula may be incorrect; the relevant fourth-moment variability must be represented.

Common mistake and practical interpretation

A common mistake is to compute a covariance matrix in a population-aligned basis, rotate the coordinates into the estimated eigenbasis, and then declare the rotated coordinates jointly Gaussian by the ordinary delta method. That reasoning silently treats the rotation as fixed. With repeated within-group eigenvalues, the estimated orientation can itself have a non-degenerate random limit and can depend on the perturbation being rotated. A data-dependent transformation does not in general preserve a Gaussian limit.

The correction is to state the delta-method limit in a fixed population basis, where the derivative and covariance are well defined, and to use an ambient projector covariance for an implementable, basis-invariant region. Estimated-basis coordinates can be used with a Gaussian Wald calibration when additional conditions ensure a deterministic limiting basis—for example, distinct eigenvalues within each group together with consistent sign conventions—or when a separate argument establishes the required limit. A boundary eigengap is essential in all cases: without it, the leading \(k\)-dimensional projector itself may not have the regular first-order behavior used here.

An asymptotic confidence region describes repeated-sampling uncertainty under the specified model and regularity conditions. It is not a guarantee that every direction inside a plotted region is equally plausible, and its accuracy can be poor in finite samples when the estimated boundary gap is small. Its value is that it represents uncertainty about the subspace as a whole, rather than assigning unwarranted certainty to a particular sign or within-subspace orientation.