arXiv CSopen access

Feature Bagging Provides Stability

Feature Bagging Provides Stability Yuheng Ma†1,3 , Qiang Sun‡2,3 1 School of Statistics, East China Normal University, 2 University of Toronto, 3 MBZUAI

We study feature bagging through the lens of algorithmic stability. Feature bagging is an ensemble strategy that aggregates base learners trained on randomly subsampled feature subsets, possibly in a data-dependent manner. We introduce feature instability (FI), the feature-axis analogue of instance instability (II), which measures sensitivity to removing a single feature. Smaller values of II or FI correspond to stronger stability, and our experiments show that FI captures generalization-relevant information complementary to II. Within this framework, we analyze feature bagging in both a parametric linear model and a model-free setting inspired by recursive feature subsampling in random forests. In both settings, we establish formal guarantees showing that feature bagging improves the relevant stability relative to its non-bagged counterpart, with larger improvements under more aggressive subsampling. We further show that a modest number of bagging rounds is sufficient to approach the infinite-bagging stability level. Keywords: stability, bagging, feature instability, random forests, random forward selection. Date: July 29, 2026

Keywords: stability, bagging, feature instability, random forests, random forward selection.

Contents 1

Introduction 1.1 Related Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 3

2

Basic Concepts 2.1 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Why Do We Care about Feature Stability? . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Instance and Feature Bagging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 4 6 6

3

Stability under Model Assumptions 3.1 Basic Settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7 7 8

4

Model-free Feature Stability 4.1 Model-free Stability in General Spaces . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Recursive Subsampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Random Forward Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10 11 12 14

† ‡

E-mail: [email protected] E-mail: [email protected]

1

4.3

4.2.2 Random Forests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Finite Bagging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

Conclusion

1

Introduction

16 19 20

Stability is a fundamental prerequisite for trustworthy machine learning (Murdoch et al., 2019; Yu and Kumbier, 2020; Xing et al., 2021; Zhou et al., 2023), and is closely connected to generalization (Bousquet and Elisseeff, 2002), uncertainty quantification (Wang et al., 2023), and model selection (Meinshausen and Bühlmann, 2010; Nogueira et al., 2018). For a predictor f trained on n samples in dimension d, an informal way to capture instance-wise instability is to measure how much its prediction changes when one training instance is removed: n

2 1 X f (x) − f −i,: (x) n i=1

(1.1)

where f −i,: denotes the predictor obtained after removing the i-th training instance. Smaller values indicate that the algorithm is more stable to single-instance perturbations. We use this quantity only as motivation; Definition 2.1 formalizes ϕ-instance stability and names the associated instance-instability measure. Regularization (Bousquet and Elisseeff, 2002) and bagging (Soloff et al., 2024a) are two standard approaches for improving instance stability. Feature bagging has been recognized as a potential driver of ensemble success (Breiman, 2001; Sutton et al., 2006; LeJeune et al., 2020), but the mechanism behind its empirical effectiveness remains debated. One promising explanation is that feature bagging acts as a form of regularization (Mentch and Zhou, 2020; Curth et al., 2024). This viewpoint makes stability a natural lens for studying feature bagging. Instance stability is a natural starting point because it measures sensitivity to changes in the training instances. However, it does not directly capture sensitivity to changes in the feature set, precisely the axis on which feature bagging operates; see Section 2.2 for empirical evidence that feature-axis perturbations carry generalization-relevant information that is complementary to instance stability. These considerations motivate our central question: Does feature bagging improve stability, and if so, how? Answering this question requires a feature-axis analogue of classical instance instability. We therefore introduce feature instability (FI), the leave-one-feature-out analogue of (1.1): d

2 1 X f (x) − f :,− j (x) , d j=1

(1.2)

where f :,− j denotes the predictor obtained after removing the j-th feature. This informal quantity is the feature-axis analogue of instance instability; the formal definitions of feature-instability and ϕ-stability are given later in Definition 2.2. We answer this question through the following contributions: 1. We introduce FI as the feature-axis analogue of classical II. This gives a two-axis view of stability tailored to feature bagging: II measures sensitivity to removing training instances, whereas FI measures sensitivity to removing features. We further show empirically that FI carries

2

generalization-relevant information that is complementary to II, motivating FI as a distinct object of study. 2. In linear regression, we give an exact asymptotic characterization of how bagging affects II and FI. The analysis explains how ensemble averaging improves stability and how the stability gains depend on the sampling of both instances and features. It also reveals that feature-level perturbations have their own scaling behavior, governed by the ambient dimension, which is invisible from II alone. 3. Beyond parametric models, we develop a model-free theory for FI under feature bagging, with particular emphasis on recursive feature subsampling. This theory applies to algorithms such as random forward selection and random forests, where features are resampled repeatedly and data-dependently. It shows that more aggressive feature subsampling yields stronger featurestability guarantees, identifies regimes in which feature bagging is provably more stable than its non-bagged counterpart, and establishes that a number of bagging rounds proportional to the dimension is enough to approach the infinite-bagging stability level. The rest of the paper proceeds as follows. Section 2 formalizes II, FI, and bagging. Section 3 analyzes the linear-model setting, while Section 4 develops the model-free theory for recursive feature subsampling. Proofs, derivations, and additional experiments are deferred to the appendix.

1.1

Related Works

Stability. Algorithmic stability formalizes the idea that a learned predictor should not change too much after a small change to the training data. Classical formulations usually modify the instance axis, either by removing one sample or by replacing it with an independent copy (Bousquet and Elisseeff, 2002; Elisseeff et al., 2005; Kutin and Niyogi, 2002; Mukherjee et al., 2006; Liu et al., 2017). These notions are central to stability-based generalization analysis (Bousquet and Elisseeff, 2002; Feldman and Vondrak, 2018, 2019; Bousquet et al., 2020). A closely related line studies averaging stability, which averages the sensitivity over all single-sample changes (Shalev-Shwartz et al., 2010; Lei and Ying, 2020; Soloff et al., 2024a). Recent work has used this viewpoint to analyze why bagging can stabilize learning algorithms (Adrian et al., 2024; Soloff et al., 2024b,a,c; Liang et al., 2025). Our definitions follow the removal-based convention. For instance stability, we compare the predictor trained on D with the predictor trained on D−i,: . Replacement-based variants can be related to the removal-based definition by a triangle-inequality argument (Bousquet and Elisseeff, 2002; Soloff et al., 2024a). The same issue appears on the feature axis. We use feature removal as the primary convention because feature replacement requires specifying how the replacement coordinate is generated, for −j example unconditionally or conditionally on (Yi , xi ). Bagging. Instance bagging is widely used across tree-based models (Breiman, 1996; Geurts et al., 2006), nearest-neighbor methods (Biau et al., 2010; Cai et al., 2025), and neural-network ensembles (Hansen and Salamon, 1990; Perrone and Cooper, 1995; Lakshminarayanan et al., 2017; Zaidi et al., 2021). Feature bagging is most prominent in tree-based methods, especially random forests, where each split is chosen after subsampling the candidate features (Breiman, 1996, 2001; Geurts et al., 2006). Recent high-dimensional linear-model analyses provide a more controlled view of feature subsampling and ensembling. Sketched ridgeless regression studies how downsampling changes interpolation behavior (Chen et al., 2023); related reference-panel estimators analyze the exact risk of regularization induced by an auxiliary feature panel (Su et al., 2024); and bagged linear interpolators show how averaging sketched 3

𝑌

𝑋

𝑋

𝑌

𝑋

𝑌

Instance subsampling

Instance subsampling

Tree B

Tree 1 Feature subsampling

Feature subsampling

Root

Root

Node

Split

Feature subsampling

Node

(a) Instance subsam- (b) Feature subsam- (c) Instance and feature subsampling. pling. pling.

Split

Split

Feature subsampling

Node

Node

Split

Feature subsampling

Feature subsampling

Split

Split

(d) Bagging scheme in random forests.

Figure 1: Illustration of subsampling and bagging schemes. In the random forest example (d), all splits within a single tree share the same instance subsample (yellow lines), while each split uses an independently drawn feature subsample (green columns). estimators controls interpolation-driven variance (Wu and Sun, 2025). These works connect feature subsampling to prediction risk. We build on this literature by asking the corresponding stability question: whether feature bagging also stabilizes predictions along the feature axis.

2

Basic Concepts

2.1

Stability

We begin by fixing notation for the two perturbation axes considered in this paper. Let D = (X, y) be a dataset with n instances and d features. The data matrix X ∈ Rn×d has rows xi ∈ X = Rd , and the labels are collected in y ∈ Yn . For integer K, let [K] = {1, . . . , K}. For i ∈ [n], denote by D(i),: and D−i,: the datasets obtained by replacing or removing the i-th sample from D, with corresponding data matrices X (i),: ∈ Rn×d and X −i,: ∈ R(n−1)×d . Similarly, for j ∈ [d], let X :,( j) and X :,− j denote the data matrices with the j-th feature replaced or removed, and write D:,( j) and D:,− j for the associated datasets. Throughout the paper, our primary convention is feature or instance removal. Replacementbased variants are also common in the stability literature, and we discuss their relationship to our removal-based convention in Section 1.1. Following the standard algorithmic-stability setup, we let a learning algorithm A map a dataset and a random seed to a predictor f = A(D, ξ) : X → Y, where ξ collects all sources of algorithmic randomness. Conditional on ξ, the algorithm is deterministic; for example, ξ may encode random initialization, data shuffling, or subsampling. Classical instance stability quantifies the sensitivity of A to instance-wise perturbations of D. We recall this notion first, and defer a discussion of related formulations to Section 1.1.

Definition 2.1 An algorithm A is ϕ-instance stable if n

2  1 X  II := E f (x) − f −i,: (x) ≤ ϕ2 , n i=1

(2.1)

where f = A(D; ξ) and f −i,: = A(D−i,: ; ξ) are evaluated using the same algorithmic randomness ξ.

The expectation is interpreted according to the stability convention being used. In the uniform convention, the expectation is over algorithmic randomness and the bound holds for fixed D and x. In distributional variants, the expectation may also average over the randomness in D and x. With some abuse of notation, we call the left-hand side of (2.1) the instance instability of A, with smaller 4

values indicating greater instance stability. The uniform convention is the standard instance-stability convention in the algorithmic-stability literature, while the distributional convention is weaker and is used in some of our average-case results, particularly in Section 3. The average over i makes the criterion invariant to permutations of the dataset, even when A itself is not permutation-invariant. If A is permutation-invariant, then (2.1) reduces to E[( f (x) − f −i,: (x))2 ] ≤ ϕ2 for any i ∈ [n], under the same expectation convention. To study feature bagging, we need the analogous notion along the feature axis. Rather than modifying training instances, we remove one feature from the input representation and measure the resulting change in prediction. Definition 2.2 An algorithm A is ϕ-feature stable if FI :=

d 2  1 X  E f (x) − f :,− j (x) ≤ ϕ2 , d j=1

(2.2)

  where f = A(D; ξ) and f :,− j = A D:,− j ; ξ are evaluated using the same algorithmic randomness. By convention, we regard f :,− j as a predictor on the original feature space that ignores the j-th feature.

The expectation follows the same convention as in Definition 2.1. We call the left-hand side of (2.2) the feature instability of A, with smaller values indicating greater feature stability. The preceding definition fixes a shared seed across the full-data and feature-removed runs. We next allow these two runs to use coupled randomness. This is useful when removing a feature changes the admissible resampling space, as in the linear-model analysis of Section 3. Definition 2.3 (Coupled feature stability) For a dataset D, let Ξ(D) denote the seed space used to run A on D. For each feature j, let Γ j (D) be a coupling on Ξ(D) × Ξ(D:,− j ), that is, a joint law for the seeds (ξ, ξ j ) used in the full-data and feature-removed runs. An algorithm A is ϕ-feature stable under the coupling family Γ = {Γ j (D)}dj=1 if d 2  1 X  E A(D; ξ)(x) − A(D:,− j ; ξ j )(x) ≤ ϕ2 , d j=1

(2.3)

where, in the j-th summand, (ξ, ξ j ) ∼ Γ j (D), and any additional averaging follows the same convention as above. The left-hand side is the coupled feature instability under Γ.

Definition 2.2 is recovered as the shared-randomness special case of Definition 2.3: after identifying the two seed spaces, Γ j (D) is the diagonal coupling ξ j = ξ. Unless otherwise stated, FI refers to this shared-randomness convention. Section 3 is an exception. In the linear-model analysis, the reducedfeature run shares the unaffected sketching randomness but redraws the feature sketch from the reduced feature set, corresponding to a affected-axis resampling coupling with ξ j , ξ. In contrast, Section 4 returns to the shared-randomness coupling ξ j = ξ, because its goal is a conditional algorithmic stability guarantee. The FI criterion in (2.2) differs from the II criterion in (2.1) in two important respects. First, FI perturbs only the input matrix and leaves the response vector y unchanged. Second, the average over 5

features is essential: features are generally non-exchangeable, and learning algorithms are typically not permutation-invariant along the feature dimension.

2.2

Why Do We Care about Feature Stability?

Feature stability is motivated by several practical and methodological concerns. First, relevant features may be absent, unavailable, or deliberately excluded; feature inclusion can therefore depend on preprocessing decisions, acquisition constraints, and analyst judgment (Jeng et al., 2024; Ma et al., 2025; Shen and Xiu, 2025; Poudel et al., 2025). The FI measure quantifies how sensitive an algorithm is to this feature-level variation. Second, feature entries may themselves be perturbed because of corruption (McWilliams et al., 2014), missingness (Little and Rubin, 2019; Chen and Xu, 2023), or measurement error (Hou et al., 2026). Algorithms with stronger feature stability are therefore more robust to such perturbations. Beyond these practical considerations, we empirically examine whether FI carries information about generalization, paralleling a central motivation for the literature on instance stability (Bousquet and Elisseeff, 2002; Hardt et al., 2016). In a controlled synthetic regression benchmark, specifically the MARSadd random-forest benchmark (Friedman, 1991; Mentch and Zhou, 2020; Curth et al., 2024), where the signal combines a weak component distributed across many features with a concentrated nonlinear component, we vary model complexity and subsampling behavior in random forest models, and measure how FI and II account for variation in the generalization gap. For each configuration of max features, max samples, and max depth, we compute the train-test MSE gap together with II and FI by retraining after removing one instance or one feature. We then quantify whether FI explains generalization variation beyond II using linear R2 , mutual information, and residual-based partial mutual information, with Gaussian and permuted-noise controls. Section S.1 provides the detailed experiments and analysis. The results indicate that FI captures generalization-relevant variation that is not explained by II alone. Across Tables S.2–S.5, the conditional contribution of FI after accounting for II remains positive in all displayed synthetic and real-data settings. The noise-control results in Tables S.6–S.9 further show that replacing FI with unrelated Gaussian or permuted noise does not replicate this signal. These findings indicate that FI and II provide complementary measures for capturing generalization-relevant variation.

2.3

Instance and Feature Bagging

We now formalize the bagging schemes analyzed in the paper. For both instances and features, we use sub-bagging: each base learner is trained on a uniformly sampled fixed-size subset without replacement. Other bagging schemes can lead to similar theoretical conclusions (Soloff et al., 2024a), but fixed-size subsampling is especially natural for features. Classical bagging (Breiman, 1996) and Poissonized bagging (Oza and Russell, 2001) sample with replacement, which can select the same feature multiple times and create non-identifiability issues. Bernoulli sub-bagging (Harrington, 2003; Wu and Sun, 2025) produces random subset sizes, which can complicate models that require a fixed input dimension, such as neural networks. Instance bagging. We first consider instance bagging, the classical bagging setting. In each sub-bagging round, m of the n instances are sampled uniformly without replacement. Let p = m/n denote the instance(b) ⊤ subsampling ratio, and write µ(b) = (µ(b) 1 , . . . , µm ) for the indices selected in the b-th subsample. For (b) any algorithm A0 , we define its evaluation on the b-th subsample as A(D, ξ(b) ) := A0 (Dµ ,: ), where (b) Dµ ,: = {(xµ(b) , yµ(b) ), . . . , (xµ(b) , yµ(b) )}. See Figure 1(a) for an illustration. This construction defines a 1

1

m

m

randomized algorithm A, with randomness induced by instance subsampling. Let f (b) := A(D, ξ(b) ) be

6

the predictor trained in the b-th round. The bagged predictor averages B such predictors: B

f B (x) :=

1 X (b) f (x). B b=1

(2.4)

We refer to B as the number of bagging rounds. The infinite-bagging limit is the corresponding expectation over the subsampling randomness: A∞ (D)(·) := Eξ(b) [ f (b) (·)] = Eξ(b) [A(D, ξ(b) )(·)].

(2.5)

We write its output as f ∞ , so f ∞ (x) = A∞ (D)(x) = Eξ(b) [ f (b) (x)]. This predictor is the large-B limit of f B. Feature bagging. Feature bagging is defined analogously. Let s ≤ d denote the number of selected features and q = s/d the subsampling ratio. In each round, a subset ν (b) = (ν1(b) , . . . , ν s(b) ) is drawn uniformly at random without replacement. The corresponding randomized algorithm is A(D, ξ(b) ) := (b) (b) A0 (D:,ν ). Here, D:,ν denotes the dataset restricted to the selected features and is illustrated in Figure 1(b). The finite and infinite bagged predictors are defined as in (2.4) and (2.5). At prediction time, the same test point is evaluated by all base learners, but the b-th learner only receives the coordinates (b) xν that match its sampled feature subset. The final prediction is the average of these base-learner predictions. This differs from instance bagging, where every base learner receives the full test point. General bagging. Instance and feature bagging can be combined by jointly subsampling both axes, yielding an ensemble whose randomness arises from both sources; see Figure 1(c). Section 3 explores their interaction. In random forests, feature bagging is implemented recursively rather than once per tree: as Figure 1(d) illustrates, all splits in a tree share the same instance subsample, while feature subsampling is performed independently at each node. Section 4 analyzes this recursive feature-subsampling strategy.

3

Stability under Model Assumptions

This section studies the effect of bagging on the II and FI measures in linear regression, where explicit model assumptions allow sharp comparisons.

3.1

Basic Settings

Generating Distribution. We consider the standard linear regression model (LeJeune et al., 2020; Chen et al., 2023): ∗ yi = x⊤ i β + εi ,

(3.1)

where xi ∈ Rd is the feature vector, yi ∈ R is the response, εi ∼ N(0, σ2 ) is Gaussian noise, and (xi , εi )ni=1 are i.i.d. copies of (x, ε). Additionally, we assume that x ∼ N(0, I), β ∗ ∼ N(0, I/d), and that x, β ∗ , and the noise ε are mutually independent. The Gaussian assumptions on x and ε can be generally relaxed to bounded moment conditions under which the same theoretical results hold, a phenomenon known as universality (Hastie et al., 2022; Chen et al., 2023; Wu and Sun, 2025). To avoid introducing technical overhead, we adopt the Gaussian setting. Let X = (x1 , . . . , xn )⊤ denote the feature matrix, y = (y1 , . . . , yn )⊤ the response vector, and ε = (ε1 , . . . , εn )⊤ the noise vector. We work under the proportional asymptotic regime, where d, n → ∞ with d/n → γ for some constant γ. This regime has been widely adopted in recent studies analyzing the exact risk behavior of linear models; see, e.g., Hastie et al. (2022); Chen et al. (2023); Wu and Sun (2025). 7

Table 1: Closed-form asymptotic limits for the four ∆ℓ -terms in Theorem 3.1.

Setting and Case

∆V= ℓ

n∆V, ℓ

∆ℓB=

n∆ℓB,

II (ℓ = −i), γq < p

2γq(1−p) (p−γq)(1−γq) 2γq(p−p2 ) (γq−p)(γq−p2 ) 2γpq(1−q) (p−γq)(p−γq2 ) 2γp(1−q) (γq−p)(γ−p)

γq2 (1−γq2 )2 γp2 (γ−p2 )2 q2 (1−γq2 )2 p2 (γ−p2 )2

2γq(1−q)(1−p) (p−γq)(1−γq) 2p(1−p)(γ2 q−2γpq+p2 ) γ(γq−p)(γq−p2 ) 2pq(1−q)(p+γ−2γq) (p−γq)(p−γq2 ) 2p(1−q) γq−p

γq2 (1−q)2 (1−γq2 )2 p2 (γ−p)2 γ(γ−p2 )2 q2 (1+γ−2γq) γ(1−γq2 )2 p2 (1−2p+γ) (γ−p2 )2

II (ℓ = −i), γq > p FI (ℓ = − j), γq < p FI (ℓ = − j), γq > p

Bagged least square estimator. Following Wu and Sun (2025), we reformulate the bagged least square estimator as an average of sketched least square estimators. As described in Section 2.3, let µ(b) and ν (b) denote the instance and feature indices selected in the b-th subsample. Define the sketching matrices Ub ∈ Rn×m and Vb ∈ Rd×s as the corresponding selection matrices: Ub has ones at entries (µ(b) i , i) (b) for i ∈ [m], and zeros elsewhere; similarly, Vb has ones at entries (ν j , j) for j ∈ [s]. The matrix Ub⊤ XVb ∈ Rm×s extracts a submatrix of X formed by the sampled rows and columns. For each subsample, we compute the sketched minimum-norm least squares estimator  † β (b) = Vb Ub⊤ XVb Ub⊤ y, where (·)† denotes the Moore–Penrose pseudoinverse. The final bagged estimator is obtained by averaging PB β (b) . over B such subsamples β := B1 b=1 Affected-axis resampling coupling. This section uses an affected-axis resampling coupling for the sketching randomness: after an instance or feature is removed, we redraw the sketching matrix on the perturbed axis from the reduced index set, while keeping the sketching matrix on the other axis fixed. Thus the reduced run is neither a shared-seed run nor a fully independent rerun. After removing the i-th instance from the dataset, we independently resample a sketching matrix Ub′ from the reduced index set [n] \ {i}. The resulting bagged estimator after instance removal is β

−i,:

B † 1 X  ′⊤ := Vb Ub XVb Ub′⊤ y. B b=1

(3.2)

Similarly, removing the j-th feature corresponds to resampling Vb′ from [d] \ { j}, yielding B

β

3.2

:,− j

† 1 X ′ ⊤ := Vb Ub XVb′ Ub⊤ y. B b=1

(3.3)

Stability Analysis

P Under this affected-axis resampling coupling, we study the expected instability measures 1n ni=1 Ex,β,β−i,: [(β ⊤ x− P (β −i,: )⊤ x)2 ] for II and d1 dj=1 Ex,β,β:,− j [(β ⊤ x−(β :,− j )⊤ x)2 ] for FI. Since the features are isotropic, taking expectation over x reduces both quantities to   ℓ 2 Instability = E β − β 2 , ℓ ∈ {−i, − j}, (3.4)

8

where the subscript of Eβ,βℓ is omitted for brevity, ℓ = −i denotes i-th instance removal, and ℓ = − j denotes j-th feature removal. We slightly abuse notation by writing β ℓ for either β −i,: or β :,− j . Rather than upper bounding this instability, we characterize the exact limiting value of (3.4) under proportional asymptotics. Our first result gives a finite-B decomposition and the corresponding asymptotic limits.

Theorem 3.1 Assume the Gaussian linear model in Section 3.1. Then, for ℓ ∈ {−i, − j},   V= h i   ∆ℓ ∆B= B − 1 B, B − 1 V, ℓ 2 2 + ∆ℓ  + ℓ + ∆ℓ . E ∥β − β ∥2 = σ  B B B B

(3.5)

Moreover, as n, d → ∞ with d/n → γ, 1 B= V, B, ∆V= ℓ , ∆ℓ → Θ(1), and ∆ℓ , ∆ℓ → Θ( ), n

(3.6)

V, B= where Θ(·) denotes asymptotic equality up to constants. The exact asymptotic limits of ∆V= ℓ , n∆ℓ , ∆ℓ , B, and n∆ℓ are given in Table 1 for instance removal and feature removal.

The proof of Theorem 3.1 is given in Appendix S.2. The theorem decomposes the expected instance and feature instability of the bagged least-squares estimator into noise-variance terms and signal-bias terms. We use the superscripts V and B for these two contributions, respectively; the superscript B denotes bias and should not be confused with the number of bagging rounds. The = terms compare paired estimators built from the same subsampling realization and receive weight 1/B, whereas the , terms compare different submodels and receive weight (B − 1)/B. This decomposition separates the effects of subsampling and bagging. When B = 1, the estimator is a single subsampled minimum-norm least-squares fit, so only the same-submodel terms remain. These terms inherit the interpolation behavior of an individual subsampled estimator and become large near the effective interpolation threshold γq = p, as shown in Figure 2. Both II and FI have a double-descentshaped instability curve, with a peak near γq = p, where the subsampled design matrix is nearly singular. Changing the subsampling ratios p and q moves this threshold through the effective aspect ratio γq/p, but does not average away the same-submodel instability. This is the main contrast with the bagged curves in Figure 3. Increasing B downweights these same-submodel terms and transfers mass to the cross-submodel B, terms ∆V, ℓ and ∆ℓ . In the large-B limit, the cross terms dominate and are of order 1/n, yielding the dashed curves in Figure 3. Thus subsampling changes the effective interpolation threshold and the limiting constants, while bagging averages away the instability carried by individual subsampled estimators. Classical stability-based generalization bounds (Bousquet and Elisseeff, 2002; Elisseeff et al., 2005; Mukherjee et al., 2006) do not directly apply to linear regression because the output is unbounded. Nonetheless, the generalization-error curves observed in (LeJeune et al., 2020; Wu and Sun, 2025) show a related qualitative pattern: subsampling shifts the interpolation threshold, while averaging over subsampled estimators suppresses the associated peak. The closed forms further show that feature removal has a different natural scale from instance removal. In the large-B limit, most cross terms are reported on the 1/n scale, but the feature-removal bias term also reflects the number of available coordinates. In the underparameterized regime γq < p, n∆−B,j →

q2 (1 + γ − 2γq) q2 ∼ γ γ(1 − γq2 )2 9

(γ ↓ 0).

p=0.4, q=0.4 p=0.6, q=0.6 p=0.8, q=0.8 p=1, q=1

8 6 4 2 00.0

0.5

1.0

1.5

2.0

2.5

10 Instance Instability

Instance Instability

10

8 6 4 2 00.0

3.0

p=0.4, q=0.4 p=0.6, q=0.6 p=0.8, q=0.8 p=1, q=1

0.5

1.0

1.5

2.0

0.5

1.0

1.5

2.0

2.5

3.0

2.5

3.0

(b) II with instance subsampling

2.5

Feature Instability ( n 1)

Feature Instability ( n 1)

(a) II with both subsampling

14 12 10 8 6 4 2 0 0.0

p=0.4, q=1 p=0.6, q=1 p=0.8, q=1 p=1, q=1

3.0

(c) FI with both subsampling

14 12 10 8 6 4 2 0 0.0

p=1, q=0.4 p=1, q=0.6 p=1, q=0.8 p=1, q=1

0.5

1.0

1.5

2.0

(d) FI with feature subsampling

Figure 2: Instability of the single subsampled minimum-norm OLS estimator as a function of the aspect ratio γ. Colors indicate subsampling ratios (p, q); solid lines denote theoretical predictions and crosses denote empirical averages over 300 realizations with σ = 1, n = 400, and d ∈ {100, 200, . . . , 900}. Vertical dashed lines indicate the interpolation threshold γq = p. Since d = γn, the unscaled term behaves as ∆−B,j ≍ q2 /d. Thus this component of FI is controlled by the feature dimension, not only by the sample size: removing one feature is more consequential when d is small relative to n, while the same feature-specific contribution vanishes as γ → ∞. Finally, Theorem 3.1 quantifies the finite-B averaging gain, as illustrated in Figure 3. For small B, B= the 1/B-weighted same-submodel terms ∆V= ℓ and ∆ℓ remain visible. As B increases, the instability approaches its order-1/n large-B limit; away from interpolation, taking B of order n is already sufficient to make further bagging rounds have diminishing returns. Near the interpolation boundary, however, the same-submodel constants can be large, so larger ensembles may still be needed to flatten the finite-B peak. The limiting formulas also show that p and q do not enter symmetrically: in the underparameterized regime γq < p, the leading large-B constants are primarily controlled by q, whereas in the overparameterized regime γq > p, they are primarily controlled by p.

4

Model-free Feature Stability

In this section, we investigate FI under model-free settings. We first develop a general-space featurestability tool, with the scalar R result appearing as a special case (Section 4.1), and then extend it to

10

6 4 2 00.0

0.5

1.0

1.5

2.0

2.5

20 15 10 5 0 0.0

3.0

(a) II with both subsampling p=0.4, q=1 p=0.6, q=1 p=0.8, q=1 p=1, q=1

Instance Instability ( n 1)

8 6 4 2 00.0

0.5

1.0

1.5

2.0

2.5

0.5

1.0

1.5

2.0

2.5

600

3.0

(b) II with instance subsampling

30

400 300 200

40

100

20

0 00.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

3.0

(e) II across bagging rounds

p=1, q=0.4 p=1, q=0.6 p=1, q=0.8 p=1, q=1

25 20 15 10 5 0 0.0

0.5

1.0

1.5

2.0

2.5

B=1 B=5 B=10 B=20 B=50 B=100 B=200 B=

500

(c) FI with both subsampling

Feature Instability ( n 1)

10

p=0.4, q=0.4 p=0.6, q=0.6 p=0.8, q=0.8 p=1, q=1

25

3.0

(d) FI with feature subsampling

600 Feature Instability ( n 1)

Instance Instability ( n 1)

8

30

Instance Instability ( n 1)

p=0.4, q=0.4 p=0.6, q=0.6 p=0.8, q=0.8 p=1, q=1

Feature Instability ( n 1)

10

B=1 B=5 B=10 B=20 B=50 B=100 B=200 B=

500 400 300 200

40

100

20

0 00.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

(f) FI across bagging rounds

Figure 3: Instability E[∥β − β ℓ ∥22 ] of the bagged least-squares estimator. Panels (a)–(d) use B = 200 bagging rounds and vary the aspect ratio γ; solid lines denote theoretical predictions, crosses mark empirical averages, and dashed curves show the infinite-bagging limit. Panels (e)–(f) compare different numbers of bagging rounds on an n−1 scale; colors indicate B, solid lines denote theoretical predictions, and crosses mark empirical averages. All panels use the same simulation setting as in Figure 2. recursive subsampling (Section 4.2). These results are applied to random forward selection and random forests in Sections 4.2.1 and 4.2.2, and extended to finite bagging rounds in Section 4.3. Unlike linear regression analysis in Section 3, the model-free results use the shared-randomness, conditional-on-thedataset coupling in Definitions 2.1–2.2, namely ξ j = ξ. The aim is therefore a uniform algorithmic guarantee rather than an exact distributional perturbation calculation.

4.1

Model-free Stability in General Spaces

This subsection develops the general-space stability tool used later for recursive feature subsampling, used in random forests (Breiman, 2001; Zhang et al., 2024). As illustrated in Figure 1(d), each tree uses a fixed subset of instances, while features are subsampled recursively at each node. This recursive, data-dependent feature subsampling makes feature-bagging analysis more challenging than instancebagging analysis. To isolate the feature-subsampling effect, we focus on FI and first present the scalar R guarantee before giving its Hilbert-space generalization. Proposition 4.1 (FI of feature bagging) Assume the base algorithm A0 produces predictors satisfying f (x) ∈ Y ⊆ [−M, M] for every x ∈ X, and let q < 1. Then the infinitely feature-bagged algorithm A∞ based on A0 is ϕ-feature stable whenever ϕ2 ≥

M2 q . d−11−q

Proposition 4.1 is the feature-wise dual of the infinite instance-bagging guarantee of Soloff et al. (2024a). Under the same bounded-output assumption, their result states that infinite instance bagging 11

2

p M with subsampling ratio p < 1 is ϕ-instance stable whenever ϕ2 ≥ n−1 1−p . Thus, their guarantee controls II under instance subsampling, whereas Proposition 4.1 controls FI under feature subsampling; formally, (n, p, II) is replaced by (d, q, FI). This feature-side guarantee is the scalar starting point for the general-space and recursive results below. We now lift the scalar feature-bagging result to general output spaces. Consider algorithms of the form w = A(D, ξ) whose outputs lie in a Hilbert space H with inner product ⟨·, ·⟩ and induced norm ∥ · ∥H . We assume that the output belongs to a convex, bounded set W ⊂ H and that the zero element of H lies in W. The associated feature instability over W is measured using ∥ · ∥H .

Definition 4.2 An algorithm A with output in W is ϕ-feature stable if, for all datasets D = {(xi , yi )}ni=1 , FIH :=

d i 1X h E ∥w − w− j ∥2H ≤ ϕ2 , d j=1

where w = A(D; ξ) and w− j = A(D:,− j ; ξ). (b)

The bagged estimator is defined analogously. Let w(b) = A(D, ξ(b) ) = A0 (D:,ν ) be the output of PB the b-th feature subsample. Since W is convex, w B = B−1 b=1 w(b) also lies in W. The infinite-bagged ∞ estimator is w := E[A(D, ξ)]. This framework extends the scalar case W ⊆ R above. We use the following Hilbert-space analogue of the bounded-range assumption of Soloff et al. (2024c).

Assumption 1 Let W be a closed and convex subset of a Hilbert space with norm ∥ · ∥H . We assume that its radius rad(W) := inf w∈W supw′ ∈W ∥w − w′ ∥H is finite. When H = R, ∥ · ∥H = | · |, and W = [−M, M], Assumption 1 recovers the bounded-range condition used in Proposition 4.1.

Proposition 4.3 (General-space FI of feature bagging) Assume Assumption 1 holds, and let q < 1. For any base algorithm A0 with output in W, infinite feature bagging A∞ is ϕ-feature stable whenever ϕ2 ≥

rad2 (W) q . d−1 1−q

Proposition 4.3 generalizes Proposition 4.1 by replacing the scalar range bound M with the Hilbertspace radius rad(W). This general-space form is the key input for the recursive procedures analyzed next, including random forward selection and random forests.

4.2

Recursive Subsampling

The previous subsection studies feature bagging where each base learner only subsamples features once. Recursive algorithms are more complex: each base learner performs multiple rounds of feature subsampling, and the subsampled features at each round are data-dependent. This recursive structure 12

makes the analysis more challenging because the effect from one feature removal can propagate in later updates. We model this propagation through a stochastic process {wt }t≥0 , where each state wt lies in a set Wt ⊆ H and H is equipped with the norm ∥ · ∥H . At iteration t, the update is wt = A(D, (wt−1 , ξt )) , where ξt denotes the randomness injected at step t. The admissible set Wt may depend on wt−1 , allowing the state space itself to be data-dependent and recursive. We initialize the process at w0 , the zero element of H, and assume that w0 ∈ Wt for all t. For a fixed dataset, define the update maps −j

S ξ (w) := A(D, (w, ξ)),

S ξ (w) := A(D, (w, ξ, j)).

Here the extra argument j means that the same seed ξ first generates the candidate feature set as in the original run, and then feature j is removed from that set if it appears. Using the same randomness ξ1:T , the original trajectory is A(T ) (D, ξ1:T ) = wT = S ξT ◦ · · · ◦ S ξ1 (w0 ),

(4.1)

whereas the trajectory after removing feature j is −j

−j

−j

A(T ) (D, (ξ1:T , j)) = wT = S ξT ◦ · · · ◦ S ξ1 (w0 ).

(4.2)

Independent copies of this recursive trajectory give the finite- and infinite-bagged outputs B

) A(T B (D) =

1 X (b) w , B b=1 T

) A(T ∞ (D) = Eξ1:T [wT ],

with their feature-removed counterparts denoted by a superscript − j. This framework covers random forward selection (Mentch and Zhou, 2020) and random forests (Breiman, 2001), both described in the following two subsections. −j −j To state the recursive guarantee, let ∆t := wt −wt . We isolate the part of the step-t perturbation that j −j −j −j comes from propagating the previous feature-removal discrepancy: pitht := E[S ξt (wt−1 ) − S ξt (wt−1 )]. The following one-step condition requires this propagated discrepancy to grow by at most a factor 1 + δt : for each step t, there exists δt ≥ 0 such that, for every j ∈ [d], j

pitht

H

h −j i ≤ (1 + δt ) E ∆t−1

H

.

(4.3)

This is a Lipschitz-type bound on the expected transition map, analogous to stepwise stability assumptions for iterative algorithms (Hardt et al., 2016; Lei and Ying, 2020).

Proposition 4.4 Let Assumption 1 hold. Suppose that for each recursion step t, the one-step condition (4.3) holds for every j ∈ [d]. Then d h − ji 2 MT q 1X E [wT ] − E wT ≤ , H d j=1 (d − 1)(1 − q)

13

(4.4)

q=0.05 q=0.1 q=0.3 q=0.5

0.14 0.12

0.4

0.1

0.3 2 2

0.08

0

w w0 22

q=0.05 q=0.1 q=0.3 q=0.5

0.06

0.2

0.04 0.1

0.02 0

1

5

10 T

15

0.0

20

1 Pd

(a) d

−j 2 j=1 ∥w − w ∥2

1

5

(b)

10 15 T 1 Pd :,− j 2 ∥2 j=1 ∥β − β d

20

Figure 4: Left: Feature instability of feature-bagged random forward selection and non-bagged baseline (forward selection). Solid curves with crosses show empirical averages over 100 repetitions for different q, dashed curves show the upper bound in (4.5), and the solid black curve gives the non-bagged baseline. Right: Parameter estimation difference of feature-bagged random forward selection for different q. where

 T  T  2 X  Y      MT :=   (1 + δt′ ) sup rad(Wt ) . ξ t=1 t′ =t+1

1:(t−1)

The bound preserves the same feature-subsampling factor q/((d − 1)(1 − q)) as the one-step1 result in Proposition 4.3. The new factor MT is the price of recursion: it accumulates the radii of the intermediate state spaces and the one-step inflation factors. In the simple case rad(Wt ) ≡ r and δt ≡ δ, we obtain MT = O(r2 T 2 ) when δ = 0, and MT = O(r2 (1 + δ)2T ) for fixed δ > 0. Proposition 4.4 is stated for infinite bagging to isolate the mechanism that improves feature stability. Section 4.3 shows that using B finite bagging rounds adds only an O(1/B) concentration error to the infinite-bagging instability term. Thus, for moderately large B, finite bagging inherits the same stability behavior up to this vanishing error. 4.2.1

Random Forward Selection

We first apply the recursive bound to random forward selection. The algorithm builds a feature set one coordinate at a time: at each step, it samples a fraction q of the available features and adds the sampled feature that gives the largest reduction in residual error. This places random forward selection within the class of data-dependent recursion algorithms covered by Proposition 4.4. To represent the recursive state, encode the selected variables by w = (w1 , . . . , wd ), where w j ∈ {0, 1} and w j = 1 means that feature j is included in the fitted linear model. In this specialization, Wt ⊆ {0, 1}d ⊂ H, with H = Rd and ∥ · ∥H = ∥ · ∥2 . For feature-subsampled ordinary least squares, w j = 1 exactly when j ∈ ν. Let Ew denote the diagonal matrix with diagonal w. The b-th subsampled ordinary least-squares estimator is
β (b) = Ew(b) XEw(b) † y. where the inverse is understood as the Moore–Penrose pseudoinverse in the overparameterized regime. Let β = X † y. We consider the orthogonal-design case, where X ⊤ X = nId . Then the feature-subsampled estimator 1

Here one step means that each base learner only subsamples features once.

14

Algorithm 1 Randomized Forward Selection for b ∈ [B] do (b) n Draw bootstrap sample D(b) = {(x(b) i , yi )}i=1 from original data D. Initialize empty active set w0 = 0. for t ∈ [T ] do Select subset of s features uniformly at random, denoted νt(b) . Select jt ∈ argmin j∈ν (b) ∥y (b) − X (b) Ewt−1 +e j (X (b) Ewt−1 +e j )† y (b) ∥22 , with deterministic tiet breaking; if every candidate is already selected, set e jt = 0. Update wt = wt−1 + e jt . Update coefficient estimates by fitting least squares on the selected coordinates: b(b) = argmin ∥y (b) − X (b) Ew β∥2 , Ewt β t 2

b(b) = 0. (I − Ewt )β

β

end for end for b = B−1 PB β b(b) . Compute final coefficient estimates β b=1 b b = X β. Compute predictions y reduces to β (b) = Ew(b) β and, with wB := B−1

PB

b=1 w

(b) ,

B B 1 X (b) 1 X β = E (b) β = EwB β. βB = B b=1 B b=1 w

Thus, in the orthogonal setting, the stability of the averaged estimator is governed by the selectionfrequency vector wB . Under uniform feature subsampling, each feature is selected with probability q, so E[wB ] = (q, . . . , q). Under random forward selection, wB is data-dependent and records how frequently each feature is selected, thereby capturing feature importance. We measure discrepancies between selection vectors by the ℓ2 distance. The following theorem bounds the resulting feature instability. The full randomized forward-selection procedure is shown in Algorithm 1. At step t of the b-th (b) bagging round, starting from the current selection vector wt−1 , the algorithm first subsamples a fraction q of the d available features according to the randomness ξt(b) , yielding a candidate set νt(b) . It then selects the feature that yields the greatest reduction in residual sum of squares: j(b) t = arg min min y − XEw(b) +e j β j∈νt(b) β

t−1

2 2

,

where e j denotes the j-th standard basis vector in Rd and Ew is the diagonal matrix that projects onto the coordinates selected by w. The selected variable is then added to the active set, (b) (b) A(wt−1 ) = wt(b) = wt−1 + e j(b) . t

If all indices in the candidate set νt(b) have already been selected, we set j(b) t = None and eNone = 0, so the selection vector remains unchanged.

Theorem 4.5 Let A be the random forward selection algorithm in Algorithm 1 with feature-subsampling ratio q < 1. Then rad(Wt ) ≤ 1 for t ∈ [T ]. Moreover, the one-step condition (4.3) in Proposition 4.4

15

holds with δt = 4−1 max{q(1 + q), q2 t}. When T ≲ q−1 , we have d h − ji 2 T2 1X q ≲ E [wT ] − E wT · . d j=1 d−1 1−q

(4.5)

The theorem specializes Proposition 4.4 to the selection path of random forward selection. The bound in (4.5) gives a stronger stability guarantee for smaller feature-subsampling ratios q, showing that more aggressive feature subsampling improves stability. Its quadratic dependence on T reflects the other side of the recursive analysis: instability can accumulate along a long selection path. For comparison, deterministic forward selection without bagging has feature instability d h − j i 2 2T 1X = E [wT ] − E wT , d j=1 d

whose derivation is deferred to the appendix. Hence, feature bagging is provably more stable whenever 2T T2 q ≳ · , d d−1 1−q

or equivalently,

T≲

1−q . q

√ √ This is consistent with the condition T ≲ q−1 ; for example, it is satisfied when T = o( d) and q = 1/ d. Figure 4 evaluates this prediction under the experimental settings detailed in Appendix S.3.1. Figure 4(a) shows that feature bagging consistently improves empirical feature stability, with smaller subsampling ratios q giving larger stability gains. Figure 4(b) shows a corresponding decrease in parameter estimation difference, suggesting that the stability gain does not come at the expense of estimation accuracy. When qT ≤ 1, the theoretical curve upper bounds the averaged empirical instability and certifies better stability than the non-bagged baseline. For small T , it also tracks the empirical curve closely, indicating that the bound captures the early-stage behavior of random forward selection. As T grows, the gap widens, consistent with the quadratic dependence on T and the conservative one-step condition used in the analysis. Overall, the experiment supports the same message as the theorem: feature bagging stabilizes random forward selection most clearly in the shallow-recursion regime. 4.2.2

Random Forests

Random forests provide the second application. In a full classification-and-regression tree, both the split feature and the split threshold are data-dependent. A decision tree’s partition can be encoded by two vectors: one that records the feature chosen at each internal node and another that stores the corresponding split threshold; see the feature and threshold attributes in Pedregosa et al. (2011). Both vectors have length 2T − 1, where T is the max depth parameter, i.e., the maximum number of splits from the root node to any leaf node. We analyze a simplified dyadic forest: each split is made at the midpoint of the selected coordinate, so the local cell containing a fixed test point x is determined only by the sequence of split features along the path to x. We further consider the max-edge version of the dyadic tree. The corresponding localpath procedure is summarized in Algorithm 2. Dyadic splitting is a simplified model of decision-tree partitions and is widely used in nonparametric statistical estimation (e.g. Blanchard et al., 2007; Perchet and Rigollet, 2013; Cai et al., 2023; Ma et al., 2023; Cai et al., 2024; Ma and Yang, 2024). 16

Figure 5: Two possibilities of dyadic tree partition. The black dot represents the target test sample x. The two partitions, though different in other areas, have the same representation w that corresponds to the grey area. Algorithm 2 Max-edge dyadic randomized forest, local path at x Let Er,k ∈ RT ×d be the matrix with one at entry (r, k) and zero elsewhere. for b ∈ [B] do (b) n Draw bootstrap sample D(b) = {(x(b) i , yi )}i=1 from original data D. Initialize local encoding w0 = 0T ×d and local cell L0 = [0, 1]d . for t ∈ [T ] do Select subset of s features uniformly at random, denoted νt(b) . k = ∥w k ∥ and M(w , ν (b) ) = argmin Define Nt−1 Nk . t−1 t t−1 1 k∈ν (b) t−1 t

Select jt ∈ argmaxk∈M(wt−1 ,ν (b) ) ∆k (wt−1 ; D(b) ), with deterministic tie-breaking. t Update wt = wt−1 + EN jt +1, jt . t−1 Split Lt−1 along feature jt at the midpoint and set Lt to the child cell containing x. end for Compute the prediction at x as the average of labels in LT . end for Average the B predictions.

Rather than analyzing the entire tree, we study a local property: the leaf node that contains a given query point x. This node is uniquely determined by the sequence of features used to split the path from the root to that leaf; see Figure 5. To encode the node, define the feature-weight matrix wt = (wt1 , . . . , wtd ) ∈ RT ×d . If, up to the t-th step, the leaf node that contains x has been split along feature j for t j times, then j the first t j entries of wt are one, while the remaining entries are zero. We write Ntk := ∥wtk ∥1 and − j,k −j Nt := ∥(wt )k ∥1 for the corresponding split counts in the full and feature-removed runs. At each local split t, draw a candidate feature set νt ⊂ [d] with |νt | = s and q = s/d uniformly without replacement. Define the candidate-set max-edge set M(wt , νt ) := argmin Ntk . k∈νt

The tree then chooses kt ∈ argmax ∆k (wt ; D), k∈M(wt ,νt )

17

−j

with deterministic tie-breaking. For the feature-removed run, use the coupled candidate set νt := νt \ { j}. −j If νt = ∅, the removed update returns the current encoding. Otherwise define −j

−j

− j,k

M− j (wt , νt ) := argmin Nt

,

−j k∈νt

and choose the feature in this set with the largest split decrease, using the same deterministic tie-breaking rule. In this dyadic case, only one vector with length 2T − 1 is enough to store the local tree path because the thresholds are uniquely defined. For instance, the x in Figure 5 has     0 1 1 0     w1 = 0 0 or 0 0 ,     0 0 0 0

  1 1   w2 = 0 0 ,   0 0

  1 1   w3 = 0 1 .   0 0

Moreover, given a fixed wt , the grey region in Figure 5 is the unique leaf node compatible with that encoding. Denote this region by Lwt . We now justify the choice of the encoding w. First, the encoding scheme is closely related to feature importances (Pedregosa et al., 2011), which measures a feature’s global importance by the proportion of splits that use that feature. Analogously, for a specific test point x, we can quantify the local importance of feature j by j

∥wT ∥1 j′ j′ =1 ∥wT ∥1

Pd

.

This point-specific view is aligned with recent work on individual variable importance, which moves beyond population-level summaries (Dai et al., 2025). We also demonstrate a direct relationship between the prediction at x and the matrix w through the following proposition. Consider data generated from the additive model Ey|x [y] = f ∗ (x) =

d X

f j∗ (x j ),

x ∼ Unif([0, 1]d ).

(4.6)

j=1

Proposition 4.6 Let Lw be the rectangle of the node associated with x under w, and let f w be the decision-tree predictor associated with w. Since we only care about x, we view all decision trees with the same w as equivalent. Based on the above definitions, there exists a map K such that   K(w) = Ey|X,w f w (x) . Moreover, the map K is affine: there exist a constant cx and a matrix Θx ∈ RT ×d such that K(w) = cx + ⟨Θx , w⟩. Consequently, K(w1 ) − K(w2 ) = ⟨Θx , w1 − w2 ⟩ ≲ ∥w1 − w2 ∥F .

18

Given Proposition 4.6, the random forest predictor is fully determined by the ensemble feature-weight matrix Eξ [w]:  h     i Ey|X f (x) = Eξ Ey|X,w f w (x) = Eξ [K(w)] = K Eξ [w] , where the last equality follows from the affine representation of K. Moreover, for the feature-removed forest, the affine representation gives  h i2   −j 2 Ey|X f (x) − f − j (x) = K(Eξ [wT ]) − K(Eξ [wT ]) D E −j 2 −j 2 = Θx , Eξ [wT ] − Eξ [wT ] ≲ Eξ [wT ] − Eξ [wT ] . F

Thus, it suffices to analyze Eξ [wT ] with respect to the Frobenius norm. We next give the instability bound of the encoding.

Theorem 4.7 Let A be the max-edge dyadic random forest algorithm described in Section 4.2.2. Then √ √ P rad(Wt ) ≤ 1, and the one-step condition (4.3) holds with δt = 2. Let AT = Tt=1 (1 + 2)T −t . Then d h − ji 2 q · A2T 1X . E [wT ] − E wT ≤ F d j=1 (d − 1)(1 − q)

(4.7)

Theorem 4.7 shows that recursive feature bagging stabilizes the local partition induced by a dyadic tree. Its dependence on the feature-subsampling ratio matches the one-step bound in Proposition 4.3: smaller q yields a stronger feature-stability guarantee. The factor A2T is the price of recursion, capturing how perturbations accumulate along the tree path. For comparison, a non-bagged depth-T dyadic tree has feature instability 2T/d, as derived in Appendix S.3.4. Thus, for shallow trees and sufficiently small q, feature bagging provides a strict stability improvement.

4.3

Finite Bagging

The preceding model-free guarantees focus on infinite bagging, whereas practical ensembles use a finite number B of bagging rounds. We now show that the infinite-bagging bounds remain representative when B is moderately large.

Proposition 4.8 (Feature instability of finite bagging) Let Assumption 1 hold. Then for a bagged algorithm A with B bagging rounds and outputs w(b) ∈ W, b ∈ [B], there holds   √ 2 √ X d B B d 6 3 rad(W)2 log d+1 1 X 1 X (b) 1 X (b),− j 3 2 δ w − w ≤ Eξ [w] − Eξ [w− j ] + d j=1 B b=1 B b=1 d j=1 B with probability 1 − δ.

The proof of Proposition 4.8 is given in Appendix S.3.8; it is a direct application of Hoeffding’s inequality and the union bound. The finite-bagging feature instability exceeds the infinite-bagging 19

B=1 B=4 B=16 B=64 B=256

Feature Instability (q = 0.05)

0.12

B=1 B=4 B=16 B=64 B=256

0.14 0.12 Feature Instability (q = 0.1)

0.14

0.1 0.08 0.06 0.04 0.02

0.1 0.08 0.06 0.04 0.02

0

1

5

10 T

15

0

20

1

5

3

3

2

2

1

1

0 1 2

B=1 B=4 B=16 B=64 B=256

3 4 5

1

2

3

4 T

5

15

20

(b) RFS, q = 0.1

(log) Feature Instability (q = 0.2)

(log) Feature Instability (q = 0.1)

(a) RFS, q = 0.05

10 T

6

0 1 2

B=1 B=4 B=16 B=64 B=256

3 4 5

7

1

(c) RF, q = 0.1

2

3

4 T

5

6

7

(d) RF, q = 0.2

Figure 6: Comparison of feature instability across different numbers of bagging rounds B. Each cross corresponds to 100 repetitions. The RFS and RF simulation settings are the same as those described in Appendix S.3.1. The black lines are instability bounds from Theorems 4.5 and 4.7. feature instability by a concentration term that decays linearly with the number of bagging rounds B. On the constant, the radius of the overall algorithm A is no larger than the sum of the suprema of rad(Wt ). Combining this finite-B concentration term with Proposition 4.4, it suffices to set B≳d·

1−q . q

(4.8)

In Figure 6, empirical feature stability improves as B increases. For larger q, fewer bagging rounds are needed to reach the theoretical stability regime predicted by (4.8). As observed by Soloff et al. (2024a), finite instance bagging is also harmless when B ≳ n(1 − p)/p. Though the results are dual, the feature side can be much more practical when d ≪ n. For instance, a typical tabular scale (n, d) ∼ (105 , 102 ) leads to a requirement of order 102 bagging rounds for feature stability, as opposed to order 105 rounds for instance stability.

5

Conclusion

We showed that feature bagging provides a principled route to algorithmic stabilization. By introducing FI as the feature-side counterpart to II, we made it possible to study how bagging improves stability along both the instance and feature axes. In linear regression, we derived sharp asymptotic characterizations of II and FI; in a model-free setting inspired by recursive feature subsampling in random forests, we obtained tight upper bounds for FI in early steps. Together, these results show that feature bagging consistently improves stability, that more aggressive feature subsampling can yield stronger stabilization, 20

and that a modest number of bagging rounds is enough to approach the infinite-bagging stability level.

Acknowledgments We thank the anonymous reviewers for their constructive comments and suggestions, which have helped improve the clarity and presentation of this work. The bulk of the work was carried out while YM was a visiting student at MBZUAI and University of Toronto. QS was partially supported by NSERC Grant RGPIN-2026-06888, Compute Canada, and MBZUAI.

Impact Statement This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.

References Adrian, M., Soloff, J. A., and Willett, R. (2024). Stabilizing black-box model selection with the inflated argmax. arXiv preprint arXiv:2410.18268. Biau, G., Cérou, F., and Guyader, A. (2010). On the rate of convergence of the bagged nearest neighbor estimate. Journal of Machine Learning Research, 11(22):687–712. Blanchard, G., Schäfer, C., Rozenholc, Y., and Müller, K.-R. (2007). Optimal dyadic decision trees. Machine Learning, 66:209–241. Boucheron, S., Lugosi, G., and Massart, P. (2013). Concentration Inequalities: A Nonasymptotic Theory of Independence. Oxford University Press, Oxford. Bousquet, O. and Elisseeff, A. (2002). Stability and generalization. Journal of Machine Learning Research, 2:499–526. Bousquet, O., Klochkov, Y., and Zhivotovskiy, N. (2020). Sharper bounds for uniformly stable algorithms. In Proceedings of the 33rd Conference on Learning Theory. Breiman, L. (1996). Bagging predictors. Machine Learning, 24(2):123–140. Breiman, L. (2001). Random forests. Machine Learning, 45(1):5–32. Cai, C., Cai, T. T., and Li, H. (2024). Transfer learning for contextual multi-armed bandits. The Annals of Statistics, 52(1):207–232. Cai, Y., Ma, Y., Dong, Y., and Yang, H. (2023). Extrapolated random tree for regression. In Proceedings of the 40th International Conference on Machine Learning. Cai, Y., Yang, H., Ma, Y., and Hang, H. (2025). Bagged regularized k-distances for anomaly detection. Journal of Machine Learning Research, 26(178):1–59. Chen, S. and Xu, C. (2023). Handling high-dimensional data with missing values by modern machine learning techniques. Journal of Applied Statistics, 50(3):786–804.

21

Chen, X., Zeng, Y., Yang, S., and Sun, Q. (2023). Sketched ridgeless linear regression: The role of downsampling. In Proceedings of the 40th International Conference on Machine Learning. Curth, A., Jeffares, A., and van der Schaar, M. (2024). Why do random forests work? Understanding tree ensembles as self-regularizing adaptive smoothers. arXiv preprint arXiv:2402.01502. Dai, G., Shao, L., and Chen, J. (2025). Moving beyond population variable importance: concept, theory and applications of individual variable importance. Journal of the Royal Statistical Society Series B: Statistical Methodology, 87(3):816–832. Elisseeff, A., Evgeniou, T., and Pontil, M. (2005). Stability of randomized learning algorithms. Journal of Machine Learning Research, 6:55–79. Feldman, V. and Vondrak, J. (2018). Generalization bounds for uniformly stable algorithms. In Advances in Neural Information Processing Systems. Feldman, V. and Vondrak, J. (2019). High probability generalization bounds for uniformly stable algorithms with nearly optimal rate. In Proceedings of the 32nd Conference on Learning Theory. Friedman, J. H. (1991). Multivariate adaptive regression splines. The Annals of Statistics, 19(1):1–67. Geurts, P., Ernst, D., and Wehenkel, L. (2006). Extremely randomized trees. Machine Learning, 63(1):3–42. Haff, L. (1979). An identity for the wishart distribution with applications. Journal of Multivariate Analysis, 9(4):531–544. Hansen, L. K. and Salamon, P. (1990). Neural network ensembles. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(10):993–1001. Hardt, M., Recht, B., and Singer, Y. (2016). Train faster, generalize better: Stability of stochastic gradient descent. In Proceedings of the 33rd International Conference on Machine Learning. Harrington, E. F. (2003). Online ranking/collaborative filtering using the perceptron algorithm. In Proceedings of the 20th International Conference on Machine Learning. Hastie, T., Montanari, A., Rosset, S., and Tibshirani, R. J. (2022). Surprises in high-dimensional ridgeless least squares interpolation. The Annals of Statistics, 50(2):949–986. Hou, H., Huang, W., and Zhang, Z. (2026). Nonparametric quantile regression and uniform inference with unknown error distribution. Journal of Business & Economic Statistics, 44(1):1–12. Jeng, X. J., Hu, Y., Sun, Q., and Li, Y. (2024). Weak signal inclusion under dependence and applications in genome-wide association study. The Annals of Applied Statistics, 18(1):841–857. Kutin, S. and Niyogi, P. (2002). Almost-everywhere algorithmic stability and generalization error. In Proceedings of the 18th Conference on Uncertainty in Artificial Intelligence. Lakshminarayanan, B., Pritzel, A., and Blundell, C. (2017). Simple and scalable predictive uncertainty estimation using deep ensembles. In Advances in Neural Information Processing Systems. Lei, Y. and Ying, Y. (2020). Fine-grained analysis of stability and generalization for stochastic gradient descent. In Proceedings of the 37th International Conference on Machine Learning.

22

LeJeune, D., Javadi, H., and Baraniuk, R. (2020). The implicit regularization of ordinary least squares ensembles. In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics. Liang, R., Soloff, J. A., Barber, R. F., and Willett, R. (2025). Assumption-free stability for ranking problems. arXiv preprint arXiv:2506.02257. Little, R. J. A. and Rubin, D. B. (2019). Statistical Analysis with Missing Data. John Wiley & Sons, Hoboken, NJ, 3rd edition. Liu, B. and Mazumder, R. (2025). Randomization can reduce both bias and variance: A case study in random forests. Journal of Machine Learning Research, 26(150):1–49. Liu, T., Lugosi, G., Neu, G., and Tao, D. (2017). Algorithmic stability and hypothesis complexity. In Proceedings of the 34th International Conference on Machine Learning. Ma, Y., Jia, K., and Yang, H. (2025). Locally private estimation with public features. In Proceedings of the 28th International Conference on Artificial Intelligence and Statistics. Ma, Y. and Yang, H. (2024). Optimal locally private nonparametric classification with public data. Journal of Machine Learning Research, 25(167):1–62. Ma, Y., Zhang, H., Cai, Y., and Yang, H. (2023). Decision tree for locally private estimation with public data. In Advances in Neural Information Processing Systems. McWilliams, B., Krummenacher, G., Lucic, M., and Buhmann, J. M. (2014). Fast and robust least squares estimation in corrupted linear models. In Advances in Neural Information Processing Systems. Meinshausen, N. and Bühlmann, P. (2010). Stability selection. Journal of the Royal Statistical Society Series B: Statistical Methodology, 72(4):417–473. Mentch, L. and Zhou, S. (2020). Randomization as regularization: A degrees of freedom explanation for random forest success. Journal of Machine Learning Research, 21(171):1–36. Mukherjee, S., Niyogi, P., Poggio, T., and Rifkin, R. (2006). Learning theory: stability is sufficient for generalization and necessary and sufficient for consistency of empirical risk minimization. Advances in Computational Mathematics, 25:161–193. Murdoch, W. J., Singh, C., Kumbier, K., Abbasi-Asl, R., and Yu, B. (2019). Definitions, methods, and applications in interpretable machine learning. Proceedings of the National Academy of Sciences, 116(44):22071–22080. Nogueira, S., Sechidis, K., and Brown, G. (2018). On the stability of feature selection algorithms. Journal of Machine Learning Research, 18(174):1–54. Oza, N. C. and Russell, S. J. (2001). Online bagging and boosting. In Proceedings of the 8th International Workshop on Artificial Intelligence and Statistics. Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825–2830.

23

Perchet, V. and Rigollet, P. (2013). The multi-armed bandit problem with covariates. The Annals of Statistics, 41(2):693 – 721. Perrone, M. P. and Cooper, L. N. (1995). When networks disagree: Ensemble methods for hybrid neural networks. In How We Learn; How We Remember: Toward an Understanding of Brain and Neural Systems: Selected Papers of Leon N. Cooper, volume 10 of World Scientific Series in 20th Century Physics, pages 342–358. World Scientific, Singapore. Poudel, P., Chhetri, A., Gyawali, P. K., Leontidis, G., and Bhattarai, B. (2025). Multimodal federated learning with missing modalities through feature imputation network. In Medical Image Understanding and Analysis. Shalev-Shwartz, S., Shamir, O., Srebro, N., and Sridharan, K. (2010). Learnability, stability and uniform convergence. Journal of Machine Learning Research, 11:2635–2670. Shen, Z. and Xiu, D. (2025). Can machines learn weak signals? NBER Working Paper 33421, National Bureau of Economic Research. Soloff, J. A., Barber, R. F., and Willett, R. (2024a). Bagging provides assumption-free stability. Journal of Machine Learning Research, 25(131):1–35. Soloff, J. A., Barber, R. F., and Willett, R. (2024b). Building a stable classifier with the inflated argmax. In Advances in Neural Information Processing Systems. Soloff, J. A., Barber, R. F., and Willett, R. (2024c). Stability via resampling: statistical problems beyond the real line. arXiv preprint arXiv:2405.09511. Su, B., Sun, Q., Yang, X., and Zhao, B. (2024). The exact risks of reference panel-based regularized estimators. arXiv preprint arXiv:2401.11359. Sutton, C., Sindelar, M., and McCallum, A. (2006). Reducing weight undertraining in structured discriminative learning. In Proceedings of the Human Language Technology Conference of the NAACL, Main Conference. Wang, Y., Wu, H., and Nettleton, D. (2023). Stability of random forests and coverage of random-forest prediction intervals. In Advances in Neural Information Processing Systems. Wu, M. and Sun, Q. (2025). Ensemble linear interpolators: The role of ensembling. SIAM Journal on Mathematics of Data Science, 7(2):438–467. Xing, Y., Song, Q., and Cheng, G. (2021). On the algorithmic stability of adversarial training. In Advances in Neural Information Processing Systems. Yu, B. and Kumbier, K. (2020). Veridical data science. Proceedings of the National Academy of Sciences, 117(8):3920–3929. Zaidi, S., Zela, A., Elsken, T., Holmes, C. C., Hutter, F., and Teh, Y. W. (2021). Neural ensemble search for uncertainty estimation and dataset shift. In Advances in Neural Information Processing Systems. Zhang, Y., Ji, W., and Bradic, J. (2024). Adaptive split balancing for optimal random forest. arXiv preprint arXiv:2402.11228. Zhou, S., Lei, Y., and Kabán, A. (2023). Toward better PAC-Bayes bounds for uniformly stable algorithms. In Advances in Neural Information Processing Systems. 24

Zhou, S. and Mentch, L. (2023). Trees, forests, chickens, and eggs: when and why to prune trees in a random forest. Statistical Analysis and Data Mining: The ASA Data Science Journal, 16(1):45–64.

25

Appendix Contents S.1 Experiment Details of PMI between Feature Instability and Generalization S.1.1 Data Generating Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S.1.2 Model Training and Instability Estimation . . . . . . . . . . . . . . . . . . . . . . . . S.1.3 Statistical Dependence Measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S.1.4 More Explanation about Partial Mutual Information . . . . . . . . . . . . . . . . . . . S.1.5 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S.1.6 Additional Robustness Checks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26 26 27 27 27 28 29

S.2 Contents Related to Stability of Bagged Linear Regression S.2.1 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S.2.2 Useful Lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S.2.3 Proofs for Appendix S.2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31 31 34 43

S.3 Contents Related to Model Free Stability Guarantees S.3.1 Additional Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S.3.2 Proofs of Results in Section 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S.3.3 Related Contexts and Proofs of Results in Section 4.2 . . . . . . . . . . . . . . . . . . S.3.4 Derivation of Feature Instability for Non-bagged RFS and RF . . . . . . . . . . . . . . S.3.5 Derivation of the Affine Random Forest Encoding Map . . . . . . . . . . . . . . . . . S.3.6 Proof of the Random Forward Selection Instability Bound . . . . . . . . . . . . . . . S.3.7 Proof of the Random Forest Instability Bound . . . . . . . . . . . . . . . . . . . . . . S.3.8 Proof of the Finite-Bagging Instability Bound . . . . . . . . . . . . . . . . . . . . . .

49 49 51 53 55 55 56 61 63

S.1 Experiment Details of PMI between Feature Instability and Generalization This appendix provides full experimental details and quantitative analyses supporting Section 2.2.

S.1.1

Data Generating Process

We adopt the MARSadd setting from Friedman (1991), a standard benchmark for evaluating random forests (Mentch and Zhou, 2020; Curth et al., 2024). For each experiment, we independently sample X ∈ [0, 1]n×d from the uniform distribution and generate responses according to f (x) = (1 − ν) fweak (x) + ν fstrong (x), where

d X 1 xj fweak (x) = √ d − 4 j=5

represents a linear but low signal-to-noise component, in which the signal is distributed across many features. Let 4 fstrong (x) = 0.1e4x1 + + 3x3 + 2x4 1 + e−20x2 +10 26

represents a concentrated but nonlinear signal. Independent Gaussian noise ϵ ∼ N(0, σ2 ) is added to form y = f (x) + ϵ. An independent test set of equal size is generated using the same procedure. The mixing parameter ν ∈ [0, 1] controls the dominance of nonlinear structure.

S.1.2

Model Training and Instability Estimation

The random forest models were trained over a predefined hyperparameter grid to systematically explore the effects of model complexity and ensemble diversity. The number of estimators was fixed at nestimators = 256. To control tree diversity, the maximum feature ratio (max features) and the maximum sample ratio (max samples) were varied over {0.1, 0.3, 0.5, 0.7}. Model complexity was adjusted by varying the maximum tree depth (max depth) from 1 to 7, while the minimum number of samples per leaf (min samples leaf) was fixed at 1. This design enables a systematic investigation of model complexity along both the feature-subsampling and depth dimensions. In particular, max features and max samples control the diversity of individual trees, whereas max depth governs the bias-variance trade-off of the ensemble. Fixing the number of estimators improves comparability across runs. For each configuration of hyperparameters {max features, max samples, max depth}, we train a RandomForestRegressor on the training data and compute both training and test mean squared errors (MSE). Both instance and feature instability metrics are then measured as in (1.1) and (1.2), by retraining the model on the data with the removed feature or instance. All configurations are repeated 500 times with different random seeds, ensuring robust Monte Carlo estimates. For each unique hyperparameter configuration, we compute the mean instability scores, mean squared errors, and define the generalization gap as Gap = MSEtest − MSEtrain . Results are filtered for fixed conditions (n = 500, d = 10, σ = 1). The MSEtest is computed using additional independent test samples, and MSEtrain is the in-sample training error.

S.1.3

Statistical Dependence Measures

We then investigate the statistical dependence between instability measures and the generalization gap. Three types of associations are evaluated: • Linear dependence: R2 values from linear regression models of the form Gap ∼ S feature , Gap ∼ S instance , and Gap ∼ (S feature , S instance ). • Mutual information (MI): Nonlinear associations are estimated using the mutual info regression function from scikit-learn. • Partial mutual information (PMI): To isolate the unique contribution of each instability measure, we apply a residual-based approach using random forest regressions: we regress out the conditioning variable (e.g., S instance ) from both the predictor and target, and compute the mutual information between the residuals via the MINE estimator. This yields estimates of PMI(S feature ; Gap | S instance ) and PMI(S instance ; Gap | S feature ).

S.1.4

More Explanation about Partial Mutual Information

We provide an overview of partial mutual information used in the experiments.

27

Mutual information (MI) quantifies the amount of shared information between two random variables X and Y. The partial mutual information (PMI) aims to measure the direct association between X and Y while removing the indirect influence of Z, typically by regressing out Z from both X and Y and then e and Y e denote the residuals computing the mutual information between the residuals. Formally, if X obtained from regressing X and Y on Z, respectively, then e Y). e PMI(X; Y | Z) = I(X;

(S.1)

Thus, PMI can be viewed as a nonparametric extension of the partial correlation concept to informationtheoretic measures. Consider three Gaussian random variables satisfying X = g1 (Z) + ϵX ,

(S.2)

Y = g2 (Z) + ϵY ,

(S.3)

where ϵX and ϵY are independent Gaussian noises and Z is a common latent variable. In this case, X and Y are dependent due to their shared dependence on Z, yet they contain no direct interaction. If we compute the partial mutual information PMI(X; Y | Z), it is zero because, after regressing out Z, the e and Y e are independent. However, if we slightly modify the model by introducing a direct residuals X coupling: Y = g1 (Z) + g2 (X) + ϵY , (S.4) then PMI(X; Y | Z) > 0, indicating a direct dependence between X and Y. If we treat the generalization gap as Y, and take (X, Z) to be either (FI, II) or (II, FI), then both PMIs should be positive when both measures contain generalization-relevant information.

S.1.5

Analysis

We investigate how FI and II jointly and individually explain variation in the generalization gap Y, under two configurations of the hyperparameter ν (ν = 1 and ν = 0). We report both deterministic correlations (via R2 ) and information-theoretic dependencies (via mutual information and partial mutual information). Case 1: ν = 1. When ν = 1, the individual explanatory powers of II and FI are moderate: R2 (II → Y) = 0.143 and R2 (FI → Y) = 0.120. However, their joint contribution increases substantially, achieving R2 (FI + II → Y) = 0.345. This suggests a complementary effect between II and FI: neither measure alone suffices to capture the variability of the generalization gap, but together they explain a larger fraction of variance. The mutual information (MI) values support this observation. The joint dependency between (FI, II) and Y is substantial (I((FI, II); Y) = 1.345), with asymmetric contributions: PMI(FI; Y | II) = 0.469, PMI(II; Y | FI) = 0.439, and mutual dependence between FI and II (I(FI; II) = 0.558). These results indicate that both measures convey overlapping but not redundant information about generalization. FI contributes slightly more conditional information about the generalization gap once II is known. Case 2: ν = 0. At ν = 0, the relationships become markedly stronger. R2 (II → Y) = 0.231, R2 (FI → Y) = 0.546, and R2 (FI + II → Y) = 0.675. The dominant role of FI indicates that feature-level perturbations are the primary source of generalization-relevant variation, while II contributes additional but smaller explanatory power. Information-theoretic results are consistent: I((FI, II); Y) = 1.062 (with FI dominating), PMI(FI; Y | II) = 0.498, PMI(II; Y | FI) = 0.402, and I(FI; II) = 0.748. Compared with the ν = 1 case, the 28

interdependence between FI and II increases, suggesting that when regularization or noise is reduced, the two measures become more correlated. Partial mutual information values also reveal that the unique contribution of each variable to generalization decreases slightly, consistent with the growing redundancy between FI and II. Interpretation Overall, these findings highlight two trends: (1) complementarity, where FI and II jointly enhance generalization predictability, especially when they are less correlated (as in ν = 1), and (2) a redundancy-dominance shift, where for smaller ν, FI becomes the dominant predictor while its correlation with II strengthens, reducing the independent information each measure provides. This suggests that different regimes of model regularization or noise change how II and FI interact with generalization: from complementary (when ν = 1) to partially redundant but stronger overall predictors (when ν = 0).

S.1.6

Additional Robustness Checks

The preceding analysis is a focused synthetic random-forest study at two representative values of ν. The robustness study here strictly broadens that analysis along four axes: six synthetic signal regimes ν ∈ {0, 0.2, 0.4, 0.6, 0.8, 1}, four base estimator families (RF, DT, AdaBoost, and GBRT), four PMI residualizers (RF, DT, AdaBoost, and GBRT), and four completed real tabular regression benchmarks (diabetes, abalone, cpu act, and house prices). The base estimator grids contain 105, 35, 63, and 63 hyperparameter settings, respectively, and each reported summary uses 100 bootstrap replicates, each subsampling one third of the corresponding grid. Tables S.2–S.5 report the conditional contribution of FI after accounting for II. The entries are positive across all displayed synthetic and real-data settings under both PMI(FI; gap | II) and incremental ∆R2 . Across the 96 synthetic estimator/residualizer/ν groups and the 64 real-data estimator/residualizer/dataset groups, the bootstrap positive rate is 1.0 for every group. For ∆R2 , the values repeat across residualizer columns because the residualizer is only used in the PMI calculation. Thus, the conclusion from the focused analysis above persists under a broader estimator grid, multiple dependence residualizers, and real tabular data.

29

Table S.2: Synthetic PMI(FI; gap | II) across estimator families and PMI residualizers (6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Synthetic setting

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

ν = 0.0 ν = 0.2 ν = 0.4 ν = 0.6 ν = 0.8 ν = 1.0

0.380 0.251 0.238 0.228 0.223 0.250

0.633 0.496 0.378 0.373 0.322 0.303

0.485 0.454 0.347 0.522 0.582 0.458

0.376 0.443 0.375 0.306 0.368 0.414

0.417 0.325 0.252 0.299 0.236 0.291

0.841 0.723 0.555 0.495 0.404 0.378

0.598 0.590 0.428 0.636 0.682 0.585

0.696 0.644 0.525 0.424 0.550 0.537

0.375 0.263 0.239 0.249 0.229 0.265

0.582 0.441 0.342 0.344 0.307 0.291

0.460 0.455 0.349 0.512 0.595 0.488

0.350 0.404 0.365 0.281 0.348 0.397

0.337 0.239 0.223 0.226 0.222 0.224

0.486 0.391 0.335 0.352 0.294 0.302

0.432 0.378 0.303 0.427 0.477 0.379

0.343 0.360 0.340 0.255 0.318 0.381

Table S.3: Synthetic ∆R2 across estimator families and PMI residualizers (6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Synthetic setting

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

ν = 0.0 ν = 0.2 ν = 0.4 ν = 0.6 ν = 0.8 ν = 1.0

0.152 0.064 0.014 0.028 0.011 0.016

0.485 0.402 0.320 0.310 0.133 0.095

0.246 0.300 0.204 0.510 0.496 0.316

0.231 0.273 0.334 0.332 0.420 0.454

0.152 0.064 0.014 0.028 0.011 0.016

0.485 0.402 0.320 0.310 0.133 0.095

0.246 0.300 0.204 0.510 0.496 0.316

0.231 0.273 0.334 0.332 0.420 0.454

0.152 0.064 0.014 0.028 0.011 0.016

0.485 0.402 0.320 0.310 0.133 0.095

0.246 0.300 0.204 0.510 0.496 0.316

0.231 0.273 0.334 0.332 0.420 0.454

0.152 0.064 0.014 0.028 0.011 0.016

0.485 0.402 0.320 0.310 0.133 0.095

0.246 0.300 0.204 0.510 0.496 0.316

0.231 0.273 0.334 0.332 0.420 0.454

Table S.4: Real-data PMI(FI; gap | II) across estimator families and PMI residualizers (6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Real dataset

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

Abalone CPU Act Diabetes House prices

0.219 0.264 0.358 0.216

0.334 0.387 0.368 0.310

0.384 0.278 0.549 0.404

0.292 0.283 0.437 0.365

0.232 0.330 0.412 0.222

0.387 0.438 0.444 0.313

0.452 0.296 0.597 0.399

0.573 0.361 0.538 0.436

0.227 0.279 0.364 0.219

0.316 0.372 0.332 0.288

0.402 0.281 0.524 0.373

0.283 0.266 0.404 0.355

0.216 0.252 0.339 0.213

0.390 0.471 0.328 0.294

0.363 0.293 0.505 0.368

0.276 0.269 0.398 0.347

Table S.5: Real-data ∆R2 across estimator families and PMI residualizers (6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Real dataset

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

Abalone CPU Act Diabetes House prices

0.017 0.087 0.100 0.016

0.212 0.202 0.217 0.047

0.295 0.028 0.469 0.238

0.114 0.027 0.314 0.048

0.017 0.087 0.100 0.016

0.212 0.202 0.217 0.047

0.295 0.028 0.469 0.238

0.114 0.027 0.314 0.048

0.017 0.087 0.100 0.016

0.212 0.202 0.217 0.047

0.295 0.028 0.469 0.238

0.114 0.027 0.314 0.048

0.017 0.087 0.100 0.016

0.212 0.202 0.217 0.047

0.295 0.028 0.469 0.238

0.114 0.027 0.314 0.048

30

We also performed a noise-control check by replacing FI with Gaussian and permuted-noise baselines and rerunning the same conditional-dependence analysis. Tables S.6–S.9 report the differences between the FI-based statistic and the corresponding noise-based statistic, averaged over the two noise controls. The incremental ∆R2 advantage of FI over the noise baseline is positive in all 320 noise-control aggregate entries. The PMI(FI; gap | II) advantage over PMI(noise; gap | II) is positive in 319 out of 320 entries; the only negative entry is −7.7 × 10−4 , which is numerically negligible. Thus, the observed conditional signal is not explained by simply adding another random covariate.

S.2

Contents Related to Stability of Bagged Linear Regression

S.2.1

Error Analysis

In the general setting, we may remove both instance i and feature j simultaneously. The updated sketching matrices Ub′ and Vb′ reflect these removals as defined above. The corresponding modified bagged estimator is: B 1 X ′ ′⊤ ′ β := V (U XVb′ )† Ub′⊤ y. B b=1 b b For notational simplicity, we use Ub′ and Vb′ throughout to denote the sketching matrices under the appropriate removal operation. If only instance subsampling is applied (i.e., no feature subsampling), then Vb′ = Vb for all b; conversely, if only feature sampling is used, we have Ub′ = Ub . For notational simplicity, we define the weighting matrices g(X) =

B B 1 X ′ ⊤′ 1X Vb (Ub⊤ XVb )† Ub⊤ , g′ (X) = V (U XVb′ )† Ub⊤′ , B b=1 B b=1 b b

so that β = g(X)y, β ′ = g′ (X)y. We decompose the instability into   h

i 2 E β − β ′ 2 =EX,y,U ,V y ⊤ g(X) − g′ (X) ⊤ g(X) − g′ (X) y h 

i = σ2 EX,U ,V tr g(X) − g′ (X) ⊤ g(X) − g′ (X) | {z } Variance

h 

i 1 + EX,U ,V tr X ⊤ g(X) − g′ (X) ⊤ g(X) − g′ (X) X . d | {z } Bias

The two terms correspond to the impact of removal on the variance and bias terms, respectively. The following propositions give the exact values of the variance and the bias term.

Proposition S.2.1 (Variance term of instability) Assume the above assumptions hold. Then the expectation of the variance term is ! V= h 
⊤
i (B − 1)∆V, ′ ′ 2 ∆ σ EX,U ,V tr g(X) − g (X) g(X) − g (X) = σ + , B B 2

31

Table S.6: Synthetic noise-control: average PMI difference PMI(FI; gap | II)−PMI(noise; gap | II) across estimator families and PMI residualizers (averaged over Gaussian and permuted noise; 6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Synthetic setting

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

ν = 0.0 ν = 0.2 ν = 0.4 ν = 0.6 ν = 0.8 ν = 1.0

0.167 0.042 0.028 0.020 0.008 0.042

0.371 0.239 0.115 0.095 0.050 0.039

0.241 0.206 0.103 0.262 0.325 0.221

0.116 0.198 0.134 0.049 0.135 0.162

0.201 0.119 0.040 0.079 0.026 0.079

0.565 0.454 0.280 0.241 0.139 0.098

0.357 0.344 0.181 0.382 0.424 0.338

0.452 0.394 0.270 0.173 0.323 0.298

0.170 0.045 0.027 0.042 0.021 0.047

0.325 0.185 0.100 0.094 0.039 0.019

0.216 0.204 0.104 0.247 0.344 0.241

0.101 0.160 0.122 0.037 0.120 0.150

0.118 0.027 0.007 0.011 0.006 0.011

0.191 0.128 0.074 0.063 0.029 0.029

0.188 0.112 0.053 0.174 0.209 0.131

0.096 0.120 0.099 0.012 0.086 0.123

Table S.7: Synthetic noise-control: average ∆R2 difference ∆R2 (FI) − ∆R2 (noise) across estimator families and PMI residualizers (averaged over Gaussian and permuted noise; 6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Synthetic setting

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

ν = 0.0 ν = 0.2 ν = 0.4 ν = 0.6 ν = 0.8 ν = 1.0

0.152 0.062 0.011 0.027 0.011 0.013

0.477 0.408 0.294 0.292 0.137 0.083

0.240 0.306 0.194 0.516 0.448 0.316

0.235 0.254 0.333 0.335 0.422 0.463

0.152 0.061 0.011 0.027 0.011 0.013

0.475 0.407 0.295 0.290 0.141 0.085

0.241 0.307 0.195 0.517 0.447 0.318

0.234 0.255 0.334 0.335 0.423 0.465

0.152 0.061 0.011 0.027 0.011 0.013

0.475 0.407 0.295 0.289 0.142 0.082

0.240 0.307 0.195 0.515 0.447 0.319

0.235 0.253 0.333 0.335 0.423 0.463

0.153 0.061 0.011 0.027 0.011 0.013

0.476 0.411 0.296 0.290 0.140 0.085

0.240 0.305 0.195 0.516 0.448 0.317

0.234 0.254 0.334 0.336 0.423 0.463

Table S.8: Real-data noise-control: average PMI difference PMI(FI; gap | II)−PMI(noise; gap | II) across estimator families and PMI residualizers (averaged over Gaussian and permuted noise; 6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Real dataset

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

Abalone CPU Act Diabetes House prices

0.007 0.052 0.140 0.005

0.068 0.114 0.115 0.045

0.153 0.032 0.294 0.150

0.051 0.040 0.200 0.120

0.022 0.119 0.192 0.004

0.096 0.149 0.172 0.043

0.224 0.045 0.337 0.143

0.325 0.119 0.302 0.194

0.020 0.067 0.151 0.009

0.047 0.093 0.062 0.030

0.167 0.027 0.272 0.131

0.043 0.021 0.149 0.108

0.000 0.039 0.118 0.001

0.025 0.039 0.063 0.007

0.116 0.042 0.243 0.122

0.033 0.017 0.152 0.104

Table S.9: Real-data noise-control: average ∆R2 difference ∆R2 (FI) − ∆R2 (noise) across estimator families and PMI residualizers (averaged over Gaussian and permuted noise; 6 seeds/config, 100 bootstrap repetitions). RF residualizer

DT residualizer

AdaBoost residualizer

GBRT residualizer

Real dataset

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

RF

DT

Ada

GBRT

Abalone CPU Act Diabetes House prices

0.011 0.083 0.092 0.010

0.184 0.173 0.201 0.035

0.312 0.018 0.430 0.238

0.103 0.020 0.327 0.046

0.013 0.084 0.092 0.011

0.186 0.171 0.200 0.035

0.312 0.017 0.430 0.236

0.103 0.021 0.328 0.046

0.012 0.082 0.092 0.010

0.184 0.175 0.201 0.033

0.311 0.019 0.431 0.238

0.103 0.021 0.327 0.047

0.013 0.082 0.091 0.009

0.183 0.172 0.199 0.036

0.311 0.017 0.431 0.239

0.103 0.021 0.327 0.046

Table S.10: Noise-control consistency against the Gaussian baseline. “FI wins on” is the average, within each scope, of the bootstrap win rate for a fixed completed row; that row-level win rate is the fraction of bootstrap replicates in which the FI-based statistic exceeds its Gaussian-noise counterpart. Scope

mean FI wins mean FI wins ∆R2 diff on ∆R2 PMI diff on PMI

Real Synthetic subset

0.143 0.252

0.935 0.972

32

0.103 0.155

0.852 0.907

where  2γq(1−p)     (p−γq)(1−γq)    2γpq(1−q)     (p−γq)(p−γq2 ) V= lim ∆ :=   2γq(p−p2 )  n,d→∞    (γq−p)(γq−p2 )      2γq(1−p)

if γq < p and i-th instance removed

 γq2     (1−γq2 )2    q2     (1−γq2 )2 V, lim n∆ :=   γp2  n,d→∞    (γ−p2 )2   2    p

if γq < p and i-th instance removed

(γq−p)(γ−p)

if γq < p and j-th feature removed if γq > p and i-th instance removed if γq > p and j-th feature removed

and

(γ−p2 )2

if γq < p and j-th feature removed if γq > p and i-th instance removed if γq > p and j-th feature removed .

Proposition S.2.2 (Bias term of instability) Assume the above assumptions hold. Then, the expectation of the bias term is h 

i B − 1 B, 1 B= 1 EX,U ,V tr X ⊤ g(X) − g′ (X) ⊤ g(X) − g′ (X) X = ∆ + ∆ , d B B where  2γq(1−q)(1−p)     (p−γq)(1−γq)    2γpq(1−q)2 2pq(1−q)     (p−γq)(p−γq2 ) + p−γq2 B= lim ∆ :=   2γpq(γ−p)(1−p) 2(1−p)p2  n,d→∞   2 ) − γ(γq−p2 )  γ(γq−p)(γq−p      2p(1−q) γq−p

if γq < p and i-th instance removed if γq < p and j-th feature removed if γq > p and i-th instance removed if γq > p and j-th feature removed .

and  γq2 (1−q)2     (1−γq2 )2   2 (1+γ−2γq)  q     γ(1−γq2 )2 B, lim n∆ :=   p2 (γ−p)2  n,d→∞    γ(γ−p2 )2   2    p (1−2p+γ) (γ−p2 )2

if γq < p and i-th instance removed if γq < p and j-th feature removed if γq > p and i-th instance removed if γq > p and j-th feature removed .

Proof of Theorem 3.1. Theorem 3.1, including all entries of Table 1, follows by applying Propositions S.2.1 and S.2.2 to the two perturbations ℓ = −i and ℓ = − j, and then collecting the corresponding variance and bias terms. □

33

S.2.2

Useful Lemmas

Lemma S.2.3 (Expectation of matrix under MP-law) Let X be a n × d matrix with each row sampled from N(0, Σ), where Σ is invertible. Then, there holds h i EX X ⊤ X = nΣ,

h i EX (X ⊤ X)−1 =

1 Σ−1 . n−d−1

Proof of S.2.3. The conclusions follow from standard properties of inverse Wishart distribution, see e.g. Haff (1979). □

Lemma S.2.4 (Lemma A.1 of LeJeune et al. (2020)) Let V and V c be the selection matrix corresponds to ν and ν c , respectively. Then for any random matrix X ∈ Rn×d whose entries are sampled i.i.d. from the standard Gaussian distribution, we have the following holds true h i h i EXV c V ⊤ X † = (XS)† if n > d, and EU c⊤ X X † T = (U ⊤ X)† if n < d. Moreover, if h(XS) and XV c are independent, there holds h i h i EXV c V c⊤ X ⊤ h(XS)V ⊤ X † = 0 if n > d, and EU c⊤ X X ⊤ U c h(U ⊤ X)X † T = 0 if n < d.

h i h i |µ | 1 The assumption ET U1⊤ U1 = ET U2⊤ U2 = tr(Λ) Λ appears less natural, but is important to unify ′ the results for T and T . Lemma S.2.5 (Variance) For µ1 , µ2 , ν1 , ν2 sampled independently, let U1 , U2 , V1 , V2 be their selection matrix. Λ1 , Λ2 andh Θ1 , Θ2iare diagonal matrices with 1s and 0s Then let h Suppose i h i only h on their i diagonal. |µ1 | |µ2 | |ν1 | |ν2 | ⊤ ⊤ ⊤ ⊤ ET U1 U1 = tr(Λ1 ) Λ1 , ET U2 U2 = tr(Λ2 ) Λ2 , ES V1 V1 = tr(Θ1 ) Θ1 , ES V2 V2 = tr(Θ2 ) Θ2 . Assume we always have Λ1 Λ2 ∈ {Λ1 , Λ2 } and Θ1 Θ ∈ {Θ1 , Θ2 }. Define Λ∧ = Λ1 Λ2 , Λ∨ = Λ1 + Λ2 − Λ∧ , and Θ∧ = Θ1 Θ2 , Θ∨ = Θ1 + Θ2 − Θ∧ . Then there holds  |ν1 ∩ν2 |  h  i    tr(Λ∨ )−|ν1 ∩ν2 |−1 ⊤ †⊤ ⊤ ⊤ † ⊤ EX,S,T tr U1 (U1 XV1 ) V1 V2 (U2 XV2 ) U2 =   1 ∩µ2 |   tr(Θ∨|µ )−|µ ∩µ |−1 1

2

if |ν1 | < |µ1 |, |ν2 | < |µ2 |, if |ν1 | > |µ1 |, |ν2 | > |µ2 |. (S.1)

Moreover, even if µ1 and µ2 are not independent, and ν1 and ν2 are not independent, yet we know µ1 ⊂ µ2 , ν1 ⊂ ν2 , the quantity becomes  |ν1 |  h  i    |µ2 |−|ν1 |−1 if |ν1 | < |µ1 |, |ν2 | < |µ2 |, EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ =  (S.2)  |µ1 |   |ν |−|µ if |ν | > |µ |, |ν | > |µ |. 1 1 2 2 2 1 |−1

34

Proof of Lemma S.2.5. By definition, {Λ1 , Λ2 } = {Λ∧ , Λ∨ }. W.o.l.g., we assume (Λ1 , Λ2 ) = (Λ∧ , Λ∨ ). Similarly, w.o.l.g., we assume (Θ1 , Θ2 ) = (Θ∧ , Θ∨ ). The first circumstance We begin with the first circumstances. Let V1∩2 be the selection matrix associated with ν1 ∩ ν2 . In this case, we can directly use the first argument of Lemma S.2.4 to get h  i EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ h  i ⊤ =EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V1∩2 V1∩2 V2 (U2⊤ XV2 )† U2⊤ h  i =EX,S,T tr U1 (U1⊤ XV1∩2 )†⊤ (U2⊤ XV1∩2 )† U2⊤ . Let ΠA = I − A⊤† A⊤ be the projection matrix associated to A. Then, we can break the above equation to get h  i EX,S,T tr U1 (U1⊤ XV1∩2 )†⊤ (U2⊤ XV1∩2 )† U2⊤ h  i =EX,S,T tr (U2⊤ XV1∩2 )† U2⊤ U1 (U1⊤ XV1∩2 )†⊤ h    i ⊤ ⊤ X ∨⊤ X ∨ V1∩2 )−1 V1∩2 X ∨⊤ U1 (U1⊤ XV1∩2 )†⊤ , =EX,S,T tr (U2⊤ XV1∩2 )† U2⊤ ΠX ∨ V1∩2 + X ∨ V1∩2 (V1∩2 (S.3) where we define X ∨ = Λ∨ X as a temporary notation. We scope in the first half of this quantity, where ⊤ ⊤ (U2⊤ XV1∩2 )† U2⊤ ΠX ∨ V1∩2 = (V1∩2 X ⊤ U2 U2⊤ XV1∩2 )−1 V1∩2 X ⊤ U2 U2⊤ ΠX ∨ V1∩2 .

Since the distribution of U2 U2⊤ is independent of the rest conditioned on we know |µ2 |, we have h i |µ2 | ⊤ ⊤ Λ2 ΠX ∨ V1∩2 ET (U2⊤ XV1∩2 )† U2⊤ ΠX ∨ V1∩2 =(V1∩2 X ⊤ U2 U2⊤ XV1∩2 )−1 V1∩2 X⊤ tr(Λ2 ) ⊤ ⊤ =(V1∩2 X ⊤ U2 U2⊤ XV1∩2 )−1 V1∩2 X ∨⊤ ΠX ∨ V1∩2 = 0, where we recall that we assumed Λ∨ = Λ2 . This brings (S.3) into h  i EX,S,T tr U1 (U1⊤ XV1∩2 )†⊤ (U2⊤ XV1∩2 )† U2⊤ h  i ⊤ ⊤ =EX,S,T tr (U2⊤ XV1∩2 )† U2⊤ X ∨ V1∩2 (V1∩2 X ∨⊤ X ∨ V1∩2 )−1 V1∩2 X ∨⊤ U1 (U1⊤ XV1∩2 )†⊤ h  i (i) ⊤ ⊤ =EX,S,T tr (U2⊤ XV1∩2 )† U2⊤ XV1∩2 (V1∩2 X ∨⊤ X ∨ V1∩2 )−1 V1∩2 X ⊤ U1 (U1⊤ XV1∩2 )†⊤ h  i ⊤ =EX,S,T tr (V1∩2 X ∨⊤ X ∨ V1∩2 )−1 . To show step (i), we observe the fact that U1⊤ X ∨ = U1⊤ X, and U2⊤ X ∨ = U2⊤ X.

(S.4)

This is because Λ∨ is a diagonal matrix with only 1s and 0s on its diagonal, while on its positions with 0s, U1 and U2 must have zero rows. Otherwise, the expectation of U1 and U2 can not have zero expectation on these positions. The last step is due to Lemma S.2.3, where we have h  i ⊤ EX,S,T tr (V1∩2 X ∨⊤ X ∨ V1∩2 )−1 =

35

|ν1 ∩ ν2 | . ∨ tr(Λ ) − |ν1 ∩ ν2 | − 1

To show the second condition, we notice that the operations between h  i h  i EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ = EX,S,T tr U1 (U1⊤ XV1∩2 )†⊤ (U2⊤ XV1∩2 )† U2⊤ are still lawful, where we continue to have h  i EX,S,T tr U1 (U1⊤ XV1∩2 )†⊤ (U2⊤ XV1∩2 )† U2⊤ h  i ⊤ ⊤ ⊤ =EX,S,T tr (V1∩2 X ⊤ U2 U2⊤ XV1∩2 )−1 V1∩2 X ⊤ U2 U2⊤ U1 U1⊤ XV1∩2 (V1∩2 X ⊤ U1 U1⊤ XV1∩2 )−1 h  i (ii) ⊤ ⊤ ⊤ = EX,S,T tr (V1∩2 X ⊤ U2 U2⊤ XV1∩2 )−1 V1∩2 X ⊤ U1 U1⊤ XV1∩2 (V1∩2 X ⊤ U1 U1⊤ XV1∩2 )−1 h  i |ν1 | ⊤ =EX,S,T tr (V1∩2 X ⊤ U2 U2⊤ XV1∩2 )−1 = , |µ2 | − |ν1 | − 1 where (ii) used the assumption that µ1 ⊂ µ2 , and the last step used Lemma S.2.3. We notice that ν1 ∩ ν2 = ν1 . The second circumstance Then we proceed to the second circumstance where |ν1 | > |µ1 |, |ν2 | > |µ2 |. This is an analog to the first circumstance, yet we operate on the dual dimension. Let U1∩2 be the selection matrix associated with µ1 ∩ µ2 . We use the first argument of Lemma S.2.4 to get h  i EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ h  i ⊤ =EX,S,T tr V2 (U2⊤ XV2 )† U2⊤ U1∩2 U1∩2 U1 (U1⊤ XV1 )†⊤ V1⊤ h  i ⊤ ⊤ =EX,S,T tr V2 (U1∩2 XV2 )† (U1∩2 XV1 )†⊤ V1⊤ h  i ⊤ ⊤ =EX,S,T tr (U1∩2 XV1 )†⊤ V1⊤ V2 (U1∩2 XV2 )† h    i ⊤ ⊤ ⊤ ⊤ XV1 )†⊤ V1⊤ ΠX ∨⊤ U1∩2 + X ∨⊤ U1∩2 (U1∩2 X ∨ X ∨⊤ U1∩2 )−1 U1∩2 X ∨ V2 (U1∩2 XV2 )† , =EX,S,T tr (U1∩2 where X ∨ = XΘ∨ . Following the same reasoning in the first circumstance, we get h  i h  i ⊤ EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ =EX,S,T tr (U1∩2 X ∨ X ∨⊤ U1∩2 )−1 =

|µ1 ∩ µ2 | . ∨ tr(Θ ) − |µ1 ∩ µ2 | − 1

To show the second condition, we notice that the operations between h  i h  i ⊤ ⊤ EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ = EX,S,T tr (U1∩2 XV1 )†⊤ V1⊤ V2 (U1∩2 XV2 )† are still lawful, where we continue to have h  i ⊤ ⊤ EX,S,T tr (U1∩2 XV1 )†⊤ V1⊤ V2 (U1∩2 XV2 )† h  i ⊤ ⊤ ⊤ =EX,S,T tr (U1∩2 XV1 V1⊤ X ⊤ U1∩2 )−1 U1∩2 XV1 V1⊤ V2 V2⊤ X ⊤ U1∩2 (U1∩2 XV2 V2⊤ X ⊤ U1∩2 )−1 h  i ⊤ ⊤ ⊤ =EX,S,T tr (U1∩2 XV1 V1⊤ X ⊤ U1∩2 )−1 U1∩2 XV1 V1⊤ X ⊤ U1∩2 (U1∩2 XV2 V2⊤ X ⊤ U1∩2 )−1 h  i |µ1 | ⊤ . =EX,S,T tr (U1∩2 XV2 V2⊤ X ⊤ U1∩2 )−1 = |ν2 | − |µ1 | − 1 □

36

Lemma S.2.6 (Bias) For µ1 , µ2 , ν1 , ν2 sampled independently, let U1 , U2 , V1 , V2 be their selection matrix. Moreover, let U1c be the selection matrix for µc1 , and so on for U2c , V1c , V2c . hSupposei Λ1 , Λ2 and |µ1 | Θ1 , Θ2 are diagonal matrices with only 1s and 0s on their diagonal. Then let ET U1 U1⊤ = tr(Λ Λ1 , 1) i i h h h i |µ2 | |ν1 | |ν2 | ⊤ ⊤ ⊤ ET U2 U2 = tr(Λ2 ) Λ2 , ES V1 V1 = tr(Θ1 ) Θ1 , ES V2 V2 = tr(Θ2 ) Θ2 . We want to calculate h  i (%) = EX,S,T tr X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ X . Assume we always have Λ1 Λ2 ∈ {Λ1 , Λ2 } and Θ1 Θ ∈ {Θ1 , Θ2 }. Define Λ∧ = Λ1 Λ2 , Λ∨ = Λ1 + Λ2 − Λ∧ , and Θ∧ = Θ1 Θ2 , Θ∨ = Θ1 + Θ2 − Θ∧ . Then there holds  |ν1 ∩ν2 ||ν1c ∩ν2c |     tr(Λ∨ )−|ν1 ∩ν2 |−1 + |ν1 ∩ ν2 | (%) =   |µ1 /µ2 ||µ2 /µ1 | d−|µ1 ∩µ2 |−1   tr(Θ ∨ )−|µ ∩µ | + |µ1 ∩ µ2 | tr(Θ∨ )−|µ ∩µ |−1 1 2 1 2

if |ν1 | < |µ1 |, |ν2 | < |µ2 |, if |ν1 | > |µ1 |, |ν2 | > |µ2 |.

(S.5)

Moreover, even if µ1 and µ2 are not independent, and ν1 and ν2 are not independent, yet we know µ1 ⊂ µ2 , ν1 ⊂ ν2 , the quantity becomes   |ν1 |(d−|ν2 |)    |µ2 |−|ν1 |−1 + |ν1 | if |ν1 | < |µ1 |, |ν2 | < |µ2 |, (%) =  (S.6)    |µ1 |(d−|ν2 |) + |µ1 | if |ν1 | > |µ1 |, |ν2 | > |µ2 |. |ν2 |−|µ1 |−1

Proof of Lemma S.2.6. By definition, {Λ1 , Λ2 } = {Λ∧ , Λ∨ }. W.o.l.g., we assume (Λ1 , Λ2 ) = (Λ∧ , Λ∨ ). Similarly, w.o.l.g., we assume (Θ1 , Θ2 ) = (Θ∧ , Θ∨ ). The first circumstance We first write i h  (%) = EX,S,T tr V1c V1c⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ XV2c V2c⊤ | {z } (∗)

h  i + EX,S,T tr V1 V1⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ XV2c V2c⊤ | {z } (∗∗)

h  i + EX,S,T tr V1c V1c⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ XV2 V2⊤ | {z } (∗∗∗)

h  i + EX,S,T tr V1 V1⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ XV2 V2⊤ . | {z } (∗∗∗∗)

We notice that h  i (∗∗) = EX,S,T tr V1 V1⊤ V2 (U2⊤ XV2 )† U2⊤ XV2c V2c⊤ = 0, where we used the fact that XV2c V2c⊤ is independent of the rest part and has zero mean. Similarly, (∗ ∗ ∗) = 0. It is also straightforward to see (∗ ∗ ∗∗) = |ν1 ∩ ν2 |. It remains to evaluate (∗). We continue

37

to use the notation X ∨ = Λ∨ X and X ∧ = Λ∧ X. We first notice that h i ET U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ h i =ET (ΠX ∧ V1 + X ∧ V1 (V1⊤ X ∧⊤ X ∧ V1 )−1 V1⊤ X ∧⊤ )U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ . The first term is a zero matrix since h i h i ET ΠX ∧ V1 U1 (U1⊤ XV1 )†⊤ =ET ΠX ∧ V1 U1 U1⊤ XV1 (V1⊤ X ⊤ U1 U1⊤ XV1 )−1 =ΠX ∧ V1 X ∧ V1 (V1⊤ X ⊤ U1 U1⊤ XV1 )−1 = 0, since the distribution of U1 U1⊤ is independent of the rest conditioned on we know |µ1 | and |µ1 ∩ µ2 |. The second term can be computed as h i ET (X ∧ V1 (V1⊤ X ∧⊤ X ∧ V1 )−1 V1⊤ X ∧⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ h i =ET (X ∧ V1 (V1⊤ X ∧⊤ X ∧ V1 )−1 V1⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ h i =ET (X ∧ V1 (V1⊤ X ∧⊤ X ∧ V1 )−1 V1⊤ V2 (U2⊤ XV2 )† U2⊤ , where we used (S.4) in the first equality. Via exactly the same strategy, the above quantity reduces to h i ET (X ∧ V1 (V1⊤ X ∧⊤ X ∧ V1 )−1 V1⊤ V2 (U2⊤ XV2 )† U2⊤ h i h i =ET (X ∧ V1 (V1⊤ X ∧⊤ X ∧ V1 )−1 V1⊤ V2 (V2⊤ X ∨⊤ X ∨ V2 )−1 V2⊤ X ∨⊤ = ET (V1⊤ X ∧⊤ )† V1⊤ V2 (X ∨ V2 )† . Bringing this into (∗), we get h  i (∗) = EX,S,T tr V1c V1c⊤ X ⊤ (V1⊤ X ∧⊤ )† V1⊤ V2 (X ∨ V2 )† XV2c V2c⊤ Let V1∩2 , V1∪2 , V1/2 , and V2/1 be the selection matrix associated to ν1 ∩ ν2 , ν1 ∪ ν2 , ν1 /ν2 , and ν2 /ν1 , respectively. Since order-invariant permutation is allowed, we arrive h i c c⊤ c c⊤ (∗) =EX,S ,T tr(V1∪2 V1∪2 X ⊤ (V1⊤ X ∧⊤ )† V1⊤ V2 (X ∨ V2 )† XV1∪2 V1∪2 ) h i c c⊤ ⊤ c c⊤ =EX,S ,T tr(V1∪2 V1∪2 X ⊤ (V1⊤ X ∧⊤ )† V1⊤ V1∩2 V1∩2 V2 (X ∨ V2 )† XV1∪2 V1∪2 ) h i c c⊤ ⊤ c c⊤ =EX,S ,T tr(V1∪2 V1∪2 X ⊤ (V1∩2 X ∧⊤ )† (X ∨ V1∩2 )† XV1∪2 V1∪2 ) . where in the last step we used the first argument of Lemma S.2.4. Note that this is applicable since c V c⊤ is independent of both XV XV1∪2 1/2 and XV2/1 . We first note that 1∪2 h i c c⊤ c c⊤ ⊤ c V c⊤ XV = |ν1c ∩ ν2c | · In . EXV1∪2 1∪2 V1∪2 V1∪2 V1∪2 X 1∪2 Also, there holds ⊤ ⊤ ⊤ ⊤ tr((V1∩2 X ∧⊤ )† (X ∨ V1∩2 )† ) = tr((V1∩2 X ∧⊤ X ∧ V1∩2 )−1 V1∩2 X ∧⊤ X ∨ V1∩2 (V1∩2 X ∨⊤ X ∨ V1∩2 )−1 ) ⊤ ⊤ ⊤ = tr((V1∩2 X ∧⊤ X ∧ V1∩2 )−1 V1∩2 X ∧⊤ X ∧ V1∩2 (V1∩2 X ∨⊤ X ∨ V1∩2 )−1 ) |ν1 ∩ ν2 | ⊤ = tr((V1∩2 X ∨⊤ X ∨ V1∩2 )−1 ) = . ∨ tr(Λ ) − |ν1 ∩ ν2 | − 1

38

They together yield |ν1c ∩ ν2c | · |ν1 ∩ ν2 | (∗) = . tr(Λ∨ ) − |ν1 ∩ ν2 | − 1 Thus, (%) = (∗) + (∗ ∗ ∗∗) =

|ν1c ∩ ν2c | · |ν1 ∩ ν2 | + |ν1 ∩ ν2 |. tr(Λ∨ ) − |ν1 ∩ ν2 | − 1

As for the second conclusion, it is still straightforward that (∗∗) = (∗ ∗ ∗) = 0. Also, (∗ ∗ ∗∗) = |γ1 |. For (∗ ∗ ∗∗), we have h  i (∗) =EX,S,T tr V1c V1c⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ XV2c V2c⊤ h  i =EX,S,T tr V2c V2c⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ XV2c V2c⊤ h  i =EX,S,T tr V2c V2c⊤ X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V1 (U2⊤ XV1 )† U2⊤ XV2c V2c⊤ , where the last follows from Lemma S.2.4. Lemma S.2.3 yields h i EXV2c XV2c V2c⊤ V2c V2c⊤ X ⊤ = (d − |ν2 |)In . Then, h  i (∗) =(d − |ν2 |)EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V1 (U2⊤ XV1 )† U2⊤ h  i =(d − |ν2 |)EX,S,T tr U1 (U1⊤ XV1 )†⊤ V1⊤ V1 (U2⊤ XV1 )† U1⊤ h  i =(d − |ν2 |)EX,S,T tr U1 U1⊤ XV1 (V1⊤ X ⊤ U1 U1⊤ XV1 )−1 V1⊤ V1 (V1⊤ X ⊤ U2 U2⊤ XV1 )−1 V1⊤ X ⊤ h  i |ν1 |(d − |ν2 |) =(d − |ν2 |)EX,S,T tr (V1⊤ X ⊤ U2 U2⊤ XV1 )−1 = . |µ2 | − |ν1 | − 1 The second circumstance Then we proceed to the second circumstance where |ν1 | > |µ1 |, |ν2 | > |µ2 |. This is an analog to the first circumstance, yet we operate on the dual dimension. We switch the notation to use the notation X ∨ = XΘ∨ and X ∧ = XΘ∧ . Let U1∩2 , U1∪2 , U1/2 , and U2/1 be the selection matrix associated to µ1 ∩ µ2 , µ1 ∪ µ2 , µ1 /µ2 , and µ2 /µ1 , respectively. We first write h  i (%) = EX,S,T tr X ⊤ U2c U2c⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ U1c U1c⊤ X | {z } (∗)

h  i + EX,S,T tr X ⊤ U2c U2c⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ U1 U1⊤ X | {z } (∗∗)

h  i + EX,S,T tr X ⊤ U2 U2⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ U1c U1c⊤ X | {z } (∗∗∗)

h  i + EX,S,T tr X ⊤ U2 U2⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ U1 U1⊤ X | {z } (∗∗∗∗)

h i (i) We begin with (∗). Due to a similar reasoning in the first circumstance, we know that ES (U1⊤ XV1 )†⊤ V1⊤ ΠX ∨⊤ U1 =

39

0. Thus, we can do the following transformation h  i (∗) =EX,S,T tr X ⊤ U2c U2c⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ X ∨⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ X ∨ V2 (U2⊤ XV2 )† U2⊤ U1c U1c⊤ X h  i =EX,S,T tr X ⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ X ∨ V2 (U2⊤ XV2 )† U2⊤ U1c U1c⊤ X . Via exactly the same strategy, the above quantity reduces to h  i (∗) =EX,S,T tr X ⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ X ∨ X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X h  i =EX,S,T tr X ⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ X ∧ X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X . (S.7) We also note that i h i h EXΘ∨c U2⊤ XX ⊤ U2c U2c⊤ =U2⊤ X ∨ X ∨⊤ U2c U2c⊤ + EXΘ∨c U2⊤ XΘ∨c Θ∨c X ⊤ U2c U2c⊤ h i =U2⊤ X ∨ X ∨⊤ U2c U2c⊤ + tr(Θ∨c )EXΘ∨c U2⊤ U2c U2c⊤ = U2⊤ X ∨ X ∨⊤ U2c U2c⊤ . This means we can turn the Xs on the side into h  i (∗) = EX,S,T tr X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ X ∧ X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X ∨ . We further decompose In = U2 U2⊤ + U2c U2c⊤ and have h  i (∗) =EX,S,T tr X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ U2 U2⊤ X ∧ X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X ∨ h  i +EX,S,T tr X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ U2c U2c⊤ X ∧ X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X ∨ . The first term becomes h  i EX,S,T tr X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ U2 U2⊤ X ∧ X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X ∨ h  i =EX,S,T tr X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ U2 U2⊤ U1c U1c⊤ X ∨ = 0. (S.8) ⊤ X ∨ = X ∨ , U ⊤ X ∨ = X ∨ , U ⊤ X ∨ = X ∨ , U c⊤ X ∨ = For the second term, we write U1∩2 1∩2 2/1 2/1 1/2 1/2 1∪2 X1∨c ∩2c . W.o.l.g., we assume (µ1 ∩ µ2 , µ1 /µ2 , µ2 /µ1 , µc1 ∩ µc2 ) = (1, . . . , d), which means X = ∨⊤ , X ∨⊤ , X ∨⊤ , X ∨⊤ )⊤ . Then we have (X1∩2 1c ∩2c 1/2 2/1



X ∨ X ∨⊤ 1∩2 U1⊤ X ∨ X ∨⊤ U1 =  1∩2 ∨ X ∨⊤ X1/2 1∩2

 ∨ X ∨⊤  X1∩2  . 1/2   ∨ X ∨⊤  X1/2 1/2

We can compute X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ U2c U2c⊤ X ∨ ∨⊤ ∨ =X1/2 (U1⊤ X ∨ X ∨⊤ U1 )−1 X1/2 ∨⊤ ∨ ∨⊤ ∨ ∨⊤ ∨ ∨⊤ −1 ∨ ∨⊤ −1 ∨ ∨⊤ ∨ ∨⊤ −1 ∨ ∨⊤ X =X1/2 (X1/2 X1/2 − X1/2 X1∩2 (X1∩2 X1∩2 ) X1∩2 X1/2 ) X1/2 = X1/2 (X1/2 ΠX1∩2 1/2 ) X1/2 .

(S.9) ∨ has i.i.d. standard Gaussian entries, decomposing it into X ∨⊤ = Π ∨⊤ X ∨⊤ + Π ∨c X ∨⊤ Since X1/2 X1∩2 1/2 X1∩2 1/2 1/2

40

yields the following, where the cross terms disappear. i h ∨⊤ −1 ∨ ∨⊤ ∨ ∨⊤ X ∨⊤ X EX1/2 1/2 ) X1/2 1/2 (X1/2 ΠX1∩2 i i h h ∨⊤ −1 ∨ ∨⊤ −1 ∨ ∨⊤ ∨ ∨⊤ ∨ ∨ ∨⊤ ΠX ∨c X ∨⊤ X ∨c ∨⊤ X ∨⊤ ΠX ∨ X ) X Π . ) X Π + E (X Π (X Π =EX1/2 X X X X X 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1∩2 1∩2 1∩2 1∩2 1/2 1∩2 1∩2 The first term has   ∨⊤ ΠΠ ∨⊤ EX1/2 = X ∨⊤ X1/2 1∩2

|µ1 /µ2 | Π ∨ , tr(Θ∨ ) − |µ1 ∩ µ2 | X1∩2

∨⊤ ∨ ∨ X which is by the fact that ΠX1∩2 1/2 is distributed as N(0, ΠX1∩2 ). The second term has

i h ∨⊤ −1 ∨ ∨⊤ ∨ ∨⊤ ΠX ∨c X ∨⊤ X ∨c ) X Π EX1/2 (X Π X X 1/2 1/2 1/2 1/2 1∩2 1∩2 1∩2 " # I |µ /µ | 1 2 ∨⊤ ∨ ∨⊤ ΠX ∨c X =EX1/2 X Π ∨c 1/2 1∩2 tr(Θ∨ ) − |µ1 ∩ µ2 | − |µ1 /µ2 | − 1 1/2 X1∩2 |µ1 /µ2 | ∨c =ΠX1∩2 ∨ tr(Θ ) − |µ1 ∩ µ2 | − |µ1 /µ2 | − 1 due to Lemma S.2.3. This means the second term of (∗) is (·) := =

|µ1 /µ2 | |µ1 /µ2 | ∨ + ΠX ∨c ΠX1∩2 ∨ ∨ 1∩2 tr(Θ ) − |µ ∩ µ | − |µ /µ | − 1 tr(Θ ) − |µ1 ∩ µ2 | 1 2 1 2

|µ1 /µ2 | |µ1 /µ2 | ∨c Itr(Θ∨ ) + ΠX1∩2 . tr(Θ∨ ) − |µ1 ∩ µ2 | (tr(Θ∨ ) − |µ1 ∩ µ2 | − |µ1 /µ2 | − 1)(tr(Θ∨ ) − |µ1 ∩ µ2 |) (S.10)

Bringing this into (∗), we get h  i (∗) = EX,S,T tr (·)X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X ∨ . ∨c = 0, and thus Note that U1c U1c⊤ X ∨ ΠX1∩2

(∗) =

h  i |µ1 /µ2 | ∧⊤ ⊤ ∧ ∧⊤ −1 ⊤ c c⊤ ∨ E tr X U (U X X U ) U U U X . X,S,T 2 2 2 2 1 1 tr(Θ∨ ) − |µ1 ∩ µ2 |

Note that we can change the last X ∨ into X ∧ since trace operator can permutate and X ∨ X ∧⊤ = X ∧ X ∧⊤ . Following the same notations, we have X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X ∧ ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ ∧ ∧⊤ ∧ ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ −1 ∧ = − X1∩2 (X1∩2 X1∩2 ) X1∩2 X2/1 (X2/1 X2/1 − X2/1 X1∩2 (X1∩2 X1∩2 ) X1∩2 X2/1 ) X2/1 ∧⊤ ∧ ∧⊤ ∧ ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ −1 ∧ (X2/1 X2/1 − X2/1 X1∩2 (X1∩2 X1∩2 ) X1∩2 X2/1 ) X2/1 + X2/1 ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ X ∧⊤ X =ΠX1∩2 2/1 (X2/1 ΠX1∩2 2/1 ) X2/1 ,

whose trace is h  i EX,S,T tr X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1c U1c⊤ X ∧ h  i ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ X ∧⊤ X ∧⊤ = |µ2 /µ1 |. =EX,S,T tr ΠX1∩2 2/1 (X2/1 ΠX1∩2 2/1 ) X2/1 ΠX1∩2

41

Thus, we reach (∗) =

|µ1 /µ2 ||µ2 /µ1 | . tr(Θ∨ ) − |µ1 ∩ µ2 |

(ii) For (∗∗), notice that the derivations till (S.8) and (S.10) are still valid, except we need to replace U1c U1c⊤ by U1 U1⊤ , namely h  i (∗∗) =EX,S,T tr X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ U2 U2⊤ U1 U1⊤ X ∨ h  i |µ1 /µ2 | ∧⊤ ⊤ ∧ ∧⊤ −1 ⊤ ⊤ ∨ + E tr X U (U X X U ) U U U X . X,S,T 2 2 1 2 2 1 tr(Θ∨ ) − |µ1 ∩ µ2 | − |µ1 /µ2 | − 1 Under the same notation, we can compute ∨⊤ ∨ ∨⊤ −1 ∨ ∨⊤ X ∨c⊤ X ∨⊤ U2c U2c⊤ U1 (U1⊤ X ∨ X ∨⊤ U1 )−1 U1⊤ U2 U2⊤ X ∨ = −X1/2 (X1/2 ΠX1∩2 1/2 ) X1/2 ΠX1∩2

and X ∧⊤ U2 (U2⊤ X ∧ X ∧⊤ U2 )−1 U2⊤ U1 U1⊤ X ∨ ∧⊤ ∧ ∧⊤ −1 ∧ =X1∩2 (X1∩2 X1∩2 ) X1∩2 ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ X +X1∩2 (X1∩2 X1∩2 ) X1∩2 X2/1 (X2/1 ΠX1∩2 2/1 ) X2/1 X1∩2 (X1∩2 X1∩2 ) X1∩2 ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ X −X2/1 (X2/1 ΠX1∩2 2/1 ) X2/1 X1∩2 (X1∩2 X1∩2 ) X1∩2 ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ ∧ ∧⊤ −1 ∧ ∧⊤ X ∧⊤ X ∧c⊤ . =X1∩2 (X1∩2 X1∩2 ) X1∩2 + ΠX1∩2 2/1 (X2/1 ΠX1∩2 2/1 ) X2/1 ΠX1∩2

Bringing in the same calculations as in (S.9) yields (∗∗) = 0. (iii) The calculation of (∗ ∗ ∗) follows symmetrically to (∗). (iv) We directly apply Lemma S.2.4 to have h  i EX,S,T tr X ⊤ U2 U2⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V2 (U2⊤ XV2 )† U2⊤ U1 U1⊤ X h  i ⊤ ⊤ ⊤ =EX,S,T tr X ⊤ U1∩2 (U1∩2 XV1 )†⊤ V1⊤ V2 (U1∩2 XV2 )† U1∩2 X , which is followed by the transformations similar to (S.7) to have h  i ⊤ ⊤ (∗ ∗ ∗∗) =EX,S,T tr X ⊤ U1∩2 (U1∩2 X ∨ X ∨⊤ U1∩2 )−1 U1∩2 X h  i ⊤ ⊤ =EX,S,T tr X ∨⊤ U1∩2 (U1∩2 X ∨ X ∨⊤ U1∩2 )−1 U1∩2 X∨ h  i ⊤ ⊤ +EX,S,T tr X ⊤ U1∩2 (U1∩2 X ∨ X ∨⊤ U1∩2 )−1 U1∩2 X(Id − Θ∨ ) . The first term is |µ1 ∩ µ2 |, while the second is (d − tr(Θ∨ ))|µ1 ∩ µ2 | . tr(Θ∨ ) − |µ1 ∩ µ2 | − 1 They together yield (∗ ∗ ∗∗) = |µ1 ∩ µ2 |

d − |µ1 ∩ µ2 | − 1 . tr(Θ∨ ) − |µ1 ∩ µ2 | − 1

42

Thus, (%) = (∗) + (∗ ∗ ∗∗) =

d − |µ1 ∩ µ2 | − 1 |µ1 /µ2 ||µ2 /µ1 | + |µ1 ∩ µ2 | . ∨ tr(Θ ) − |µ1 ∩ µ2 | tr(Θ∨ ) − |µ1 ∩ µ2 | − 1

For the second conclusion, we notice that (∗) = (∗∗) = (∗ ∗ ∗) = 0 is obvious. For (∗ ∗ ∗∗), we have h  i (∗ ∗ ∗∗) =EX,S,T tr X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V1 (U2⊤ XV2 )† U2⊤ U1 U1⊤ X h  i =EX,S,T tr X ⊤ U1 (U1⊤ XV1 )†⊤ V1⊤ V1 (U1⊤ XV2 )† U1⊤ X by Lemma S.2.4, from where we proceed to have h  i (∗ ∗ ∗∗) =EX,S,T tr X ⊤ U1 (U1⊤ XV1 V1⊤ X ⊤ U1 )−1 U1⊤ XV1 V1⊤ X ⊤ U1 (U1⊤ XV2 V2⊤ X ⊤ U1 )−1 U1⊤ X h  i =EX,S,T tr X ⊤ U1 (U1⊤ XV2 V2⊤ X ⊤ U1 )−1 U1⊤ X =|µ1 | +

|µ1 |(d − |ν2 |) . |ν2 | − |µ1 | − 1 □

S.2.3

Proofs for Appendix S.2.1

Proof of Proposition S.2.1. Recall that 

 tr g(x) − g′ (x) ⊤ g(x) − g′ (x) B    1 X  ′ † ′⊤ ⊤ ⊤ † ⊤ ′ ′⊤ ′ † ′⊤ = 2 tr Vb1 (Ub⊤1 XVb1 )† Ub⊤1 − Vb′1 (Ub′⊤ XV ) U V (U XV ) U − V (U XV ) U . b b 2 2 b1 b1 b2 b2 b2 b2 b2 b2 1 B b ,b =1 | {z } 1 2 Vb1 ,b2

There are B(B − 1) cross terms which have b1 , b2 , and B squared terms with b1 = b2 . Underparameterized case For instance removal, the squared terms become     ⊤ † ⊤ ′⊤ † ′⊤ ⊤ ⊤ † ⊤ ′⊤ † ′⊤ Vb,b = tr Vb (Ub XVb ) Ub − Vb (Ub XVb ) Ub Vb (Ub XVb ) Ub − Vb (Ub XVb ) Ub , where Ub and Ub′ are independent. We apply (S.2) in Lemma S.2.5 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ Vb (Ub⊤ XVb )† Ub⊤   ⊤   =EX,S,T tr Vb (Ub′⊤ XVb )† Ub′⊤ Vb (Ub′⊤ XVb )† Ub′⊤ =

|νb | . |µb | − |νb | − 1

For the other two cross terms, we treat them as only have νb columns, and apply (S.1) in Lemma S.2.5 to have   ⊤   |νb | EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ Vb (Ub′⊤ XVb )† Ub′⊤ = . n − |νb | − 1 Combined together, this is   EX,S,T Vb,b = 2

|νb | |νb | −2 . |µb | − |νb | − 1 n − |νb | − 1

43

Taking the limit, this is   γq γq 2γq(1 − p) lim EX,S,T Vb,b = 2 −2 = . n,d→∞ p − γq 1 − γq (p − γq)(1 − γq) For feature removal, the squared term becomes  ⊤   Vb,b = tr Vb (Ub⊤ XVb )† Ub⊤ − Vb′ (Ub⊤ XVb′ )† Ub⊤ Vb (Ub⊤ XVb )† Ub⊤ − Vb′ (Ub⊤ XVb′ )† Ub⊤ , where Vb and Vb′ are independent. We apply (S.2) in Lemma S.2.5 to have      ⊤ † ⊤ ⊤ ⊤ † ⊤ EX,S,T tr Vb (Ub XVb ) Ub Vb (Ub XVb ) Ub   ⊤   =EX,S,T tr Vb′ (Ub⊤ XVb′ )† Ub⊤ Vb′ (Ub⊤ XVb′ )† Ub⊤ =

|νb | . |µb | − |νb | − 1

For the other two cross terms, we treat them as only have µb rows, and apply (S.1) in Lemma S.2.5 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ Vb′ (Ub⊤ XVb′ )† Ub⊤ =

|νb ∩ νb′ | . np − |νb ∩ νb′ | − 1

Combined together, this is   EX,S,T Vb,b = 2

|νb ∩ νb′ | |νb | . −2 |µb | − |νb | − 1 np − |νb ∩ νb′ | − 1

Taking the limit, this is   γq γq2 2γpq(1 − q) lim EX,S,T Vb,b = 2 −2 = . 2 n,d→∞ p − γq p − γq (p − γq)(p − γq2 ) For the cross terms (b1 , b2 ) for instance removal, we have " #   |ν1 ∩ ν2 | |ν1 ∩ ν2 | lim n · EX,S,T Vb1 ,b2 = lim n · EX,S,T − n,d→∞ n,d→∞ n − 1 − |ν1 ∩ ν2 | − 1 n − |ν1 ∩ ν2 | − 1   2 s s2   γq2 d d  = = lim n ·  − 2 2  (1 − γq2 )2 . n,d→∞ n − sd − 2 n − sd − 1 As for the cross terms (b1 , b2 ) for feature removal, we have " # |ν1′ ∩ ν2′ |   |ν1 ∩ ν2 | − lim n · EX,S,T Vb1 ,b2 = lim n · EX,S,T n,d→∞ n,d→∞ n − |ν1′ ∩ ν2′ | − 1 n − |ν1 ∩ ν2 | − 1   s2 s2   q2 d−1 d  = = lim n ·  − . 2 2 s n,d→∞ (1 − γq2 )2 n − d−1 − 1 n − sd − 1 Overparameterized case For instance removal, the squared term becomes  ⊤   Vb,b = tr Vb (Ub⊤ XVb )† Ub⊤ − Vb (Ub′⊤ XVb )† Ub′⊤ Vb (Ub⊤ XVb )† Ub⊤ − Vb (Ub′⊤ XVb )† Ub′⊤ ,

44

where Ub and Ub′ are independent. We apply (S.2) in Lemma S.2.5 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ Vb (Ub⊤ XVb )† Ub⊤   ⊤   =EX,S,T tr Vb (Ub′⊤ XVb )† Ub′⊤ Vb (Ub′⊤ XVb )† Ub′⊤ =

|µb | . |νb | − |µb | − 1

For the other two cross terms, we treat them as only have νb columns, and apply (S.1) in Lemma S.2.5 to have      ⊤ † ⊤ ⊤ ′ ⊤ ′ † ⊤ EX,S,T tr Vb (Ub XVb ) Ub Vb (Ub XVb ) Ub =

|µb ∩ µ′b | . qd − |µb ∩ µ′b | − 1

Combined together, this is |µb ∩ µ′b | |µb | EX,S,T Vb,b = 2 . −2 |νb | − |µb | − 1 qd − |µb ∩ µ′b | − 1 



Taking the limit, this is   lim EX,S,T Vb,b = 2

n,d→∞

p p2 2γq(p − p2 ) −2 = . γq − p γq − p2 (γq − p)(γq − p2 )

For feature removal, the squared terms become  ⊤   Vb,b = tr Vb (Ub⊤ XVb )† Ub⊤ − Vb′ (Ub⊤ XVb′ )† Ub⊤ Vb (Ub⊤ XVb )† Ub⊤ − Vb′ (Ub⊤ XVb′ )† Ub⊤ , where Vb and Vb′ are independent. We apply (S.2) in Lemma S.2.5 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ Vb (Ub⊤ XVb )† Ub⊤   ⊤   =EX,S,T tr Vb′ (Ub⊤ XVb′ )† Ub⊤ Vb′ (Ub⊤ XVb′ )† Ub⊤ =

|µb | . |νb | − |µb | − 1

For the other two cross terms, we treat them as only have µb rows, and apply (S.1) in Lemma S.2.5 to have   ⊤   |µb | EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ Vb′ (Ub⊤ XVb′ )† Ub⊤ = . d − |µb | − 1 Combined together, this is   EX,S,T Vb,b = 2

|µb | |µb | −2 . |νb | − |µb | − 1 d − |µb | − 1

Taking the limit, this is   lim EX,S,T Vb,b = 2

n,d→∞

p p 2γp(1 − q) −2 = . γq − p γ − p (γq − p)(γ − p)

45

For the cross terms (b1 , b2 ) with instance removal, we have " # |µ′1 ∩ µ′2 |   |µ1 ∩ µ2 | lim n · EX,S,T Vb1 ,b2 = lim n · EX,S,T − n,d→∞ n,d→∞ d − |µ′1 ∩ µ′2 | − 1 d − |µ1 ∩ µ2 | − 1   m2 m2   γp2 n−1 n  = = lim n ·  . − 2 2 m n,d→∞ (γ − p2 )2 d − n−1 − 1 d − mn − 1 For the cross terms (b1 , b2 ) with feature removal, we have " #   |µ1 ∩ µ2 | |µ1 ∩ µ2 | lim n · EX,S,T Vb1 ,b2 = lim n · EX,S,T − n,d→∞ n,d→∞ d − 1 − |µ1 ∩ µ2 | − 1 d − |µ1 ∩ µ2 | − 1   m2 m2   p2 n n  =  = lim n ·  − 2 2  (γ − p2 )2 . n,d→∞ d − 1 − mn − 1 d − mn − 1 □ Proof of Proposition S.2.2. Recall that

i h

Eβ∗ β ∗⊤ X ⊤ g(x) − g′ (x) ⊤ g(x) − g′ (x) Xβ ∗ B ⊤   1 X  ′ † ′⊤ ⊤ † ⊤ ′ ′⊤ ′ † ′⊤ = 2 tr Vb1 (Ub⊤1 XVb1 )† Ub⊤1 X − Vb′1 (Ub′⊤ XV ) U X V (U XV ) U X − V (U XV ) U X b2 b2 b1 b1 b2 b2 b2 b2 b2 b2 1 dB b ,b =1 | {z } 1 2 Bb1 ,b2

There are B(B − 1) terms which have b1 , b2 , and B terms with b1 = b2 . Underparameterized case For instance removal, the squared terms become  ⊤   Bb,b = tr Vb (Ub⊤ XVb )† Ub⊤ X − Vb (Ub′⊤ XVb )† Ub′⊤ X Vb (Ub⊤ XVb )† Ub⊤ X − Vb (Ub′⊤ XVb )† Ub′⊤ X , where Ub and Ub′ are independent. We apply (S.6) in Lemma S.2.6 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ X Vb (Ub⊤ XVb )† Ub⊤ X    ⊤  |νb |(d − |νb |) + |νb |. =EX,S,T tr Vb (Ub′⊤ XVb )† Ub′⊤ X Vb (Ub′⊤ XVb )† Ub′⊤ X = |µb | − |νb | − 1 For the other two cross terms, we treat them as only have νb columns, and apply (S.5) in Lemma S.2.6 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ X Vb (Ub′⊤ XVb )† Ub′⊤ X =

|νb ||νbc | n − |νb | − 1

+ |νb |.

Combined together, this is |νb ||νbc |   |νb |(d − |νb |) + |νb |) − 2( + |νb |). EX,S,T Bb,b = 2( |µb | − |νb | − 1 n − |νb | − 1 Taking the limit, this is   1 γq(1 − q) γq(1 − q) 2γq(1 − q)(1 − p) EX,S,T Bb,b = 2 −2 = . n,d→∞ d p − γq 1 − γq (p − γq)(1 − γq) lim

46

For feature removal, the squared term becomes  ⊤   Bb,b = tr Vb (Ub⊤ XVb )† Ub⊤ X − Vb′ (Ub⊤ XVb′ )† Ub⊤ X Vb (Ub⊤ XVb )† Ub⊤ X − Vb′ (Ub⊤ XVb′ )† Ub⊤ X , where Vb and Vb′ are independent. We apply (S.6) in Lemma S.2.6 to have   ⊤   ⊤ † ⊤ ⊤ † ⊤ EX,S,T tr Vb (Ub XVb ) Ub X Vb (Ub XVb ) Ub X   ⊤   |νb |(d − |νb |) =EX,S,T tr Vb′ (Ub⊤ XVb′ )† Ub⊤ X Vb′ (Ub⊤ XVb′ )† Ub⊤ X = + |νb |. |µb | − |νb | − 1 For the other two cross terms, we treat them as only have νb rows, and apply (S.5) in Lemma S.2.6 to have   ⊤   |νb ∩ νb′ ||νbc ∩ νb′c | + |νb ∩ νb′ |. EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ X Vb′ (Ub⊤ XVb′ )† Ub⊤ X = np − |νb ∩ νb′ | − 1 Combined together, this is EX,S,T



! ! |νb ∩ νb′ ||νbc ∩ νb′c | |νb |(d − |νb |) ′ + |νb | − 2 + |νb ∩ νb | . Bb,b = 2 |µb | − |νb | − 1 np − |νb ∩ νb′ | − 1 

Taking the limit, this is   1 γq(1 − q) γq2 (1 − q)2 2γpq(1 − q)2 2pq(1 − q) EX,S,T Bb,b = 2 −2 + 2q(1 − q) = + . 2 2 n,d→∞ d p − γq p − γq (p − γq)(p − γq ) p − γq2 lim

For the cross terms (b1 , b2 ) for instance removal, we have n

  · EX,S,T Bb1 ,b2   c c |νb1 ∩ νb2 ||νbc1 ∩ νbc2 |  |νb1 ∩ νb2 ||νb1 ∩ νb2 |  n = lim · EX,S,T  + |νb1 ∩ νb2 | − − |νb1 ∩ νb2 | ∨ ∨ n,d→∞ d tr(Λ ) − |νb1 ∩ νb2 | − 1 tr(Λ ) − |νb1 ∩ νb2 | − 1   2 2 2 2 s (d−s) s (d−s) s2 s2  γq2 (1 − q)2 n  d d d d  = = lim ·  + − − . 2 2 n,d→∞ d (1 − γq2 )2 n − 1 − sd − 1 d n − sd − 1 d lim

n,d→∞ d

As for the cross terms (b1 , b2 ) for feature removal, we have n

  · EX,S,T Bb1 ,b2  ′  |νb1 ∩ νb2 ||νbc1 ∩ νbc2 |   |νb1 ∩ νb′ 2 ||νb′c1 ∩ νb′c2 | n ′ ′ = lim · EX,S,T  + |νb1 ∩ νb2 | − − |νb1 ∩ νb2 | ′ ′ ∨ ∨ n,d→∞ d tr(Λ ) − |νb1 ∩ νb2 | − 1 tr(Λ ) − |νb1 ∩ νb2 | − 1  s2 (d−1−s)2  2 s2 (d−s) n  d−1 ( d−2 + 1) s2 s2   = lim ·  + − d d−1 −  2 2 s n,d→∞ d d−1 n− s −1 d n − d−1 −1 d lim

n,d→∞ d

=

γq4 (1 − q2 ) − 2q3 (1 − q) q2 + + γ−1 q2 . (1 − γq2 )2 1 − γq2

Overparameterized case For instance removal, the squared terms become  ⊤   Bb,b = tr Vb (Ub⊤ XVb )† Ub⊤ X − Vb (Ub′⊤ XVb )† Ub′⊤ X Vb (Ub⊤ XVb )† Ub⊤ X − Vb (Ub′⊤ XVb )† Ub′⊤ X ,

47

where Ub and Ub′ are independent. We apply (S.6) in Lemma S.2.6 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ X Vb (Ub⊤ XVb )† Ub⊤ X   ⊤   |µb |(d − |νb |) =EX,S,T tr Vb (Ub′⊤ XVb )† Ub′⊤ X Vb (Ub′⊤ XVb )† Ub′⊤ X = + |µb |. |νb | − |µb | − 1 For the other two cross terms, we treat them as only have νb columns, and apply (S.5) in Lemma S.2.6 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ X Vb (Ub′⊤ XVb )† Ub′⊤ X =

d − |µb ∩ µ′b | − 1 |µb /µ′b ||µ′b /µb | ′ + |µ ∩ µ | . b b qd − |µb ∩ µ′b | − 1 qd − |µb ∩ µ′b | − 1

Combined together, this is |µb /µ′b ||µ′b /µb | d − |µb ∩ µ′b | − 1   |µb |(d − |νb |) ′ EX,S,T Bb,b = 2( + |µb |) − 2( + |µ ∩ µ | ). b b |νb | − |µb | − 1 qd − |µb ∩ µ′b | − 1 qd − |µb ∩ µ′b | − 1 Taking the limit, this is 2   p(1 − q) γ−1 p2 (1 − p)2 1 −1 2 γ − p EX,S,T Bb,b =2 + 2γ−1 p − 2 p − 2γ . n,d→∞ d γq − p γq − p2 γq − p2

lim

For feature removal, the squared term becomes  ⊤   Bb,b = tr Vb (Ub⊤ XVb )† Ub⊤ X − Vb′ (Ub⊤ XVb′ )† Ub⊤ X Vb (Ub⊤ XVb )† Ub⊤ X − Vb′ (Ub⊤ XVb′ )† Ub⊤ X , where Vb and Vb′ are independent. We apply (S.6) in Lemma S.2.6 to have   ⊤   ⊤ † ⊤ ⊤ † ⊤ EX,S,T tr Vb (Ub XVb ) Ub X Vb (Ub XVb ) Ub X   ⊤   |µb |(d − |νb |) =EX,S,T tr Vb′ (Ub⊤ XVb′ )† Ub⊤ X Vb′ (Ub⊤ XVb′ )† Ub⊤ X = + |µb |. |νb | − |µb | − 1 For the other two cross terms, we treat them as only have µb rows, and apply (S.5) in Lemma S.2.6 to have   ⊤   EX,S,T tr Vb (Ub⊤ XVb )† Ub⊤ X Vb′ (Ub⊤ XVb′ )† Ub⊤ X = |µb |. Combined together, this is EX,S,T



! |µb |(d − |νb |) Bb,b = 2 + |µb | − 2|µb |. |νb | − |µb | − 1 

Taking the limit, this is lim

1

n,d→∞ d

  2p(1 − q) EX,S,T Bb,b = . γq − p

48

As for the cross terms (b1 , b2 ) for instance removal, we have n

  · EX,S,T Bb1 ,b2 " # |µ′1 /µ′2 ||µ′2 /µ′1 | d − |µ′1 ∩ µ′2 | − 1 n ′ ′ = lim · EX,S,T + |µ1 ∩ µ2 | n,d→∞ d tr(Θ∨ ) − |µ′1 ∩ µ′2 | tr(Θ∨ ) − |µ′1 ∩ µ′2 | − 1 " # |µ1 /µ2 ||µ2 /µ1 | n d − |µ1 ∩ µ2 | − 1 − lim · EX,S,T + |µ1 ∩ µ2 | n,d→∞ d tr(Θ∨ ) − |µ1 ∩ µ2 | tr(Θ∨ ) − |µ1 ∩ µ2 | − 1  m2 (n−1−m)2  2 m2 (n−m) n  n−1 n−2 m2 m2    = lim ·  + − n n−12 − 2 m n,d→∞ d n−1 n  d − n−1 d − mn lim

n,d→∞ d

=

γ−1 p4 (1 − p2 ) − 2p3 (1 − p) + γ−1 p2 . (γ − p2 )2

As for the cross terms (b1 , b2 ) for feature removal, we have n

  · EX,S,T Bb1 ,b2 " # n |µ1 /µ2 ||µ2 /µ1 | d − |µ1 ∩ µ2 | − 1 = lim · EX,S,T + |µ1 ∩ µ2 | n,d→∞ d d − 1 − |µ1 ∩ µ2 | d − 1 − |µ1 ∩ µ2 | − 1 " # n |µ1 /µ2 ||µ2 /µ1 | d − |µ1 ∩ µ2 | − 1 − lim · EX,S,T + |µ1 ∩ µ2 | n,d→∞ d d − |µ1 ∩ µ2 | d − |µ1 ∩ µ2 | − 1 −1 2 2 −1 2 γ p (1 − p) γ p = + . 2 2 (γ − p ) γ − p2 lim

n,d→∞ d

□

S.3

Contents Related to Model Free Stability Guarantees

This appendix proves stronger stability guarantees in which the expectation is taken only over the algorithmic randomness.

S.3.1

Additional Experiments

We adopt the MARSadd setting from Friedman (1991), now a standard benchmark for random forests (Mentch and Zhou, 2020; Curth et al., 2024). We generate X from Unif([0, 1]d ) and set 1

y = 0.1e4X +

4 2 1 + e−20(X −0.5)

+ 3X 3 + 2X 4 + X 5 + N(0, σ2 ).

(S.1)

We choose n = 250, d = 10, and σ = 1. We set B = 256, q ∈ {0.1, 0.2, 0.3, 0.4, 0.5}, and T ∈ [7]. In Figures S.8 and S.9, we conduct a parallel set of experiments to those for RF, which brings similar conclusions. We adopt the following simulation setup. An orthogonal design matrix X ∈ Rn×d is generated by applying QR decomposition to a matrix with i.i.d. standard normal entries. The true coefficient vector is set to β ∗ = d−1/2 (1, . . . , 1) ∈ Rd , and the noise terms are drawn independently as εi ∼ N(0, σ2 ). We fix the parameters as follows: n = 250, d = 200, and σ = 1. The number of bagging rounds is set to B = 256. We vary the feature subsampling ratio q ∈ {0.05, 0.1, 0.3, 0.5} and the number of forward selection steps T ∈ [20]. Larger values of T are not considered, as none of the chosen q values satisfy the regime qT ≲ 1 beyond this range. 49

0.007

0.008

q=0.05 q=0.1 q=0.3 q=0.5

0.006

q=0.05 q=0.1 q=0.3 q=0.5

0.007 0.006

0.005

0.006

0.003

2 2

2 2

0.004

0

0

0

2 2

0.005 0.004

0.004

0.003

0.002

0.002

0.001

0.002

0.001 0.000

0.000

1

5

10 T

15

20

0.000

1

5

(a) Weak signal

15

20

1

0.004 0.003 0.002

0.005

0.004 0.003 0.002

5

10 T

15

20

1

5

(d) Weak signal

10 T

15

0.003 0.002

20

1

0.01425

5

10 T

15

0.0135 0.0130

0.01400

0.0125

MSE

0.0145

MSE

0.01375 0.01350 0.01325

0.0140

0.0120 0.0115

0.01300

0.0135

5

10 T

15

0.01250

20

q=0.05 q=0.1 q=0.3 q=0.5

0.0110

0.01275

1

1

5

10 T

15

0.0105

20

1

(h) Mixed signal

(g) Weak signal q=0.05 q=0.1 q=0.3 q=0.5

0.013

15

20

0.011 0.010 0.009

q=0.05 q=0.1 q=0.3 q=0.5

0.012

Empirical MSE

Empirical MSE

0.010

10 T

0.013

0.012

0.011

5

(i) Strong signal q=0.05 q=0.1 q=0.3 q=0.5

0.013

0.012

Empirical MSE

20

(f) Strong signal

q=0.05 q=0.1 q=0.3 q=0.5

0.01450

0.0150

MSE

0.004

(e) Mixed signal

q=0.05 q=0.1 q=0.3 q=0.5

0.0155

20

0.000

0.000

1

15

0.001

0.001

0.001

10 T

q=0.05 q=0.1 q=0.3 q=0.5

0.006

Generalization gap

0.005

0.005

5

(c) Strong signal

q=0.05 q=0.1 q=0.3 q=0.5

0.006

Generalization gap

Generalization gap

0.006

10 T

(b) Middle signal

q=0.05 q=0.1 q=0.3 q=0.5

0.007

0.000

q=0.05 q=0.1 q=0.3 q=0.5

0.008

0.011 0.010 0.009 0.008

0.009 0.008

0.007

0.008

1

5

10 T

(j) Weak signal

15

20

1

5

10 T

15

(k) Mixed signal

20

1

5

10 T

15

20

(l) Strong signal

Figure S.7: The feature instability, generalization gap, generalization error, and empirical risk of RFS under different signal types. Each cross corresponds to 100 repetitions. We set σ = 0.1.

50

∗ We next explore the connection between feature instability and generalization. Define βstrong = −1/2 −1/2 200 ,...,5 , 0, . . . , 0) ∈ R to be a sparse coefficient with 5 nonzero positions, which represents (5 ∗ ∗ ∗ ∗ ∗ ∗ a strong signal. Let βweak = d−1/2 (1, . . . , 1) and βmixed = (βweak + βstrong )/∥βweak + βstrong ∥2 represent a weak and a mixed signal, respectively. Note that this sense of strong and weak depicts a different aspect of learning hardness opposed to the SNR (Zhou and Mentch, 2023; Liu and Mazumder, 2025). We compare the relationship between feature instability, generalization gap, generalization error, and empirical risk in Figure S.7, respectively, in four rows.

• Although there is no theoretical result establishing a relationship between the generalization gap and feature instability, we observe a similarity in their trends, regardless of whether the signal is concentrated or spread across the features. • The generalization-risk curves, however, have different patterns under different signal types. Observing the second last row of Figure S.7, when the signal is weak, it is suggested to use a smaller subsampling ratio. In contrast, when the signal is strong, using a larger q and a smaller T is a better choice. This observation aligns with past experience that the bagging provides a regularization effect, which is more promising when the signal is weak (LeJeune et al., 2020; Mentch and Zhou, 2020). • In practice, the generalization error is unobservable. Practitioners see only the empirical risk in the last row. In all signal types, the empirical risk decreases as T grows, i.e., the model capacity becomes larger. Thus, in practice, hyperparameter tuning involves balancing between the empirical risk and the generalization gap, where the latter can take evidence from feature instability. 1.0

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

0.5 0.4

0.6

MSE

w w0 22

0.8

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

0.6

0.3

0.4

0.2 0.2

0.0

0.1

1

2

3

4 T

1 Pd

(a) d

j=1 ∥w − w

5

6

0.0

7

1

2

3 1 Pd

−j 2 ∥2

(b) d

4 T

j=1 | f (x) − f

5 −j

6

7

(x)|22

Figure S.8: Feature-instability metrics for RF. Left: comparison of the feature importance between the feature-bagged RF and the non-bagged RF. The real lines with crosses correspond to the simulated values over 100 repetitions under different feature subsampling ratios. The dashed lines are the theoretical upper bound (4.7). The black solid line is the feature instability of non-bagged RF, i.e., a single tree. Right: comparison of the prediction value between the feature-bagged RF with different subsampling ratios. We compute over an independent test set of size n.

S.3.2

Proofs of Results in Section 4.1

Proof of Proposition 4.1. The proposition is a special case of Proposition 4.3, where rad(W) = rad([−M, M]) = M. □ 51

Proof of Proposition 4.3. The proof is a direct extension of that of Proposition 4.1, except that here we are considering the norm in general spaces. We use ν to denote a subsampled set of features. h Let eventi − j − j b−w b , for which we similarly define w b = Eν A(D− j ; ν) . E j,ν = { j < ν}. Denote the quantity L j = w b = Eν [A(D; ν)]. We also denote wν = A(D; ν). As a result, we can express L j Recall that we have w as h i L j =Eν A(D; ξ) − A(D− j ; ξ) h   i =Eν Eν wν − wν | E j,ν h i  
1 = Eν Eν wν − wν 1{E j,ν } . P(E j,ν ) Note that we took uniform sampling without replacement. Thus, P(E j,ν ) = 1−q. Since Eν [Eν [wν ] − wν ] = 0, we have h i  
1 Eν Eν wν − wν 1{E j,ν } P(E j,ν ) i h  
 1 Eν Eν wν − wν 1{E j,ν } − (1 − q) . = 1−q

Lj =

Then, the sum of norms has d X

b−w b− j ∥2H = ∥w

j=1

d X

∥L j ∥2H

j=1

=

d * X

L j,

j=1

h i+  
 1 Eν Eν wν − wν 1{E j,ν } − (1 − q) 1−q H

d X (i) =

hD E i
   1 Eν L j , Eν wν − wν 1{E j,ν } − (1 − q) H 1−q j=1

d X

hD   E i   1 Eν L j 1{E j,ν } − (1 − q) , Eν wν − wν H 1−q j=1 * +  d  X      1 Eν  L j 1{E j,ν } − (1 − q) , Eν wν − wν  = 1−q j=1

=

H

Here, (i) holds since L j is a constant and thus is independent of ν. Applying the Cauchy-Schwarz inequality, we have v u u u t

 2    X d h i   1   b−w b− j ∥2H ≤ · Eν ∥Eν [wν ] − wν ∥2H ∥w Eν  L j 1{E j,ν } − (1 − q)     1−q j=1 j=1 H v u u u 2  t  X     d rad(W) . ≤ Eν  L j 1{E j,ν } − (1 − q)   1−q

d X

j=1

H

Here, the last step is true because       2 2 Eν Eν wν − wν H = inf Eν w − wν H ≤ inf w∈W

 sup

w∈W w′ ∈W

52

 2 w − w′ H ≤ rad2 (W),

(S.2)

where the first equality is due to the fact that expectation is the minimizer of squared loss. We calculate the remaining term by noticing  s(s−1)  h i    d(d−1) if j1 , j2 , c c Eν 1{E j1 ,ν }1{E j2 ,ν } =    s if j1 = j2 d

and thus   s(s−1) s2  h i h i   d(d−1) − d2 c c c c 2 Cov 1{E j1 ,ν }, 1{E j2 ,ν } = Eν 1{E j1 ,ν }1{E j2 ,ν } − q =  2    ds − ds 2

if j1 , j2 , if j1 = j2 .

Then we expand (S.2) by  2   X d d    X h i  c c  = Eν  L j 1{E j,ν } − (1 − q) ⟨L , L ⟩ Cov 1{E }, 1{E } j j 1 2 j1 ,ν j2 ,ν   j=1

j1 , j2

H

! X ! d s(s − 1) s2 s s2 = ⟨L j1 , L j2 ⟩ − − + ⟨L j , L j ⟩ d(d − 1) d2 d d2 j ,j j=1 X 1

=

2

d X

2

Lj

j=1

≤

H

! ! X d s(s − 1) s2 s s(s − 1) − + ⟨L j , L j ⟩ − d(d − 1) d2 d d(d − 1) j=1

d X d q(1 − q) ∥L j ∥2H . d−1 j=1

Bringing back this into (S.2), we have d X

rad(W) ∥L j ∥2H ≤ 1−q j=1

v u t

d X d q(1 − q) ∥L j ∥2H , d−1 j=1

which leads to d

1X rad2 (W) q . ∥L j ∥2H ≤ d j=1 d−1 1−q □

S.3.3

Related Contexts and Proofs of Results in Section 4.2

Proof of Proposition 4.4. The proof is done by recursively applying equation (S.7) and Proposition 4.3. □ Our final goal is to bound the feature instability measure: d h − ji 2 1X Eξ [wT ] − Eξ wT , H d j=1

where T is some preset horizon length. We first consider a fixed j ∈ [d] for the convenience of derivation.

53

−j

−j

Let ∆wt := wT − wT and write h i h −j i 2 h i −j −j 2 . = Eξ S ξt ◦ · · · ◦ S ξ1 (w0 ) − Eξ S ξt ◦ · · · ◦ S ξ1 (w0 ) Eξ ∆wt H

H

(S.3)

−j

In the above quantity, the discrepancy between wt and wt accumulates as t grows as each step introduces −j additional discrepancy due to the difference between S ξt and S ξt . The recursive nature of the procedure indicates that we should analyze the instability using a peeling strategy. Specifically, at the t-th recursive step, we decompose the instability into two terms. The first term, referred to as the peel at step t, deals with the heterogeneity caused by the removal in the t-th step, formally defined as h i h −j i j peelt = Eξ S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) − Eξ S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) .

(S.4)

The second term is the error accumulated along 1, . . . , t − 1, referred to as the pith at step t, is defined as h −j i h −j i j −j −j pitht = Eξ S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) − Eξ S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) .

(S.5)

−j

Compared to ∆wt−1 , the pith at step t can be expressed as h −j i j −j pitht := Eξ S ξt ◦ ∆wt−1 , −j

−j

which applies the random operator S ξt to the change ∆wt−1 from the previous step. We will then show the one-step condition (4.3): j

pitht

H

h i −j ≤ (1 + δt ) Eξ ∆wt−1

H

,

(S.6)

where the inflation factor 1 + δt depends on the algorithm A and the geometry of the space H. Applying the triangle inequality, we bound the recursive quantity from (S.3) as d d h i 1X 1X j j 2 −j 2 Eξ ∆wt peelt + pitht = H H d j=1 d j=1 d

≤

2 1 X j j peelt + pitht H H d j=1

d h i 2 1 X j −j ≤ peelt + (1 + δt ) Eξ ∆wt−1 . H H d j=1

Then, by Minkowski’s inequality, we obtain the recurrence bound: v u t

d

h i 1X −j 2 Eξ ∆wt ≤ H d j=1

v u t

d

1X 2 peel j H + (1 + δt ) d j=1 | {z } | Proposition 4.3

v u t

d

h i 1X −j 2 Eξ ∆wt−1 . H d j=1 {z }

(S.7)

step t−1

The first term on the right-hand side corresponds to the instability due to one-step fitting, which can be bounded via Proposition 4.3. The inequality (S.7) forms an inhomogeneous linear recurrence, which admits a closed-form solution. This is formalized in the following proposition. We conclude with a brief remark on the peeling strategy. For instance stability, a closely related

54

analysis exists for stochastic gradient descent (SGD) (Hardt et al., 2016; Lei and Ying, 2020), which was noted earlier in this section as a key motivating example. That analysis relies on assumptions analogous to Assumption 1 and the inflation condition (4.3). As a result, the corresponding instability bounds for SGD also exhibit cumulative error growth with the number of update steps.

S.3.4

Derivation of Feature Instability for Non-bagged RFS and RF

We derive the feature instability without bagging. We first consider RFS. Notice that without bagging, every wT includes T features without doubt. Thus, Eξ [wT ] has T positions being one and others being zero. Moreover, the selected features are the T features with the largest y ⊤ X j X j⊤ y, as suggested by h(S.8).i Thus, as long as j is in the group of T −j features with largest y ⊤ X j X j⊤ y, the term Eξ [wT ] − Eξ wT has one position being 1 and the other h − ji 2 being -1. Thus, Eξ [wT ] − Eξ wT = 2. There are in total T out of d such terms, and thus d h − j i 2 2T 1X . = Eξ [wT ] − Eξ wT d j=1 d

For random forest, things are the same. Eξ [wT ] has in total T positions being one and others being zero. Removing j-th feature results in a Eξ [wT ] with the positions in j-th column being zero, and h − ji 2 P j equally many entries in other columns being one. Thus, Eξ [wT ] − Eξ wT = 2 · Ti=1 (Eξ [wT ])i , which leads to d d T X h − ji 2 1 X 1X 2T j . Eξ [wT ] − Eξ wT 2· (Eξ [wT ])i = = d j=1 d j=1 d i=1

S.3.5

Derivation of the Affine Random Forest Encoding Map

Proof of Proposition 4.6. Fix the query point x. Let Ik,r (x) be the dyadic interval in coordinate k that contains xk after r midpoint splits along coordinate k. Define Z Ak,r (x) := 2r fk∗ (u) du, r = 0, . . . , T. Ik,r (x)

If rk (w) = ∥wk ∥1 , then under the additive model and the uniform design, d X

Ey|X,w [ f (x)] = w

Ak,rk (w) (x).

k=1

Because the encoding is column-filled, Ak,rk (w) = Ak,0 +

T X

wr,k {Ak,r − Ak,r−1 }.

r=1

Therefore K(w) =

d X k=1

Ak,0 +

d X T X k=1 r=1

55

wr,k {Ak,r − Ak,r−1 }.

Thus K is affine in w. In particular, K(w1 ) − K(w2 ) =

d X T X (w1,r,k − w2,r,k ){Ak,r − Ak,r−1 }, k=1 r=1

which is linear in w1 − w2 . Since this is a fixed finite-dimensional linear functional, it is bounded by a constant times ∥w1 − w2 ∥F . □

S.3.6

Proof of the Random Forward Selection Instability Bound

Proof of Theorem 4.5. (i) We first bound the radius of Wt . Note that the space of Wt depends on the input wt−1 . Recall that wt−1 ∈ [0, 1]d , and suppose wt−1 has e t − 1 positive positions, where the collection e where of positive positions is denoted as W. Then there is a natural choice of w ej = w

1 1( j < W) + 1( j ∈ W). d −e t+1

For any wt , there holds 1 2 1 − +1=1− ≤ 1. 2 (d − e t + 1) d −e t+1 d −e t+1 j<W

e − wt ∥2H = ∥w

X

inf

sup

Thus, w − w′ H ≤ sup

w∈Wt w′ ∈Wt

w′ ∈Wt

e − w′ H ≤ 1 w

holds uniformly for all t and ξ. (ii) We next prove the contraction bound. Throughout this part, without loss of generality, we assume that the feature scores are ordered as y ⊤ X j1 X j1 ⊤ y ≥ y ⊤ X j2 X j2 ⊤ y

if j1 ≥ j2 .

Under the orthogonal-design assumption X ⊤ X = nId , let bw rt−1 = y − XEwt−1 β t−1 be the residual after fitting the currently active set. For every unselected feature j < supp(wt−1 ), adding feature j decreases the residual sum of squares by RSS(wt−1 ) − RSS(wt−1 + e j ) =

2 1 1  j⊤ X rt−1 = y ⊤ X j X j⊤ y. n n

Therefore, conditional on the candidate set, RFS selects jt ∈

argmax

y ⊤ X j X j⊤ y,

j∈νt \supp(wt−1 )

with deterministic tie-breaking. We first record the only structural fact about the pathwise difference that is used below. Auxiliary fact. For every fixed realization of ξ1:t and every j ∈ [d], define −j

j

−j

Dt := S ξt ◦ · · · ◦ S ξ1 (w0 ) − S ξt ◦ · · · ◦ S ξ1 (w0 ). 56

(S.8)

Then j

Dt ∈ {0, e j } ∪ {e j − ea : a ∈ [d] \ { j}}. j

Equivalently, Dt is nonnegative on the j-th coordinate, nonpositive on all other coordinates, and has at most one negative coordinate. To prove the auxiliary fact, use induction on t. The statement is trivial at t = 0. Suppose it holds at j step t − 1. If Dt−1 = 0, then the two paths have the same active set except that the reduced procedure cannot select feature j. Therefore, either both paths select the same feature, or the full path selects j while the reduced path selects the best available alternative. The new difference is thus either 0, e j , or j e j − ea . If Dt−1 = e j , then the two paths have the same active set outside coordinate j, and the reduced −j update S ξt never operates along j; hence the non- j selection is the same and the difference remains j

e j . Finally, if Dt−1 = e j − ea , then the full path contains j whereas the reduced path contains a, with all other active coordinates identical. At the next step, the only possible negative discrepancy is either still attached to a, transferred to a lower-ranked feature b < a, or removed if the full path selects a. No second negative coordinate can be created, because both paths select at most one new feature at the current step. This proves the auxiliary fact. Recall that the pith term in the proof of Proposition 4.4 is h −j i (@) :=Eξ S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) h −j i −j −j − Eξ S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) . We need to show h i −j −j ∥(@)∥2 ≤ (1 + δt ) Eξ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) − S ξt−1 ◦ · · · ◦ S ξ1 (w0 )

2

.

(S.9)

Define the pathwise previous-step difference −j

j

−j

Dt−1 := S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) − S ξt−1 ◦ · · · ◦ S ξ1 (w0 ), and define the pathwise pith difference j,pith

Dt

−j

−j

−j

−j

:= S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) − S ξt ◦ S ξt−1 ◦ · · · ◦ S ξ1 (w0 ).

Then h j,pith i (@) = Eξ Dt . j

By the auxiliary fact, on the event Dt−1 , 0, the difference is either e j or e j − ea for a unique a , j. Let I j = [d] \ { j}. For a ∈ I j , define j

ua := Pξ1:(t−1) (Dt−1 = e j − ea ), and define j

ρ := Pξ1:(t−1) (Dt−1 = e j ),

α := ρ +

X

ua .

a∈I j

Let u = (ua )a∈I j , viewed as a row vector indexed by I j , and embedded in Rd by putting zero on the j-th

57

coordinate. Then h j i Eξ Dt−1 = αe j − u.

(S.10)

j

If Eξ [Dt−1 ] = 0, then α = 0 and u = 0, and the desired contraction is trivial. Hence we assume below j that Eξ [Dt−1 ] , 0. For a, b ∈ I j , define the conditional transition kernel  j,pith  (t, j) j Mab := P Dt = e j − eb | Dt−1 = e j − ea ,

(S.11)

(t, j)

with Mab = 0 if the conditioning event has probability zero. Write M = M (t, j) for short. This is a conditional kernel. It is not the unconditional matrix E[1(Ea,b )]. Therefore, the following representation follows from the law of total probability and does not require any independence between the previous discrepancy event and the current transition event: h j,pith i Eξ Dt = αe j − uM.

(S.12)

−j

j

−j

Indeed, if Dt−1 = e j , the two inputs to S ξt differ only on coordinate j, and S ξt does not operate along j, −j

j

so the pith difference remains e j . If Dt−1 = e j − ea , then after applying S ξt , the negative coordinate is either transferred to some b ∈ I j , or disappears. This gives exactly (S.12). We now bound the matrix norms of M. First, for every a ∈ I j , X  j,pith  j Mab = P Dt , e j | Dt−1 = e j − ea ≤ 1. (S.13) b∈I j

Hence ∥M∥∞ ≤ 1. Second, Mab = 0 whenever b > a. To see this, suppose the previous negative coordinate is a. A transition from a to b means that, under the same candidate set at step t, the mixed path selects a, while the reduced path selects b. If b > a, then b has score no smaller than a. Since b is available to the mixed path whenever such a transition is considered, the mixed path would not select a before b. Thus the negative coordinate can only stay at a, move to some b < a, or disappear. It remains to control the column leakage. Define ηt := max{q(1 + q), q2 t}. We claim that for every b ∈ I j , d X

Mab ≤ ηt .

(S.14)

a=b+1 j

Fix a > b. Condition on both the event Dt−1 = e j − ea and the past σ-field generated by ξ1:(t−1) . Let A −j and A− j be the active sets of the two inputs to S ξt . Since the two paths differ only by replacing a with j, their union has cardinality at most t. Moreover, a ∈ A− j ∪ A, and the transition a 7→ b requires a and b to be included in the fresh candidate set. After fixing a and b, the remaining s − 2 sampled features must avoid every unselected feature r > b. Otherwise, the reduced path would select a feature with index larger than b, or the mixed path would not select a. Therefore, the remaining s − 2 sampled features must

58

lie among the b − 1 lower-ranked features and the at most t − 1 active features in (A ∪ A− j ) \ {a}. Thus, b+t−2 s−2

d  ,

Mab ≤

(S.15)

s

with the convention that the numerator is zero if b + t − 2 < s − 2. For s = 1, a transition a 7→ b with a , b is impossible, so the column-leakage bound is trivial. Hence assume s ≥ 2. Let k = d − b. Summing (S.15) over a > b gives d−k+t−2 d X s−2 (S.16) Mab ≤ k d . a=b+1

s

If k ≤ t, then k

d−k+t−2

d−2

s−2 d  s

≤ k d  = k

s−2 s

s(s − 1) ≤ q2 t. d(d − 1)

If k > t, then d−k+t−2 k

s−2 d  s

d−k+t−2 s(s − 1) s−2 = k d−2 d(d − 1) s−2

s(s − 1) k−t ≤ k 1− d(d − 1) d−2

! s−2 .

The one-dimensional function k−t h(k) := k 1 − d−2

! s−2 ,

k > t,

is maximized either at the boundary k = t, which yields the already bounded term q2 t, or at the stationary point d+t−2 k∗ = . s−1 Substituting k∗ and using k∗ > t gives the elementary bound k∗ − t s(s − 1) k∗ 1 − d(d − 1) d−2

! s−2 ≤ q(1 + q).

Consequently, (S.14) holds. Combining (S.13) and (S.14), we have ∥M∥∞ ≤ 1, Indeed, for a fixed column b, X a∈I j

Mab = Mbb +

X

∥M∥1 ≤ 1 + ηt .

Mab ≤

a>b

X c∈I j

59

Mbc + ηt ≤ 1 + ηt .

Therefore, ∥M∥22 ≤ ∥M∥1 ∥M∥∞ ≤ 1 + ηt .

(S.17)

We now complete the contraction argument. From (S.10) and (S.12), j

Eξ [Dt−1 ] and j,pith

Eξ [Dt

]

2 2 2 2

= α2 + ∥u∥22 , = α2 + ∥uM∥22 .

By (S.17), ∥uM∥22 ≤ ∥M∥22 ∥u∥22 ≤ (1 + ηt )∥u∥22 . Hence j,pith

Eξ [Dt

]

2

≤ α2 + (1 + ηt )∥u∥22

2

j

= Eξ [Dt−1 ]

2 2

+ ηt ∥u∥22 .

Since X

∥u∥2 ≤ ∥u∥1 =

ua ≤ α,

a∈I j

we have

∥u∥22 α2 + ∥u∥22

1 ≤ . 2

Therefore, j,pith

Eξ [Dt

j

]

Eξ [Dt−1 ] Taking square roots and using

j,pith

Since (@) = Eξ [Dt

2

ηt 2 ≤1+ . 2 2 2

√ 1 + x ≤ 1 + x/2 for x ≥ 0, we obtain  ηt  j,pith j Eξ [Dt ] ≤ 1 + Eξ [Dt−1 ] . 2 2 4

], this proves (S.9) with 1 1 δt = ηt = max{q(1 + q), q2 t}. 4 4

Finally, we apply Proposition 4.4. Since rad(Wt ) ≤ 1, when T ≲ q−1 , we have for every t ≤ T , q2 t ≤ q(1 + q), and thus δt =

q(1 + q) . 4

60

Therefore, (4.4) gives  T !T −t 2 d X  2   1X q(1 + q) q  −j  Eξ [wT ] − Eξ [wT ] ≤  1+ 2 d j=1 4 (d − 1)(1 − q) t=1 2    1 + q(1+q) T − 1  q 4   =    q(1+q)  (d − 1)(1 − q)  4

 2 !T   q(1 + q) 16 =  1 + . − 1 2 4 q(1 + q) (d − 1)(1 − q) When T ≲ q−1 , this simplifies to d

T2 q 1X −j 2 Eξ [wT ] − Eξ [wT ] ≲ · . 2 d j=1 d−1 1−q This proves the theorem. □

S.3.7

Proof of the Random Forest Instability Bound

Proof of Theorem 4.7. Let Er,k ∈ RT ×d denote the matrix with one at entry (r, k) and zero elsewhere. We first prove the radius condition. Fix a feasible state wt−1 and define the one-step candidate set n o Wt (wt−1 ) := S ξt (wt−1 ) : ξt is an admissible realization of the tree randomness . Under the max-edge dyadic rule, every possible next state has the form wt−1 + EN k +1,k t−1

for some k ∈ M(wt−1 ). If the implementation allows no split for a particular realization of ξt , the next state is simply wt−1 , which only makes the following bound easier. Therefore, for every wt ∈ Wt (wt−1 ), ∥wt − wt−1 ∥F ≤ 1. Since the center in the definition of radius may be any point in the ambient space, choosing wt−1 as the center gives
rad Wt (wt−1 ) ≤ 1. Taking the supremum over feasible wt−1 yields rad(Wt ) ≤ 1. −j

−j

The same argument applies to the feature-removed transition S ξt after replacing M(wt−1 ) by M− j (wt−1 ). We next prove the contraction condition. Fix a removed feature j ∈ [d]. Let wt−1 = S ξt−1 ◦ · · · ◦ S ξ1 (w0 )

61

be the local encoding after running the full tree for t − 1 steps, and let −j

−j

−j

wt−1 = S ξt−1 ◦ · · · ◦ S ξ1 (w0 ) be the corresponding local encoding after running the tree with feature j removed. The two runs are coupled using the same randomness ξ = (ξ1 , . . . , ξt ). Define the one-step increment of the feature-removed transition by −j

−j

Uξt (w) := S ξt (w) − w. −j

For any feasible w, the increment Uξt (w) is either zero or a standard basis matrix Er,k with k , j. Hence, for any two feasible states w and w′ , −j

−j

Uξt (w) − Uξt (w′ )

F

≤

√ 2 1{w , w′ }.

If w = w′ , the two increments are identical because the split rule and the tie-breaking rule are deterministic. Now let h i −j −j ∆t−1 := Eξ wt−1 − wt−1 . We claim that   −j −j Pξ wt−1 , wt−1 ≤ ∥∆t−1 ∥F . Indeed, if the full run has never selected feature j during the first t − 1 steps, then the full run and the feature-removed run have the same local cell, the same max-edge candidates after excluding j, and the same split scores on all features in [d] \ { j}. By deterministic tie-breaking, the two local paths are identical. Therefore, n o n j o −j wt−1 , wt−1 ⊆ Nt−1 ≥ 1 . Since the encoding is column-filled and the feature-removed tree has zero in the j-th column, the (1, j) −j entry of ∆t−1 is  −j   j  ∆t−1 = Pξ Nt−1 ≥ 1 . 1, j

Thus,    j   −j  −j Pξ wt−1 , wt−1 ≤ Pξ Nt−1 ≥ 1 = ∆t−1

−j

1, j

≤ ∥∆t−1 ∥F .

Using the decomposition −j

−j

S ξt (w) = w + Uξt (w), we obtain h −j i −j −j Eξ S ξt (wt−1 ) − S ξt (wt−1 ) F h −j i −j −j −j = ∆t−1 + Eξ Uξt (wt−1 ) − Uξt (wt−1 ) F   −j −j −j −j ≤ ∥∆t−1 ∥F + Eξ Uξt (wt−1 ) − Uξt (wt−1 ) F   √ −j −j ≤ ∥∆t−1 ∥F + 2 Pξ wt−1 , wt−1 √ −j ≤ (1 + 2)∥∆t−1 ∥F . Equivalently, h −j i −j −j Eξ S ξt (wt−1 ) − S ξt (wt−1 )

F

h i −j ≤ (1 + δt ) Eξ wt−1 − wt−1

62

F

with δt =

√ 2.

This proves the contraction condition required by Proposition 4.4. Combining this contraction condition with the radius bound rad(Wt ) ≤ 1 and applying Proposition 4.4 yields (4.7). The prediction bound follows from Proposition 4.6. □

S.3.8

Proof of the Finite-Bagging Instability Bound

Proof of Proposition 4.8. Note that Hoeffding’s inequality in Hilbert space (e.g., Boucheron et al. (2013)) yields s

B

1 X (b) w − Eξ [w] ≤ B b=1

3 rad (W)2 log ((d + 1)/δ) B

(S.18)

with probability 1 − δ/(d + 1), and B 1X

B b=1

s h

w(b),− j − Eξ w− j

i

≤

3 rad (W)2 log ((d + 1)/δ) B

(S.19)

with probability 1 − δ/(d + 1) for each j ∈ [d]. Applying union bound, we have (S.18) and (S.19) hold simultaneously for all j ∈ [d] with probability at least 1 − δ. Then, with Minkowski’s inequality, we have 2 2 √ X d B B d B 1 X 1 X (b) 1 X (b),− j 3 1 X (b) w − w w − Eξ [w] ≤ d j=1 B b=1 B b=1 d j=1 B b=1 √ X d h i 2 3 + Eξ [w] − Eξ w− j d j=1 2 √ X d B h i 3 1 X (b),− j −j − Eξ w w + d j=1 B b=1 2 √ X d B 3 1 X (b) ≤ w − Eξ [w] d j=1 B b=1 √ 6 3 rad (W)2 log ((d + 1)/δ) . + B

□

63

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

0.025

0.150 0.125

|f(x) f j(x)|2

0.020

|f(x) f j(x)|2

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

0.175

0.015 0.010

0.100 0.075

2

3

4 T

5

6

7

1

2

(a) Weak signal

4 T

5

6

0.0

7

1

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

0.8

Generalization gap

0.5 0.4 0.3 0.2

3

4 T

5

6

7

6

7

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

1.4 1.2

0.6

0.4

1.0 0.8 0.6 0.4

0.2

0.1

0.2

0.0

1

2

3

4 T

5

6

0.0

7

1

(d) Weak signal

1.08

2

3

4 T

5

6

0.0

7

2.4

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

2.2

1.8

3

4 T

5

6

7

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

4

2

1

(g) Weak signal

2

3

4 T

5

6

7

1

2

(h) Mixed signal q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

1.0 0.9

5

3

1.4

1.05

4 T

5

1.6

1.06

3

6

MSE

MSE

1.07

2

2

(f) Strong signal

2.0

1

1

(e) Mixed signal

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

1.09

MSE

2

(c) Strong signal

Generalization gap

0.6

3

(b) Mixed signal

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

0.7

0.3

0.1

0.000

1

0.4

0.2

0.025

0.000

Generalization gap

0.5

0.050

0.005

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

0.6

|f(x) f j(x)|2

0.030

3

4 T

5

6

7

(i) Strong signal q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

2.25 2.00

q=0.1 q=0.2 q=0.3 q=0.4 q=0.5

6 5

0.7 0.6

Empirical MSE

Empirical MSE

Empirical MSE

1.75

0.8

1.50 1.25 1.00

0.5

0.75

0.4

0.50

1

2

3

4 T

5

(j) Weak signal

6

7

4 3 2 1

1

2

3

4 T

5

(k) Mixed signal

6

7

1

2

3

4 T

5

6

7

(l) Strong signal

Figure S.9: The feature instability, generalization gap, generalization error, and empirical risk of RF under different signal types. Each cross corresponds to 100 repetitions.

64