Asymptotic Uncertainty for SAVE Eigenspaces

Why SAVE Eigenspaces Have Sampling Uncertainty

SAVE identifies directions through changes in the conditional covariance of predictors across response slices. In a sample, both the frequency of each slice and its conditional covariance matrix vary. Those fluctuations perturb the SAVE operator, and the operator’s eigenvalue gaps determine how much its estimated directions can move. This tutorial derives that first-order uncertainty for fixed slices and shows how it carries through to an estimated eigenspace.

We work with predictors that have already been standardized using their population mean and covariance. This isolates the sampling variation from the slice estimator itself. In applications where the standardization parameters are estimated, their additional first-order variation must also be included; the formulas below do not silently treat that variation as zero. We take the slice boundaries as fixed, with positive population probability in every slice.

Definition: Let $Z\in\mathbb{R}^p$ satisfy $E[Z]=0$ and $\operatorname{Cov}(Z)=I_p$, and let $H$ denote the slice containing the response. For slice $h$, define $p_h=P(H=h)>0$, $m_h=E[Z\mid H=h]$, and $C_h=\operatorname{Cov}(Z\mid H=h)$. The population sliced SAVE operator is $$ M=\sum_{h=1}^{H}p_h(I_p-C_h)^2. $$ For observations $Z_1,\ldots,Z_n$, let $\hat p_h$ be the observed slice proportion and $\hat C_h$ the within-slice sample covariance, using the slice sample mean. The sample SAVE operator is $$ \hat M=\sum_{h=1}^{H}\hat p_h(I_p-\hat C_h)^2. $$

The square in these definitions is matrix multiplication. Since each $C_h$ is symmetric, $M$ and $\hat M$ are symmetric. Their eigenspaces can therefore be analyzed with ordinary symmetric-matrix perturbation arguments. The relevant randomness in this fixed-slice setup comes from the observations falling in each slice and from the predictor values within those slices.

First-Order Fluctuation of the SAVE Operator

Write $I_h=1\{H=h\}$, and for an observation $Z$ in slice $h$ define $$ Q_h=(Z-m_h)(Z-m_h)^T,\qquad A_h=(I_p-C_h)^2. $$ The first-order fluctuation of $\hat p_h$ is represented by $I_h-p_h$. The corresponding fluctuation of $\hat C_h$ is represented by $I_h(Q_h-C_h)/p_h$. The factor $1/p_h$ matters: estimating a conditional covariance from a slice uses only the fraction of observations in that slice.

Theorem: Suppose the observations are independent and identically distributed, the number of fixed slices is finite, $p_h>0$ for every $h$, and $E[\|Z\|^4]<\infty$. Then $$ \sqrt{n}(\hat M-M)\xrightarrow{d}G, $$ where $G$ is a mean-zero, symmetric Gaussian matrix. Its distribution is determined by the influence matrix $$ \Psi(Z,H)=\sum_{h=1}^{H}\left\{(I_h-p_h)A_h-I_h\left[(I_p-C_h)(Q_h-C_h)+(Q_h-C_h)(I_p-C_h)\right]\right\}. $$ In particular, if $\operatorname{vec}$ stacks matrix columns, then $\operatorname{vec}(G)$ has covariance matrix $\operatorname{Var}(\operatorname{vec}\Psi)$.

Proof. For each slice, $\hat p_h=n^{-1}\sum_{i=1}^n1\{H_i=h\}$, so its observation-level influence is $I_h-p_h$. The conditional covariance is a ratio of within-slice moments. Expanding that ratio around $p_h$ and $m_h$, the first-order contribution from estimating the conditional mean cancels, leaving $$ \hat C_h-C_h=\frac{1}{n}\sum_{i=1}^n\frac{1\{H_i=h\}}{p_h}\left\{(Z_i-m_h)(Z_i-m_h)^T-C_h\right\}+o_p(n^{-1/2}). $$ The finite fourth moment and positive slice probabilities ensure the required sample moments obey the multivariate central limit theorem and the ratio expansion has a remainder of smaller order than $n^{-1/2}$.

For a symmetric perturbation $D$, the derivative of $A(C)=(I_p-C)^2$ at $C_h$ is $$ D A(C_h)[D]=-(I_p-C_h)D-D(I_p-C_h). $$ This follows by expanding $(I_p-C_h-tD)^2$, retaining the terms linear in $t$, and dividing by $t$. A first-order product expansion of $\hat p_h A(\hat C_h)$ therefore has the two contributions $$ (I_h-p_h)A_h \quad\text{and}\quad p_hD A(C_h)\left[\frac{I_h}{p_h}(Q_h-C_h)\right]. $$ Substituting the derivative and summing over slices gives $\Psi$. Its expectation is zero: $E[I_h-p_h]=0$, and $E[I_h(Q_h-C_h)]=p_hE[Q_h-C_h\mid H=h]=0$. The multivariate central limit theorem applied to the independent influence matrices now gives the stated Gaussian limit. $\square$

Worked Example: The Slice-Probability Contribution

Consider two slices with probabilities $p_1=p_2=1/2$. In a one-dimensional predictor problem, let the within-slice conditional variances be $C_1=0.4$ and $C_2=1.2$. The population SAVE operator is a scalar, and its value is $$ M=\frac12(1-0.4)^2+\frac12(1-1.2)^2 =\frac12(0.36)+\frac12(0.04)=0.20. $$

To isolate the effect of estimating the slice probabilities, hold the two conditional variances fixed. The two slice contributions are $A_1=0.36$ and $A_2=0.04$. Since $\hat p_2=1-\hat p_1$, the probability-only fluctuation is $$ (\hat p_1-p_1)A_1+(\hat p_2-p_2)A_2 =(\hat p_1-p_1)(A_1-A_2). $$ The variance of $\sqrt n(\hat p_1-p_1)$ is $p_1(1-p_1)=0.25$. Thus the asymptotic variance of this probability-only part is $$ 0.25(0.36-0.04)^2=0.25(0.1024)=0.0256. $$ This is not the full variance of $\hat M$: the within-slice covariance estimates also fluctuate and contribute through the second term of $\Psi$.

The theorem provides an asymptotic covariance for the operator, not yet for a selected direction. To obtain direction uncertainty, one must account for the separation between the relevant eigenvalues. A small gap allows modest operator noise to produce substantial eigenvector movement.

From Operator Fluctuation to Eigenspace Fluctuation

Let the eigenvalues of $M$ be ordered as $\lambda_1\geq\cdots\geq\lambda_p$, with corresponding orthonormal eigenvectors $u_1,\ldots,u_p$. For an individual eigenvector $u_j$, assume its eigenvalue is simple: $\lambda_j\neq\lambda_k$ for every $k\neq j$. Choose the sign of $\hat u_j$ so that $\hat u_j^Tu_j\geq0$.

Theorem: Under the operator limit above, if $\lambda_j$ is simple, then $$ \sqrt n(\hat u_j-u_j)\xrightarrow{d} \sum_{k\neq j}u_k\frac{u_k^TGu_j}{\lambda_j-\lambda_k}. $$ For the leading $d$-dimensional eigenspace, let $U=(u_1,\ldots,u_d)$, $P=UU^T$, and $\hat P=\hat U\hat U^T$. If the boundary eigenvalue gap satisfies $\lambda_d>\lambda_{d+1}$, then $$ \sqrt n(\hat P-P)\xrightarrow{d} \sum_{j=1}^{d}\sum_{k=d+1}^{p} \frac{u_k^TGu_j}{\lambda_j-\lambda_k} \left(u_ku_j^T+u_ju_k^T\right). $$

Proof. Put $E_n=\hat M-M$, so $\sqrt nE_n\xrightarrow{d}G$ and $E_n=O_p(n^{-1/2})$. For a simple eigenvalue, the eigenvector perturbation equation, projected onto $u_k$ for $k\neq j$, gives to first order $$ u_k^T(\hat u_j-u_j)=\frac{u_k^TE_nu_j}{\lambda_j-\lambda_k}+o_p(n^{-1/2}). $$ The component along $u_j$ is $o_p(n^{-1/2})$ under the sign convention and unit-length normalization. Summing the other components and multiplying by $\sqrt n$ yields the eigenvector limit.

For the projector, differentiate $P=UU^T$: a first-order change in a basis vector contributes both $u_k u_j^T$ and its transpose. Rotations within the target space change the basis $U$ but not $P$; hence only pairs with $j\leq d$ and $k>d$ appear. The gap $\lambda_d-\lambda_{d+1}>0$ separates the target space from its complement, so these denominators are nonzero and the projector has a stable first-order expansion. Substituting $\sqrt nE_n\xrightarrow{d}G$ proves the displayed limit. $\square$

Worked Example: A Three-Dimensional Leading Eigenspace

Suppose $M$ has eigenvalues $2.0,\ 1.2,\ 0.2$, with the leading one-dimensional eigenspace spanned by $u_1$. Its gaps from the other eigenvalues are $2.0-1.2=0.8$ and $2.0-0.2=1.8$. Suppose the limiting operator fluctuation $G$ has $$ \operatorname{Var}(u_2^TGu_1)=0.09,\qquad \operatorname{Var}(u_3^TGu_1)=0.04. $$ The two perpendicular components of the limiting eigenvector error have variances $$ \frac{0.09}{0.8^2}=\frac{0.09}{0.64}=0.140625, \qquad \frac{0.04}{1.8^2}=\frac{0.04}{3.24}. $$

The projector perturbation has both the off-diagonal entry and its transpose for each component. Therefore its expected squared Frobenius norm is twice the sum of those component variances: $$ 2\left(\frac{0.04}{3.24}+\frac{0.09}{0.64}\right)\simeq 0.30594. $$ The smaller gap, $0.8$, amplifies the corresponding fluctuation more strongly. The calculation does not require the two entries of $G$ to be independent; the expected squared Frobenius norm is the sum of the expected squares of these entries.

Worked Example: A Two-Dimensional Target Subspace

Let $M$ have eigenvalues $3.0,\ 2.0,\ 0.5$, and estimate the span of the first two eigenvectors. The gap at the boundary is $2.0-0.5=1.5$. For this target, the first-order projector fluctuation involves only the two cross-boundary entries $u_3^TGu_1$ and $u_3^TGu_2$, with denominators $3.0-0.5=2.5$ and $2.0-0.5=1.5$. Thus $$ \sqrt n(\hat P-P)\xrightarrow{d} \frac{u_3^TGu_1}{2.5}(u_3u_1^T+u_1u_3^T) +\frac{u_3^TGu_2}{1.5}(u_3u_2^T+u_2u_3^T). $$

There is no denominator involving $3.0-2.0$: rotations between the first and second eigenvectors do not alter their two-dimensional span. If the second eigenvalue were tied with the third, however, the boundary gap would vanish and this separated-subspace approximation would not apply.

Conditions, Interpretation, and a Common Pitfall

The assumptions specify which uncertainty the formulas describe. Fixed slice boundaries and $p_h>0$ allow each conditional covariance to be estimated from an increasing number of observations. The finite fourth moment ensures the within-slice covariance contributions have finite variance and a central limit theorem applies. If a slice probability is very small, the formulas remain asymptotic when it is positive, but the factor $1/p_h$ warns that finite-sample uncertainty may be large. Data-dependent boundaries, estimated standardization, or a number of slices that changes with $n$ require additional terms or different conditions.

Likewise, a separated eigenspace is essential for the stated projector limit. Individual eigenvectors within a group of equal eigenvalues are not uniquely defined, so their individual limits cannot be reported using the simple-eigenvalue formula. A whole cluster can still define a meaningful eigenspace when it is separated from the rest, which is why the projector result requires a boundary gap rather than distinct eigenvalues within the target group.

Common mistake: treating the slice proportions as fixed. This drops the $(I_h-p_h)A_h$ contribution and generally understates uncertainty. The correction is to include the probability influence in $\Psi$, even when the conditional covariance estimates are also being analyzed.
Common mistake: treating a small eigenvalue gap as harmless. The denominators $\lambda_j-\lambda_k$ show directly why this fails. A small gap magnifies operator noise; if the target boundary gap is zero, the stated eigenspace limit does not apply.
Common mistake: reporting a unique direction inside a tied eigenspace. The eigenvector may rotate within that space without changing the subspace. Report and analyze the subspace projector when the scientific target is the span.

The practical implication is that uncertainty has two linked sources: the sampling behavior of the sliced covariance summaries and the geometry of the resulting spectrum. The influence matrix describes the first; eigenvalue separation controls how it becomes uncertainty in directions or subspaces. In applications using estimated standardization or estimated slice boundaries, those estimation steps must be accounted for before treating these fixed-slice, population-standardized formulas as a complete variance calculation.

Check Your Understanding

Use the definitions and first-order formulas above to answer the following questions.

Why does the influence of a within-slice covariance estimate contain a factor $1/p_h$?
Which part of the SAVE operator influence function accounts for random slice proportions?
For an individual eigenvector, what condition is needed for the displayed first-order limit to have nonzero denominators?
Why does the projector limit for a target space omit terms describing rotations within that space?
What additional sources of variation are outside the fixed-slice, population-standardized setup used here?

Asymptotic Uncertainty for SAVE Eigenspaces

What You'll Learn