Bootstrap Confidence Intervals for Eigenvalues

Why Bootstrap an Eigenvalue?

A confidence interval for a population covariance eigenvalue must account for how the entire sample covariance matrix varies. An eigenvalue is not an ordinary coordinate of that matrix: its uncertainty depends on the corresponding population eigenvector and on the distribution of observations along that direction. A bootstrap approximates this uncertainty by repeatedly generating samples and recalculating the ordered eigenvalue.

There are two common ways to generate those samples. A nonparametric bootstrap resamples observed data rows, preserving their empirical multivariate distribution. A parametric bootstrap generates observations from a fitted probability model. Either method can produce intervals, but their validity depends on the sampling assumptions and, for ordinary eigenvalue intervals, on the population eigenvalue being separated from the others.

Definition: Let $X_1,\ldots,X_n$ be independent and identically distributed $p$-dimensional observations with mean $\mu$ and covariance matrix $\Sigma$. Write the eigenvalues of $\Sigma$ in descending order as $\lambda_1\geq\cdots\geq\lambda_p$, and those of the sample covariance matrix $$ S=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar X)(X_i-\bar X)^T $$ as $\hat\lambda_1\geq\cdots\geq\hat\lambda_p$. A bootstrap replicate is a new sample generated conditionally on the observed data; its covariance eigenvalue of rank $j$ is denoted $\hat\lambda_j^*$.

Two Ways to Generate Bootstrap Samples

For a nonparametric bootstrap, draw $n$ rows independently with replacement from the observed rows $X_1,\ldots,X_n$. The resulting sample covariance is $S^*$; sort its eigenvalues in descending order and record $\hat\lambda_j^*$. Repeat this procedure $B$ times. Resampling whole rows is essential: it preserves the observed dependence among coordinates.

For a parametric bootstrap, specify a distributional model and estimate its parameters. In the multivariate Gaussian case, generate each bootstrap row independently from $N_p(\bar X,S)$, then calculate $S^*$ and its ordered eigenvalues. More generally, use the fitted model in place of the Gaussian distribution. A parametric bootstrap is only as credible as that model; fitting a Gaussian distribution does not make a non-Gaussian sampling distribution Gaussian.

In either method, sort eigenvalues within every replicate and use the same rank $j$ throughout. Do not identify an eigenvalue by matching its eigenvector: eigenvector signs are arbitrary, and nearly equal eigenvalues can make individual directions unstable. Here the target is the ordered eigenvalue itself.

Fit and compute.
Calculate $S$, its ordered eigenvalues, and any standard errors needed for studentization.

Resample.
Generate $B$ samples using either the empirical distribution of rows or a fitted parametric model.

Recalculate.
For each replicate, form $S^*$, sort its eigenvalues, and store the rank-$j$ value.

Form and report the interval.
Use the bootstrap distribution to obtain endpoints, and state the resampling scheme and number of replicates.

Percentile, Basic, and Studentized Intervals

Let $\hat\lambda_{j,(\alpha)}^*$ denote the empirical $\alpha$-quantile of the $B$ bootstrap values for rank $j$. The percentile interval takes the corresponding quantiles directly:

\left[\hat\lambda_{j,(\alpha/2)}^*,\ \hat\lambda_{j,(1-\alpha/2)}^*\right].

For a nominal $100(1-\alpha)\%$ interval, the lower and upper quantiles are usually $\alpha/2$ and $1-\alpha/2$. The basic interval reflects the bootstrap errors around the observed estimate:

\left[2\hat\lambda_j-\hat\lambda_{j,(1-\alpha/2)}^*,\ 2\hat\lambda_j-\hat\lambda_{j,(\alpha/2)}^*\right].

If replicate-specific standard errors are available, a studentized interval can account for changing scale. Define $T_j^*=(\hat\lambda_j^*-\hat\lambda_j)/\widehat{\operatorname{SE}}_j^*$, and let $q_\alpha^*$ be its empirical $\alpha$-quantile. With an observed-data standard error $\widehat{\operatorname{SE}}_j$, the interval is

\left[\hat\lambda_j-q_{1-\alpha/2}^*\widehat{\operatorname{SE}}_j,\ \hat\lambda_j-q_{\alpha/2}^*\widehat{\operatorname{SE}}_j\right].

Studentization can improve calibration when the estimated scale varies appreciably across samples, but it requires a reliable standard error in both the original and bootstrap samples. The percentile interval is simple to implement; comparing it with the basic or studentized interval is a useful sensitivity check.

Worked Example: Gaussian Parametric Bootstrap Approximation

Suppose $p=2$, the fitted covariance eigenvalues are $4.2$ and $1.1$, and $n=100$. Construct a first-order, normal approximation to a 95% Gaussian parametric bootstrap interval for the leading eigenvalue.

For a Gaussian population with a simple eigenvalue $\lambda_j$, the asymptotic variance of $\sqrt n(\hat\lambda_j-\lambda_j)$ is $2\lambda_j^2$. Substituting the fitted leading eigenvalue $4.2$ gives the estimated standard error $$ \widehat{\operatorname{SE}}_1\approx \frac{\sqrt{2}(4.2)}{\sqrt{100}} =\frac{5.9397}{10}=0.5940. $$

The normal critical value for a two-sided 95% interval is $1.96$. Thus the first-order interval is $$ 4.2\pm 1.96(0.5940)=4.2\pm 1.1642=[3.0358,\ 5.3642]. $$

This is a normal approximation to the fitted Gaussian bootstrap distribution, not a claim that finite-$B$ bootstrap quantiles equal these numbers. In software, the parametric bootstrap interval is obtained from the simulated eigenvalue quantiles; the calculation here supplies a transparent benchmark for those results.

Why a Simple Eigenvalue Has a Bootstrap Limit

A useful condition for ordinary scalar eigenvalue intervals is that the target population eigenvalue be simple: $\lambda_j$ differs from every other eigenvalue. Let $v_j$ be a corresponding unit eigenvector. The first-order change in this eigenvalue under a small symmetric perturbation $H$ of $\Sigma$ is $v_j^THv_j$. Consequently, the influence of one observation on the eigenvalue is determined by its squared projection onto $v_j$.

Theorem: Suppose $p$ is fixed, the observations are iid, $E\|X_i\|^{4+\delta}<\infty$ for some $\delta>0$, and $\lambda_j$ is a simple eigenvalue of $\Sigma$. Then $$ \sqrt n(\hat\lambda_j-\lambda_j)\xrightarrow{d}N(0,\tau_j^2), \qquad \tau_j^2=\operatorname{Var}\!\left(\{v_j^T(X-\mu)\}^2\right). $$ The nonparametric bootstrap consistently estimates this limiting distribution: conditionally on the observed data, the distribution of $\sqrt n(\hat\lambda_j^*-\hat\lambda_j)$ converges in probability to $N(0,\tau_j^2)$. A parametric bootstrap has this property when its fitted model is correctly specified and its parameter estimates consistently estimate the model parameters.

Proof. The sample covariance satisfies a fixed-dimensional central limit theorem under the stated moment condition. Since estimating the mean contributes only a term of smaller order to the first-order covariance fluctuation, its leading perturbation is the centered average of $(X_i-\mu)(X_i-\mu)^T-\Sigma$. The eigenvalue derivative in direction $H$ is $v_j^THv_j$. Applying this derivative to the covariance fluctuation gives the scalar summands $$ v_j^T\big((X_i-\mu)(X_i-\mu)^T-\Sigma\big)v_j =\{v_j^T(X_i-\mu)\}^2-\lambda_j. $$ Their mean is zero and their variance is $\tau_j^2$. The scalar central limit theorem therefore gives the displayed normal limit.

For the nonparametric bootstrap, the empirical distribution has finite fourth moments with probability one for each observed sample. The conditional bootstrap central limit theorem for the covariance entries gives the same first-order covariance fluctuation, centered at $S$, with its empirical distribution replacing the population distribution. The empirical variance of the projected squared observations converges in probability to $\tau_j^2$. Because the eigenvalue is simple, its derivative varies continuously near $\Sigma$; the eigenvalue expansion has a remainder of smaller order than $n^{-1/2}$. Applying that expansion to both the original and resampled covariances proves the stated conditional limit. For a correctly specified parametric model, the corresponding conditional limit follows from the model’s central limit theorem and consistent parameter estimation. $\square$

If the limiting distribution is continuous and $\tau_j^2>0$, convergence of the centered bootstrap quantiles implies first-order coverage for the percentile interval. In particular, its lower and upper endpoints are approximately $\hat\lambda_j+z_{\alpha/2}\tau_j/\sqrt n$ and $\hat\lambda_j+z_{1-\alpha/2}\tau_j/\sqrt n$, where $z_a$ is the $a$-quantile of $N(0,1)$. The asymptotic distribution of the estimation error then gives coverage $1-\alpha$. The same first-order conclusion holds for the basic interval. The moment and gap conditions matter: this result does not establish validity at a repeated population eigenvalue.

Worked Example: A Nonparametric Bootstrap Replicate

Take four observed vectors $x_1=(2,0)^T$, $x_2=(-2,0)^T$, $x_3=(0,1)^T$, and $x_4=(0,-1)^T$. Their mean is zero, so their sample covariance is $$ S=\frac{1}{3} \begin{pmatrix}8&0\\0&2\end{pmatrix} = \begin{pmatrix}8/3&0\\0&2/3\end{pmatrix}. $$ The leading sample eigenvalue is $8/3$.

Suppose one bootstrap draw contains $x_1,x_1,x_2,x_4$. Its mean is $(1/2,-1/4)^T$. Centering the four rows and dividing their cross-product sum by $n-1=3$ gives $$ S^*= \begin{pmatrix}11/3&1/6\\1/6&1/4\end{pmatrix}. $$ The trace is $47/12$, and the difference between the diagonal entries is $41/12$. The larger eigenvalue is therefore $$ \hat\lambda_1^*= \frac{1}{2}\left(\frac{47}{12}+ \sqrt{\left(\frac{41}{12}\right)^2+4\left(\frac{1}{6}\right)^2}\right) =\frac{47+\sqrt{1697}}{24}\approx 3.675. $$ This single replicate contributes one value to the bootstrap distribution; it is not itself a confidence interval. Repeating the resampling and covariance calculation supplies the quantiles needed for an interval.

Comparing Bootstrap and Asymptotic Intervals

An asymptotic Wald interval uses an estimate of $\tau_j$ in $\hat\lambda_j\pm z_{1-\alpha/2}\widehat{\tau}_j/\sqrt n$. A practical nonparametric estimate uses the sample variance of squared projected scores. If $\hat v_j$ is a unit eigenvector of $S$, set $z_i=\hat v_j^T(X_i-\bar X)$; then estimate $\tau_j^2$ by the sample variance of $z_1^2,\ldots,z_n^2$. Under Gaussian sampling, $\tau_j^2=2\lambda_j^2$, giving the familiar plug-in standard error $\sqrt{2}\hat\lambda_j/\sqrt n$.

The asymptotic interval is fast and provides a useful benchmark. Bootstrap intervals can reflect skewness and other features of the estimated sampling distribution, but they do not automatically fix a poor model or a weak approximation. A meaningful comparison should consider estimated coverage and interval length, not just whether the endpoints differ.

Worked Example: Non-Gaussian Projection and Interval Width

Let $X=(Y,Z)^T$, where $Y=0$ with probability $1/2$, $Y=\sqrt 8$ or $-\sqrt 8$ with probability $1/4$ each, and $Z\sim N(0,1)$ is independent of $Y$. Then $E[Y]=0$, $E[Y^2]=4$, and $E[Y^4]=32$. The covariance eigenvalues are $4$ and $1$, so the leading eigenvalue is simple. Its asymptotic variance is $$ \tau_1^2=\operatorname{Var}(Y^2)=E[Y^4]-\{E[Y^2]\}^2=32-16=16. $$

At $n=100$, the leading-eigenvalue standard error is approximately $4/\sqrt{100}=0.4$. If the observed estimate is $4$, the corresponding first-order 95% interval is $4\pm1.96(0.4)=[3.216,4.784]$. A Gaussian parametric calculation instead uses $\sqrt{2}(4)/10\approx0.566$, giving $[2.891,5.109]$. The difference arises because the actual projection has a different fourth moment from a Gaussian projection. A row-resampling bootstrap can reflect the empirical projected fourth moment; a Gaussian parametric bootstrap imposes Gaussian fourth moments. Which interval is appropriate depends on the plausibility of the model and the adequacy of the sample.

Assessing Coverage in Practice

Coverage is a repeated-sampling property: a nominal 95% procedure should contain the fixed population eigenvalue in about 95% of repeated samples under its assumptions. To compare procedures, simulate many independent data sets from a specified distribution with known $\lambda_j$. For each data set, construct the asymptotic, nonparametric bootstrap, and parametric bootstrap intervals, then record whether each interval contains that known value.

With $K$ simulated data sets, estimate coverage by the fraction of intervals that contain $\lambda_j$. If the estimated coverage is $\hat c$, its approximate Monte Carlo standard error is $\sqrt{\hat c(1-\hat c)/K}$. For example, an estimated coverage of $0.95$ from $K=1000$ runs has Monte Carlo standard error $\sqrt{0.95(0.05)/1000}\approx0.0069$. This simulation error should be distinguished from the procedure’s actual coverage error. Report interval length as well: a procedure can obtain high coverage simply by producing unnecessarily wide intervals.

Common pitfall: Resampling rows does not guarantee a valid interval for every ordered eigenvalue in every setting. The ordinary first-order result above uses a simple population eigenvalue and fixed dimension; it does not justify treating a repeated eigenvalue as a smooth scalar target. Also, a Gaussian parametric bootstrap is not distribution-free: if the observations are non-Gaussian, its implied fourth moments can give the wrong eigenvalue variance. State the sampling model and eigenvalue-separation assumption, and compare against an appropriate asymptotic benchmark.

Check Your Understanding

Answer each question by identifying the relevant resampling or inferential step.

How does a nonparametric bootstrap sample differ from a Gaussian parametric bootstrap sample?
Why must the eigenvalues be sorted within each bootstrap replicate before recording the rank-$j$ value?
For a simple eigenvalue, what projected observation determines its first-order influence, and what moment controls its asymptotic variance?
What are the endpoint formulas for the percentile and basic bootstrap intervals?
Why can a Gaussian parametric bootstrap give a different interval from a row-resampling bootstrap under non-Gaussian sampling?
In a coverage simulation, what does the Monte Carlo standard error quantify, and why should interval length also be reported?