Asymptotic Distribution of Sample SIR Directions

Why SIR Directions Have Sampling Uncertainty

Sample SIR directions are computed from slice means and a predictor covariance matrix, all of which vary from sample to sample. Consequently, an estimated direction is not just a fixed vector with a small numerical error: it is a random estimate whose variability depends on the slice structure, the covariance estimate, and the separation between eigenvalues. The separation matters because two nearly equal eigenvalues can have directions that rotate substantially under small perturbations.

This tutorial studies the usual fixed-slice asymptotics. The response is divided into a fixed number of slices using fixed boundaries; the predictor observations are independent and identically distributed. If slice boundaries are estimated from the data, for example by sample quantiles, their randomness can contribute additional terms and is not included in the formulas below.

Population and Sample Generalized Eigenproblems

Let $X\in\mathbb{R}^p$ be a predictor with mean $\mu=E[X]$ and positive-definite covariance matrix $\Sigma=\operatorname{Var}(X)$. Let $G\in\{1,\ldots,H\}$ identify the response slice, with $q_h=P(G=h)>0$ and conditional mean $\mu_h=E[X\mid G=h]$. Write $\eta_h=\mu_h-\mu$. The population SIR matrix and its generalized eigenproblem are

B=\sum_{h=1}^H q_h\eta_h\eta_h^T, \qquad Bv_j=\lambda_j\Sigma v_j.

Order the eigenvalues as $\lambda_1\geq\cdots\geq\lambda_p\geq0$. Choose generalized eigenvectors so that $v_j^T\Sigma v_k$ equals 1 when $j=k$ and 0 otherwise. These vectors are the original-predictor-coordinate directions; equivalently, they are obtained by transforming eigenvectors of the whitened SIR matrix back to predictor coordinates.

Definition: For observations $X_1,\ldots,X_n$, let $\hat q_h$ be the proportion in slice $h$, $\bar X$ the overall sample mean, and $\hat\mu_h$ the sample mean in slice $h$. The fixed-slice sample SIR matrix and sample covariance are $$ \hat B=\sum_{h=1}^H\hat q_h(\hat\mu_h-\bar X)(\hat\mu_h-\bar X)^T, \qquad \hat\Sigma=\frac1n\sum_{i=1}^n(X_i-\bar X)(X_i-\bar X)^T. $$ The sample SIR eigenvalues and directions solve $\hat B\hat v_j=\hat\lambda_j\hat\Sigma\hat v_j$, with $\hat v_j^T\hat\Sigma\hat v_k=1$ for $j=k$ and 0 otherwise.

The divisor $n$ in $\hat\Sigma$ is convenient here; using $n-1$ gives the same first-order asymptotics. Assume $E\|X\|^4<\infty$, so the matrix influence functions below have finite second moments. Assume also that all slice probabilities are positive and that $\Sigma$ is positive definite.

Worked Example: A Two-Slice Population Problem

Suppose $p=2$, $\Sigma=I_2$, and two equally likely slices have conditional means $\mu_1=(1,0)^T$ and $\mu_2=(-1,0)^T$. The overall mean is zero, so $\eta_1=(1,0)^T$ and $\eta_2=(-1,0)^T$. Therefore

B=\frac12 \begin{pmatrix}1&0\\0&0\end{pmatrix} \begin{pmatrix}1&0\end{pmatrix} +\frac12 \begin{pmatrix}-1&0\\0&0\end{pmatrix} \begin{pmatrix}-1&0\end{pmatrix} = \begin{pmatrix}1&0\\0&0\end{pmatrix}.

Since $\Sigma=I_2$, the generalized eigenvalues are $1$ and $0$, with directions $v_1=(1,0)^T$ and $v_2=(0,1)^T$. The leading direction is identifiable because its eigenvalue is separated from the other one. The second direction is also distinct in this particular two-dimensional example, but in a higher-dimensional SIR problem the zero eigenvalue can have multiplicity greater than one, in which case its individual directions are not identifiable.

First-Order Fluctuations of the SIR Matrix

For an observation $X$ belonging to slice $g$, define the matrix-valued influence functions

\psi_B(X)=\eta_g\eta_g^T-B+ \eta_g(X-\mu_g)^T+(X-\mu_g)\eta_g^T, \qquad \psi_\Sigma(X)=(X-\mu)(X-\mu)^T-\Sigma.

Both have expectation zero. The first term of $\psi_B$ accounts for variation in slice proportions; the two cross-products account for estimation of the slice means. The influence function $\psi_\Sigma$ is the usual one for a covariance matrix when the mean is estimated.

Theorem: With fixed slices, positive slice probabilities, positive-definite $\Sigma$, and $E\|X\|^4<\infty$, the sample matrices satisfy $$ \sqrt n(\hat B-B)=\frac1{\sqrt n}\sum_{i=1}^n\psi_B(X_i)+o_p(1), \qquad \sqrt n(\hat\Sigma-\Sigma)=\frac1{\sqrt n}\sum_{i=1}^n\psi_\Sigma(X_i)+o_p(1). $$ Consequently, the vector of their entries has a joint multivariate normal limit, with covariance determined by the joint influence functions.

To verify the first expansion, regard each slice probability and conditional mean as functionals of the distribution. An infinitesimal empirical perturbation contributed by an observation in slice $g$ changes $q_h$ by $1_{\{g=h\}}-q_h$, and changes $\mu_h$ by $1_{\{g=h\}}(X-\mu_h)/q_h$. Since $B=\sum_hq_h\eta_h\eta_h^T$ and $\sum_hq_h\eta_h=0$, differentiating the centered expression gives

\sum_h(1_{\{g=h\}}-q_h)\eta_h\eta_h^T +\sum_hq_h\left\{\dot\mu_h\eta_h^T+\eta_h\dot\mu_h^T\right\} =\eta_g\eta_g^T-B+\eta_g(X-\mu_g)^T+(X-\mu_g)\eta_g^T.

The term from estimating the overall mean vanishes to first order because it is multiplied by $\sum_hq_h\eta_h=0$. The sample quantities are smooth functions of finitely many empirical slice probabilities, slice first moments, and overall first and second moments. Their first-order Taylor expansions therefore yield the stated representations. The multivariate central limit theorem applies to these empirical moments under the fourth-moment assumption, giving the joint normal limit. The covariance expansion follows by differentiating $E[(X-\mu)(X-\mu)^T]$; its derivative is $\psi_\Sigma(X)$. $\square$

Worked Example: An Influence-Function Contribution

Let $X$ be scalar and let there be two equally likely slices. In slice 1, $X$ takes values 0 and 2 with equal probability; in slice 2, it takes values $-2$ and 0 with equal probability. The slice means are 1 and $-1$, the overall mean is 0, and

B=\frac12(1)^2+\frac12(-1)^2=1, \qquad \Sigma=E[X^2]=\frac14(0+4+4+0)=2.

For an observation $X=2$ in slice 1, $\eta_1=1$ and $X-\mu_1=1$. Thus

\psi_B(2)=1^2-1+2(1)(1)=2, \qquad \psi_\Sigma(2)=(2-0)^2-2=2.

The only positive generalized eigenvalue is $\lambda=B/\Sigma=1/2$. Its first-order influence contribution is $(\psi_B-\lambda\psi_\Sigma)/\Sigma=(2-\tfrac12(2))/2=1/2$. This calculation already shows why treating the predictor covariance as fixed would give the wrong uncertainty: the covariance contribution is part of the eigenvalue fluctuation.

Asymptotic Distribution of Separated Eigenvalues and Directions

Consider an eigenvalue $\lambda_j$ that is simple, meaning it differs from every other generalized eigenvalue. The corresponding direction is identified only up to sign; for asymptotic statements, choose the sign of $\hat v_j$ to agree locally with $v_j$. Let $\Delta_B$ and $\Delta_\Sigma$ denote small perturbations of $B$ and $\Sigma$. Differentiating the generalized eigen-equation and the normalization condition gives the following result.

Theorem: If $\lambda_j$ is simple and $v_j^T\Sigma v_j=1$, then the first-order changes in its eigenvalue and normalized direction are $$ \dot\lambda_j=v_j^T(\Delta_B-\lambda_j\Delta_\Sigma)v_j, $$ $$ \dot v_j= \sum_{k\ne j}v_k \frac{v_k^T(\Delta_B-\lambda_j\Delta_\Sigma)v_j}{\lambda_j-\lambda_k} -\frac12v_j(v_j^T\Delta_\Sigma v_j). $$ In particular, the eigenvalue and direction have root-$n$ asymptotic distributions when $(\hat B,\hat\Sigma)$ have the joint root-$n$ limit above and $\lambda_j$ is separated from the other eigenvalues.

For the eigenvalue formula, differentiate $Bv_j=\lambda_j\Sigma v_j$, then left-multiply by $v_j^T$. Symmetry and $v_j^T\Sigma v_j=1$ cancel the terms involving $\dot v_j$, leaving $\dot\lambda_j=v_j^T(\Delta_B-\lambda_j\Delta_\Sigma)v_j$.

For the direction, write its change in the generalized eigenbasis as $\dot v_j=\sum_k c_kv_k$. Differentiating the eigen-equation and left-multiplying by $v_k^T$, for $k\ne j$, gives

(\lambda_k-\lambda_j)v_k^T\Sigma\dot v_j+ v_k^T(\Delta_B-\lambda_j\Delta_\Sigma)v_j=0, \qquad c_k=v_k^T\Sigma\dot v_j =\frac{v_k^T(\Delta_B-\lambda_j\Delta_\Sigma)v_j}{\lambda_j-\lambda_k}.

The denominator is $\lambda_j-\lambda_k$: the minus sign in the first equation is essential. Finally, differentiating $v_j^T\Sigma v_j=1$ gives $2c_j+v_j^T\Delta_\Sigma v_j=0$, so $c_j=-\tfrac12v_j^T\Delta_\Sigma v_j$. Combining the coefficients proves the direction formula. $\square$

Substituting $\Delta_B=\psi_B(X)$ and $\Delta_\Sigma=\psi_\Sigma(X)$ into these derivatives gives the influence functions for the eigenvalue and direction. For example,

\psi_{\lambda_j}(X)=v_j^T\{\psi_B(X)-\lambda_j\psi_\Sigma(X)\}v_j.

The central limit theorem then gives $\sqrt n(\hat\lambda_j-\lambda_j)\xrightarrow{d}N(0,\operatorname{Var}(\psi_{\lambda_j}(X)))$. The direction has a multivariate normal limit with influence function obtained from the displayed formula for $\dot v_j$. A small eigenvalue gap makes the direction’s variance large, even when the matrix estimates themselves are precise.

Worked Example: Checking the Directional Sign

Let $\Sigma=I_2$, $B=\operatorname{diag}(1/4,0)$, and consider the leading direction $v_1=(1,0)^T$, with $\lambda_1=1/4$. Let the perturbation of $B$ have off-diagonal entries $-0.02$, with $\Delta_\Sigma=0$. The other eigenpair is $v_2=(0,1)^T$, $\lambda_2=0$. The first-order component of the direction along $v_2$ is

\frac{v_2^T\Delta_Bv_1}{\lambda_1-\lambda_2} =\frac{-0.02}{0.25}=-0.08.

The sign can be checked directly. The perturbed matrix is $\begin{pmatrix}0.25&-0.02\varepsilon\\-0.02\varepsilon&0\end{pmatrix}$. If its leading eigenvector is proportional to $(1,a\varepsilon)^T$, the second row of the eigen-equation gives $-0.02\varepsilon\approx(0.25)(a\varepsilon)$, hence $a=-0.08$. The negative component agrees with the formula.

Clusters, Covariance Estimation, and Interpretation

An individual direction is not well determined when its eigenvalue is repeated. In that case, small perturbations can rotate the eigenvectors inside the repeated-eigenvalue space, so a limit for each chosen vector is generally not meaningful. A separated cluster of eigenvalues still defines an identifiable eigenspace: its eigenvalues must be separated from those outside the cluster, though they may coincide with one another. The first-order movement of that space is described by mixing between the cluster and its complement. For an eigenvalue $\lambda$ in a repeated cluster and an outside eigenvector $v_l$, the mixing coefficient for a cluster vector $v_j$ is $v_l^T(\Delta_B-\lambda\Delta_\Sigma)v_j/(\lambda-\lambda_l)$. Rotations within the cluster change the basis, not the subspace.

Estimated predictor covariance enters in two distinct ways. It changes the eigenvalue through the term $-\lambda_j\Delta_\Sigma$, and it changes the normalization and direction through both $-\lambda_j\Delta_\Sigma$ and the final normalization term. Thus, replacing $\hat\Sigma$ by $\Sigma$ in an uncertainty calculation is not justified merely because the covariance estimate is consistent.

Worked Example: How Covariance Estimation Changes a Direction

Take $B=\operatorname{diag}(0.4,0.1)$, $\Sigma=I_2$, and the leading pair $\lambda_1=0.4$, $v_1=(1,0)^T$. The second eigenvalue is $0.1$, so the gap is $0.3$. Suppose $\Delta_B$ has off-diagonal entries $0.03$, and $\Delta_\Sigma$ has off-diagonal entries $0.02$, with zero diagonal entries. The numerator for movement toward $v_2=(0,1)^T$ is

v_2^T(\Delta_B-\lambda_1\Delta_\Sigma)v_1 =0.03-(0.4)(0.02)=0.022.

The direction component is $0.022/(0.4-0.1)=0.073\overline{3}$. If covariance estimation were incorrectly ignored, the result would be $0.03/0.3=0.1$. Here the covariance contribution reduces the first-order rotation. This is a perturbation calculation: it describes the coefficient multiplying the perturbation size, not a claim that a finite perturbation has exactly that effect.

Common Pitfall

Wrong reasoning: “The numerator in the direction formula is $v_k^T(\Delta_B-\lambda_j\Delta_\Sigma)v_j$, so the component is that numerator divided by $\lambda_k-\lambda_j$.” This skips the sign obtained when the differentiated eigen-equation is rearranged. The correct equation is $(\lambda_j-\lambda_k)v_k^T\Sigma\dot v_j=v_k^T(\Delta_B-\lambda_j\Delta_\Sigma)v_j$, so the denominator is $\lambda_j-\lambda_k$. The two-dimensional calculation above verifies the sign directly.

A second practical warning is that a tiny eigenvalue gap signals unstable directions, not necessarily unstable estimated subspaces. When several eigenvalues are close, report or analyze their combined span rather than interpreting a particular basis vector as uniquely determined. Likewise, formulas here assume fixed slice boundaries; random boundaries require incorporating their estimation into the matrix influence functions before applying eigenvalue perturbation theory.

Key takeaway: Under fixed-slice and moment regularity conditions, sample SIR matrices fluctuate at rate $n^{-1/2}$. A simple, separated generalized eigenvalue and its direction inherit that rate through explicit derivatives involving both $\hat B$ and $\hat\Sigma$. Eigenvalue gaps govern directional stability, while a separated cluster—not an arbitrary basis within it—is the appropriate target when eigenvalues coincide.

Check Your Understanding

Use the definitions and perturbation formulas to check each answer, paying particular attention to eigenvalue gaps and covariance terms.

Why does the influence function for $B$ contain both a slice-proportion term and terms involving within-slice residuals?
For a simple eigenvalue, write its first-order change when both $B$ and $\Sigma$ are perturbed.
In the direction formula, why is the denominator $\lambda_j-\lambda_k$, rather than $\lambda_k-\lambda_j$?
What happens to an individual direction when its eigenvalue is repeated, and what target can remain identifiable?
Why can treating the predictor covariance as fixed give an incorrect standard error for a sample SIR eigenvalue or direction?