Tutorials › Multivariate Statistics (Graduate) › Bootstrap Distributions of Eigenvalues

Resampling and Numerical Methods · Tutorial 961 of 975

Bootstrap Distributions of Eigenvalues

Derive when bootstrap distributions for covariance eigenvalues are reliable, and see how repeated roots and shrinking eigengaps change their behavior.

Advanced 11 min read

What You'll Learn

  • Define the empirical bootstrap distribution of sample covariance eigenvalues
  • Derive the first-order perturbation of an isolated eigenvalue
  • Establish bootstrap consistency for a simple population eigenvalue
  • Describe the nonregular limit at a repeated eigenvalue
  • Use eigengap size to diagnose near-repeated-root behavior

Why Eigenvalue Uncertainty Needs Special Care

Eigenvalues of a sample covariance matrix summarize variation along data-adaptive directions. They are used, for example, to describe principal component variances. Resampling observations can quantify their sampling uncertainty, but the validity of the usual bootstrap depends on the population spectrum. An eigenvalue separated from the others changes smoothly under small perturbations. At a repeated eigenvalue, the associated directions can rotate within a whole subspace, and the ordered eigenvalues need not have a smooth, linear first-order approximation.

The key diagnostic is the eigengap: the distance from the eigenvalue of interest to the other population eigenvalues. A fixed positive gap gives a regular setting. A zero gap gives a nonregular one, and a gap shrinking toward zero can put the problem between those regimes.

Bootstrap Setup and Eigenvalue Perturbation

Definition: Let \(X_1,\ldots,X_n\) be independent and identically distributed \(p\)-dimensional random vectors with mean \(\mu\) and covariance matrix \(\Sigma\). Write their sample covariance matrix as $$ S_n=\frac{1}{n-1}\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^\top, \qquad \bar X=\frac{1}{n}\sum_{i=1}^n X_i. $$ Order eigenvalues from largest to smallest, with multiplicity, and write \(\hat\lambda_1\ge\cdots\ge\hat\lambda_p\) for those of \(S_n\), and \(\lambda_1\ge\cdots\ge\lambda_p\) for those of \(\Sigma\). In the ordinary empirical bootstrap, draw \(X_1^*,\ldots,X_n^*\) independently from the empirical distribution assigning mass \(1/n\) to each observed vector. Let \(S_n^*\) and \(\hat\lambda_j^*\) be the resulting sample covariance and ordered eigenvalues. The bootstrap distribution for root \(j\) is the conditional distribution, given the data, of \(\sqrt n(\hat\lambda_j^*-\hat\lambda_j)\).

A basic matrix perturbation fact explains the distinction between isolated and repeated roots. For a real symmetric matrix \(A\), let \(E\) be another real symmetric matrix. If \(\lambda\) is a simple eigenvalue of \(A\), with unit eigenvector \(v\) and gap \(\gamma=\min_{\eta\ne\lambda}|\lambda-\eta|>0\), then, for sufficiently small \(\|E\|\), $$ \lambda(A+E)=\lambda+v^\top E v+O\left(\frac{\|E\|^2}{\gamma}\right). $$ Here \(\|\cdot\|\) is the operator norm. The linear term depends on the single eigenvector \(v\); the gap controls the remainder.

Result and proof: If \(\|E\|\le\gamma/4\), the eigenvalue of \(A+E\) corresponding to \(\lambda\) satisfies $$ \left|\lambda(A+E)-\lambda-v^\top E v\right| \le \frac{2\|E\|^2}{\gamma}. $$ To see this, choose coordinates with \(v\) first, so \(A\) is block diagonal with first entry \(\lambda\), and write the perturbed eigenvector as \((\alpha,z)^\top\). Weyl's eigenvalue bound gives \(|\delta|\le\|E\|\), where the perturbed eigenvalue is \(\lambda+\delta\). The lower block of the eigenvector equation is $$ (A_{-}-\lambda I-\delta I+E_{--})z=-E_{-1}\alpha. $$ The matrix on the left has smallest singular value at least \(\gamma-|\delta|-\|E_{--}\|\ge\gamma/2\). Thus \(\|z/\alpha\|\le 2\|E\|/\gamma\). The first block equation gives $$ \delta=E_{11}+E_{1-}z/\alpha. $$ Since \(E_{11}=v^\top E v\), the difference between \(\delta\) and \(v^\top E v\) is at most \(2\|E\|^2/\gamma\), as claimed. \(\square\)

Isolated Eigenvalues: A Regular Bootstrap Limit

Suppose \(\lambda_j\) is simple and separated from every other eigenvalue of \(\Sigma\) by a fixed positive gap. Let \(v_j\) be its unit eigenvector. Assume \(E\|X_i\|^4<\infty\). The covariance central limit theorem, together with the perturbation expansion, gives $$ \sqrt n(\hat\lambda_j-\lambda_j) = v_j^\top \sqrt n(S_n-\Sigma)v_j+o_p(1) \xrightarrow{d} N\left(0,\operatorname{Var}\left([v_j^\top(X_i-\mu)]^2\right)\right). $$ The population eigenvector is fixed here; its sign does not affect the expression.

The empirical bootstrap consistently estimates this limit. More precisely, conditionally on the data, the distribution of \(\sqrt n(\hat\lambda_j^*-\hat\lambda_j)\) converges in probability to the same normal distribution. This statement assumes independent observations, fixed \(p\), the positive population eigengap, and the finite fourth moment above.

Result and proof: Under those conditions, the ordinary empirical bootstrap is consistent for the centered, scaled simple eigenvalue. The covariance central limit theorem applies to the vector of sample covariance entries: its first-order contribution from observation \(X_i\) is \((X_i-\mu)(X_i-\mu)^\top-\Sigma\). Finite fourth moments ensure finite variances for these entries. The empirical bootstrap version of this central limit theorem holds conditionally in probability, with the same limiting Gaussian covariance. Replacing \(\mu\) by \(\bar X\) in the sample covariance changes the first-order expansion by \(o_p(1)\), and the analogous conditional remainder for \(S_n^*\) is also \(o_p(1)\). Finally, the perturbation bound above applies with probability tending to one: the sample eigengap converges to the positive population gap, while \(\|S_n-\Sigma\|=O_p(n^{-1/2})\) and \(\|S_n^*-S_n\|=O_p^*(n^{-1/2})\). After multiplication by \(\sqrt n\), the quadratic remainders vanish. The linear terms in both expansions therefore have the same limiting law, proving the claim. \(\square\)

Worked Example: A Separated Gaussian Eigenvalue

Let \(X_i\sim N(\mu,\Sigma)\) in two dimensions with \(\Sigma=\operatorname{diag}(4,1)\). Find the asymptotic variance of the largest sample covariance eigenvalue and describe its bootstrap approximation.

The largest population eigenvalue is \(4\), its unit eigenvector is \(v_1=(1,0)^\top\), and its eigengap is \(4-1=3\). The scalar projection \(v_1^\top(X_i-\mu)\) is normal with variance \(4\). If \(Z\sim N(0,1)\), its square is distributed as \(4Z^2\), so $$ \operatorname{Var}\left([v_1^\top(X_i-\mu)]^2\right) =\operatorname{Var}(4Z^2) =16\operatorname{Var}(Z^2) =16(2)=32. $$ Consequently, \(\sqrt n(\hat\lambda_1-4)\xrightarrow{d}N(0,32)\). Resampling the observed vectors and centering each bootstrap root at \(\hat\lambda_1\) consistently estimates this distribution. The gap makes the eigenvalue's first-order fluctuation the variance fluctuation in the fixed direction \(v_1\).

Repeated Eigenvalues: A Matrix, Not a Scalar, Drives the Limit

When an eigenvalue \(\lambda\) has multiplicity \(r>1\), there is no unique population eigenvector for it. Let the columns of \(U\) be an orthonormal basis for its eigenspace. The first-order changes of the \(r\) ordered roots in that cluster are the ordered eigenvalues of the \(r\times r\) projected perturbation \(U^\top E U\), not the diagonal entries of a perturbation in one arbitrarily chosen direction.

This follows by expressing the perturbed matrix in coordinates for the eigenspace and its orthogonal complement. Its upper-left block is \(\lambda I_r+U^\top E U\); the off-diagonal blocks are order \(\|E\|\). Because the rest of the spectrum is separated from \(\lambda\), eliminating the complement changes the upper-left block only by order \(\|E\|^2/\gamma\), where \(\gamma\) is the gap from the cluster to the rest of the spectrum. Thus the cluster roots, to first order, are \(\lambda\) plus the ordered eigenvalues of \(U^\top E U\).

For the sample covariance, the projected first-order perturbation has a Gaussian limit \(G\). The cluster's scaled root fluctuations therefore have the distribution of the ordered eigenvalues of \(G\). But bootstrap roots are centered at the already perturbed sample roots. Their conditional first-order behavior involves the eigenvalues of \(D_n+G_n^*\), centered by those of \(D_n\), where $$ D_n=\sqrt n\,U^\top(S_n-\Sigma)U, \qquad G_n^*=\sqrt n\,U^\top(S_n^*-S_n)U. $$ The matrix \(D_n\) does not converge to zero: it has a nondegenerate Gaussian limit. As a result, the conditional bootstrap distribution generally depends on the random sample perturbation, rather than converging to the fixed target distribution of the ordered eigenvalues of \(G\).

Worked Example: Two Equal Population Eigenvalues

Take \(X_i\sim N(0,I_2)\). Both population eigenvalues equal \(1\), so the first-order perturbation for the pair is a full symmetric Gaussian matrix. In the limit write it as $$ G=\begin{pmatrix}a&b\\b&d\end{pmatrix}, $$ where \(a,d\) are independent \(N(0,2)\), \(b\sim N(0,1)\), and these three entries are independent. The largest-root fluctuation has limit $$ \lambda_{\max}(G)=\frac{a+d}{2}+ \sqrt{\left(\frac{a-d}{2}\right)^2+b^2}. $$ It is not just the fluctuation of the first diagonal entry: the random off-diagonal entry also changes the root at first order.

To see what centering does to the bootstrap, suppose the realized local sample perturbation is \(D=\operatorname{diag}(c,-c)\), with \(c\ge0\), and write an independent bootstrap perturbation in the limit as \(G^*\), with entries \(a^*,b^*,d^*\) having the same distributions as \(a,b,d\). The largest bootstrap root, centered at the sample largest root, has conditional first-order form $$ \frac{a^*+d^*}{2} +\sqrt{\left(c+\frac{a^*-d^*}{2}\right)^2+(b^*)^2}-c. $$ When \(c=0\), this is \(\lambda_{\max}(G^*)\), the target form. For \(c>0\), the expression inside the square root is shifted by \(c\), so the conditional distribution depends on the observed local perturbation. Since the actual \(D_n\) has a nondegenerate limit, this dependence does not disappear as \(n\) grows. This illustrates why the ordinary centered bootstrap is generally inconsistent at a repeated root.

Nearly Repeated Roots and the Scale of the Eigengap

For a simple root with eigengap \(\gamma_n\), the perturbation remainder is of order \(\|E\|^2/\gamma_n\). Sample covariance perturbations are typically order \(n^{-1/2}\), so their remainder after scaling by \(\sqrt n\) is of order \(1/(\sqrt n\,\gamma_n)\). The usual scalar, linear approximation is justified when \(\sqrt n\,\gamma_n\to\infty\). If the gap is of order \(n^{-1/2}\), the neighboring roots interact at the same scale as sampling noise; the projected matrix description, rather than a fixed eigenvector derivative, is needed.

Worked Example: A Gap at the Sampling Scale

Consider a sequence of two-dimensional covariance matrices with eigenvalues \(1+4/\sqrt n\) and \(1\). Their gap tends to zero, but \(\sqrt n\) times the gap is \(4\), not infinity. In coordinates aligned with these population eigenvectors, the local population matrix is \(\operatorname{diag}(4,0)\). A first-order sample perturbation \(H=\begin{pmatrix}1&2\\2&-1\end{pmatrix}\) changes the local matrix to $$ \begin{pmatrix}4&0\\0&0\end{pmatrix}+H = \begin{pmatrix}5&2\\2&-1\end{pmatrix}. $$ Its largest eigenvalue is $$ \frac{5+(-1)}{2} +\sqrt{\left(\frac{5-(-1)}{2}\right)^2+2^2} =2+\sqrt{13}. $$ The local sample root is therefore affected by the off-diagonal perturbation. Treating the top population eigenvector as fixed and using only the top-left entry of \(H\) would give \(4+1=5\), which is different from \(2+\sqrt{13}\). The fixed-gap scalar approximation cannot be assumed in this regime.

By contrast, if the gap is \(n^{-1/4}\), then \(\sqrt n\) times the gap diverges. The neighboring root is increasingly well separated relative to the sampling perturbation, and the simple-root approximation can apply. The relevant comparison is not only whether a gap is numerically small, but how it compares with the \(n^{-1/2}\) sampling scale.

Practical Interpretation and a Common Mistake

For a clearly isolated covariance eigenvalue, a bootstrap histogram of \(\sqrt n(\hat\lambda_j^*-\hat\lambda_j)\) is a principled estimate of its sampling distribution under the conditions stated above. For repeated or near-repeated roots, a histogram may still be computable, but computation alone does not establish that it estimates the desired sampling law. The bootstrap's centering and the random interaction among roots matter.

A common mistake is to apply the isolated-root formula at a repeated eigenvalue by selecting any one eigenvector and using \(v^\top(S_n-\Sigma)v\). The correction is that the eigenspace, not an arbitrary basis vector, determines the first-order behavior: project the perturbation onto the whole repeated-root subspace and examine that projected matrix's eigenvalues. Near repetition, check whether the eigengap is large relative to \(n^{-1/2}\); if not, the simple-root bootstrap conclusion is not justified.

Key takeaway: A positive fixed eigengap gives a smooth, scalar first-order fluctuation and a consistent ordinary bootstrap for an isolated covariance eigenvalue. Repeated roots have matrix-valued first-order fluctuations, while gaps on the \(n^{-1/2}\) scale lead to local interactions that can make the centered bootstrap distribution depend on the observed sample perturbation.

Check Your Understanding

Use the definitions and arguments above to answer the following questions.

  1. What quantity determines the asymptotic variance of a simple sample covariance eigenvalue?
  2. Why does a repeated eigenvalue require projecting a perturbation onto its full eigenspace?
  3. For a gap \(\gamma_n\), what condition makes the scaled quadratic perturbation remainder negligible?
  4. In the two-dimensional repeated-root example, which random matrix entry would be missed by considering only a fixed coordinate direction?
  5. Why does centering bootstrap roots at the sample roots create a difficulty at a repeated population eigenvalue?
Clinical Biostats