Parametric Versus Nonparametric PCA Bootstrap

Two Ways to Represent Sampling Uncertainty

A PCA bootstrap approximates how estimated eigenvalues and eigenspaces would vary across repeated samples. The central choice is what a bootstrap sample is assumed to come from. A parametric bootstrap generates new observations from a fitted distribution, commonly a multivariate normal distribution. A nonparametric bootstrap draws observed data vectors with replacement, allowing the empirical distribution to stand in for the unknown population distribution.

The methods can agree when the Gaussian model describes the relevant features of the data, but they need not. In particular, Gaussian resampling determines fourth moments from the fitted covariance, whereas observation resampling retains the observed extremes, skewness, and other features represented in the sample. Since PCA uncertainty depends on variability in quadratic and cross-product terms, this difference can materially affect both eigenvalue and eigenspace uncertainty.

Definition: Let $X_1,\ldots,X_n$ be independent, identically distributed $p$-dimensional observations with mean $\mu$ and covariance matrix $\Sigma$. Define the sample mean and unbiased sample covariance by $$ \bar X=\frac{1}{n}\sum_{i=1}^n X_i,\qquad S=\frac{1}{n-1}\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^T. $$ A nonparametric bootstrap sample draws $X_1^*,\ldots,X_n^*$ independently from the empirical distribution that places probability $1/n$ on each observed vector. A Gaussian parametric bootstrap draws new observations independently from $N(\bar X,S)$. In either case, recompute the mean, covariance, eigenvalues, and eigenspaces in each replicate. The superscript $*$ denotes a bootstrap quantity.

The fitted Gaussian distribution above uses $S$ as its covariance. If a software implementation instead fits the maximum-likelihood covariance with divisor $n$, its finite-sample simulation scale differs. A comparison should therefore document the fitting convention as well as the resampling method.

What Each Bootstrap Preserves

The Gaussian parametric bootstrap preserves the estimated mean and covariance, and imposes multivariate normality. It does not preserve the empirical distribution’s specific tail behavior or dependence patterns beyond those encoded by $S$. The nonparametric bootstrap preserves the observed vectors as its support: each resampled vector retains all its coordinates together, so the observed dependence among variables is retained. It does not generate values outside the observed sample, and its approximation can be poor when the sample is too small to represent important population features.

Both methods require independent sampling for the ordinary versions described here. They also presume that the observations are representative of the population relevant to the inference. If the data are clustered or serially dependent, resampling individual rows independently does not respect that sampling structure; an appropriate dependence-aware scheme is needed instead.

Worked Example: A Finite-Sample Covariance Check

Consider four centered vectors in two dimensions: $(2,0)$, $(-2,0)$, $(0,1)$, and $(0,-1)$. Their sample mean is $(0,0)$. Their covariance is $$ S=\frac{1}{3} \begin{pmatrix} 8&0\\ 0&2 \end{pmatrix} = \begin{pmatrix} 8/3&0\\ 0&2/3 \end{pmatrix}. $$

For a nonparametric bootstrap sample, the four draws are independent from the empirical distribution on these vectors. Its conditional covariance as a single draw is $$ \frac{1}{4}\sum_{i=1}^4 X_iX_i^T = \begin{pmatrix} 2&0\\ 0&1/2 \end{pmatrix} =\frac{3}{4}S. $$ The expected unbiased sample covariance across bootstrap replicates is this same matrix, as proved below. By contrast, the fitted Gaussian distribution $N((0,0),S)$ has covariance $S$, so its simulated unbiased sample covariance has expectation $S$. The differing finite-sample expectations are not a claim that the bootstrap eigenvalues themselves have these expectations: eigenvalues are nonlinear functions of the covariance matrix.

Two Results That Explain the Comparison

Result 1: Conditional covariance of the nonparametric bootstrap. Given the observed data, the expected unbiased sample covariance of a nonparametric bootstrap sample is $(n-1)S/n$. For a parametric sample of size $n$ from a distribution with covariance $S$, the expected unbiased sample covariance is $S$.

Proof. Conditional on the data, each $X_i^*$ has mean $\bar X$ and covariance $$ C_n=\frac{1}{n}\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^T =\frac{n-1}{n}S. $$ For any independent sample $Y_1,\ldots,Y_n$ with a common finite covariance $C$, the identity $$ \sum_{i=1}^n(Y_i-\bar Y)(Y_i-\bar Y)^T = \sum_{i=1}^n(Y_i-\mu_Y)(Y_i-\mu_Y)^T -n(\bar Y-\mu_Y)(\bar Y-\mu_Y)^T $$ gives an expectation of $nC-C=(n-1)C$. Dividing by $n-1$ shows that its unbiased sample covariance has expectation $C$. Apply this conditionally with $C=C_n$ for the nonparametric bootstrap, and with $C=S$ for the Gaussian parametric bootstrap. $\square$

The next result connects covariance perturbations to eigenvalue uncertainty. Let $\lambda_j$ be a simple eigenvalue of a symmetric matrix $A$, with unit eigenvector $u_j$, and let $\delta=\min_{k\ne j}|\lambda_j-\lambda_k|>0$. If $E$ is symmetric and $\|E\|<\delta/4$, the corresponding eigenvalue $\tilde\lambda_j$ of $A+E$ satisfies $$ \tilde\lambda_j-\lambda_j=u_j^TEu_j+R,\qquad |R|\leq \frac{2\|E\|^2}{\delta}. $$ Here $\|\cdot\|$ is the spectral norm. The gap condition ensures that the eigenvalue remains identifiable as the one corresponding to $\lambda_j$.

Proof. Express $A+E$ in an orthonormal eigenbasis of $A$, placing $u_j$ first. In this basis it has block form $$ \begin{pmatrix} \lambda_j+e & b^T\\ b & D+F \end{pmatrix}, $$ where $e=u_j^TEu_j$, $\|b\|\leq\|E\|$, $\|F\|\leq\|E\|$, and the eigenvalues of $D$ are at least $\delta$ away from $\lambda_j$. By the eigenvalue perturbation bound $|\tilde\lambda_j-\lambda_j|\leq\|E\|$. Thus the distance from $\tilde\lambda_j$ to the spectrum of $D+F$ is at least $\delta-2\|E\|>\delta/2$. The eigenvector block equations, eliminating the lower block, yield $$ \tilde\lambda_j=\lambda_j+e+ b^T\{\tilde\lambda_j I-(D+F)\}^{-1}b. $$ The inverse has norm at most $2/\delta$, so the final term has absolute value at most $2\|E\|^2/\delta$. This proves the expansion and its remainder bound. $\square$

For PCA, take $A=\Sigma$ and $E=S-\Sigma$. To first order, the eigenvalue error is $u_j^T(S-\Sigma)u_j$. With $Y=u_j^T(X-\mu)$, this quantity is driven by the centered squared projected observations $Y^2-E[Y^2]$; estimating its variability therefore involves a fourth moment of the data. A Gaussian parametric bootstrap imposes the Gaussian value of that fourth moment. A nonparametric bootstrap estimates it from the observed projected scores. This is one reason the two procedures can give different uncertainty estimates even when their fitted covariance matrices are the same.

Worked Example: Comparing Projected-Score Variability

Suppose the observed one-dimensional centered scores are $-2,-1,0,1,2$. Their mean is zero, and their unbiased sample variance is $$ \hat\lambda=\frac{4+1+0+1+4}{4}=\frac{5}{2}. $$ The empirical mean of their squared values is $2$. For a first-order, nonparametric estimate of eigenvalue variance, the empirical variance of the squared scores, divided by $n=5$, is used. The squared scores are $4,1,0,1,4$, so their empirical variance with divisor $5$ is $$ \frac{(4-2)^2+(1-2)^2+(0-2)^2+(1-2)^2+(4-2)^2}{5} =\frac{14}{5}. $$ The resulting first-order estimate is $14/25=0.56$. Under a fitted Gaussian model with variance $\hat\lambda=5/2$, the Gaussian variance of a squared centered score is $2\hat\lambda^2=25/2$, giving the parametric first-order estimate $(25/2)/5=2.5$. These are different approximations based on different fourth-moment assumptions, not exact finite-sample bootstrap variances. The example also shows why matching the estimated variance does not force the two uncertainty estimates to match.

Eigenspaces, Not Signed Eigenvectors

For a leading $r$-dimensional eigenspace, let $U$ be a matrix whose columns are orthonormal eigenvectors spanning that space. Its orthogonal projector is $P=UU^T$. Comparing $P^*$ with $P$, for example through $\|P^*-P\|_F$, avoids treating a sign change in an eigenvector as a change in the underlying direction. For a single direction, the same principle means comparing sign-invariant quantities or aligning signs before summarizing replicates.

Both bootstraps should be evaluated using the same eigenspace summary and the same confidence-region construction. A small gap at the boundary of the leading subspace can make its orientation unstable; it is then important to distinguish uncertainty about individual eigenvectors from uncertainty about the whole subspace. Bootstrap replicates that exchange nearly tied components should not be interpreted as ordinary sign flips. The earlier material on bootstrap failure near repeated eigenvalues and eigen-gap diagnostics is relevant to that distinction.

Worked Example: A Directional Uncertainty Calculation

Suppose a two-dimensional population covariance is $\Sigma=\operatorname{diag}(4,1)$. Its leading eigenvector is $u_1=(1,0)^T$, its second eigenvector is $u_2=(0,1)^T$, and the eigenvalue gap is $4-1=3$. The first-order rotation of the leading direction is driven by the off-diagonal sample covariance error, divided by the gap.

If the observations are Gaussian, the variance of the cross-product $(X_1-\mu_1)(X_2-\mu_2)$ is $4\cdot1=4$, because the coordinates are independent and their variances are 4 and 1. Thus the first-order variance of the direction’s angular error is $$ \frac{4}{n(3)^2}=\frac{4}{9n}. $$ For non-Gaussian observations with the same covariance, the cross-product variance need not be 4. For example, let the coordinates be independent and take values $\pm2$ and $\pm1$, each sign with probability $1/2$. The cross-product is then always either $2$ or $-2$, so its mean is zero and its variance is 4. In this particular distribution the first-order variance matches the Gaussian calculation. Other dependence or tail behavior can change that variance, which is precisely the feature that observation resampling can reflect when it is represented in the sample.

Practical Comparison and Common Pitfalls

A useful comparison runs both procedures with the same sample size, covariance convention, number of replicates, and inferential summary. For each method, examine eigenvalue intervals or replicate quantiles, and examine eigenspace distances or confidence regions. Also inspect the bootstrap distributions themselves: skewness, long tails, abrupt changes in component ordering, and a substantial fraction of poorly separated replicate eigenvalues can matter even when two standard errors look similar.

Disagreement is diagnostic, not an automatic vote for one method. If the Gaussian assumption is defensible and the sample is modest, the parametric bootstrap can provide a more stable model-based approximation. If distributional shape beyond covariance is uncertain and the sample adequately represents that shape, observation resampling avoids imposing Gaussian fourth moments. Conversely, the nonparametric bootstrap cannot reliably reproduce rare extremes that were not observed. Neither method repairs a weak eigen-gap or an unrepresentative sample.

Common mistake: treating equal fitted covariance matrices as equal bootstrap uncertainty. The first-order eigenvalue fluctuation depends on squared projected scores, and directional fluctuation depends on cross-products. Their variability is not determined by covariance alone outside the Gaussian model.
Correction: Compare the relevant projected-score and cross-product variability, as well as the resulting bootstrap distributions.
Common mistake: comparing signed eigenvectors directly. A sign reversal represents the same one-dimensional eigenspace.
Correction: Use projectors, sign-invariant comparisons, or principal angles, while checking eigenvalue gaps and component ordering.
Common mistake: reading the nonparametric bootstrap’s finite-sample covariance correction as evidence of a different population covariance. Its conditional covariance calculation reflects resampling from the empirical distribution and recentering the replicate.
Correction: Keep the resampling mechanism, covariance divisor, and target of inference explicit when interpreting results.

Key takeaway: Gaussian parametric and observation-resampling PCA bootstraps answer uncertainty questions under different distributional representations. Compare their assumptions, eigenvalue gaps, projected-score variability, and sign-invariant eigenspace summaries; agreement is reassuring, while disagreement identifies sensitivity that should be explained rather than hidden.

Check Your Understanding

Use the definitions and results above to answer these questions.

What distributional feature does a Gaussian parametric bootstrap impose that an observation-resampling bootstrap estimates from the observed data?
Why is the expected unbiased covariance of a nonparametric bootstrap sample $(n-1)S/n$, conditional on the data?
For a simple eigenvalue, which projected quantity determines its first-order covariance perturbation?
Why is a projector a suitable summary for comparing bootstrap eigenspaces?
Name two reasons that agreement between the two bootstrap procedures could fail to guarantee reliable inference.