Stabilizing black-box model selection with the inflated argmax
Authors: Melissa Adrian
Jake A Soloff
Rebecca Willett
Venue: Accepted by TMLR ()

Abstract: Model selection is the process of choosing from a class of candidate models given data. For instance, methods such as the LASSO and sparse identification of nonlinear dynamics (SINDy) formulate model selection as finding a sparse solution to a linear system of equations determined by training data. However, absent strong assumptions, such methods are highly unstable: if a single data point is removed from the training set, a different model may be selected. In this paper, we present a new approach to stabilizing model selection with theoretical stability guarantees that leverages a combination of bagging and an ''inflated'' argmax operation. Our method selects a small collection of models that all fit the data, and it is stable in that, with high probability, the removal of any training point will result in a collection of selected models that overlaps with the original collection. We illustrate this method in (a) a simulation in which strongly correlated covariates make standard LASSO model selection highly unstable, (b) a Lotka–Volterra model selection problem focused on identifying how competition in an ecosystem influences species' abundances, (c) a graph subset selection problem using cell-signaling data from proteomics, and (d) unsupervised $\kappa$-means clustering. In these settings, the proposed method yields stable, compact, and accurate collections of selected models, outperforming a variety of benchmarks.

Published in Transactions on Machine Learning Research (12/2025)
Stabilizing black-box model selection with the inflated argmax
Melissa Adrian
maedrian@uchicaag.. edu
Data Science Institute, University of Chicago
Jake A. Soloff
soloff@umich..cm edu
Department of Statistics, University of Michigan
Rebecca Willett
willetf@uchicago.ede edu
Department of Statistics, University of Chicago
Department of Computer Science, University of Chicago
NSF-Simous National Institute for Theory and Mathemalics in Biology
Reviewed on OpenReview: https: openrewiew. net/ forum? id= gSD
Abstract
Model selection is the process of choosing from a class of candidate models given data. For instance, methods such as the LASSO and sparse identification of nonlinear dynamics (SINDy) formulate model selection as finding a sparse solution to a lincar system of equations determined by training data. However, absent strong assumptions, such methods are highly unstable: if single data point is removed from the training set, different model may be selected. In this paper, present theoretical stability guarantees that leverages combination of bagging and an "inflated" argmax operation. Our method selects a small collection of models that all fit the data, and it is stable in that, with high probability, the removal of any training point will result in a collection of selected models that overlaps with the original collection. We illustrate this method (a) a simulation in which strongly correlated covariates make slandard LASSO model selection highly unstable, (b) Lotka-Volterra model selection problem focused on identifying how competition in an ecosystem influences species' abundances, (c) graph subset selection problem using cell-signaling data from proteomics, and (d) unsupervised k-means clustering. In these settings, the proposed method yields stable, compact, and accurate collections of selected models, outperforming a variety of benchmarks.
approach to stabilizing model selection with
1 Introduction
Model selection is typically formulated as procedure to identify the model that "best" reprosents data from among set of candidate models. For instance, scientists may choose among many interpretable models to identify salient variables, interactions, policies, with the goal of developing fundamental understanding of the natural or social phenomenon underlying experimental data. Examples of model selection include variable selection (choosing a subset of covariates from a larger pool that best explains the response variable or label), selecting decision tree, selecting the order or memory of an autoregressive model, selecting a kernel function for a support vector machinc, selecting the number of clusters in k-mcans, seleeting the number of factors in factor analysis, and more. The issue of model selection is common problem in variety of domains, including bioinformatics (Saeys et al., 2007), environmental studies (Effrosynidis and Arampatzis, 2021), and psychology (Vrieze, 2012).
In this paper, we argue that constraining a procedure to relurn a single model does not adequately express any uncertainties in the model selection process, and therefore, returning a set of selected model(s) is more appropriate than only ever returning one. An idcal method would return a single model often as possible while guaranteeing some degree of stability, ie., the selection should not be too sensitive small changes in the data (Yu, 2013).

---

Published in Transactions on Machine Learning Research (12/2025)
Especially in model selection scenarios where the goal is to interpret the model to learn something funda- mental about. the studied process, trustworthiness in the results is paramount. A necessary condition for trustworthy model selection procedures is stability. Conventional methods for model selection provide ad hoc procedures that do not provide theorctical stability guarantecs without strong assumptions on the data or model. Our goal is to provide a method that is adaptive to the underlying uncertainly of the model selection process without placing strong assumptions on the data distribution or the models themselves.
In this paper, propose stabilizing model selection, loveraging ideas from stabilized classification (Soloff et al., 2024b), which exhibits desirable theoretical guarantees. We provide comparisons of our proposed approach to widely used conventional model selection methods as well as explore the tradc-offs of stability, accuracy, and interpretability of the model selection procedures.
1.1 Our contributions
In this work, we allow model selection procedures to return a set of candidate models, conveying uncertainty about which model is best while returning a single model as often as possible. We formulate a notion of model selection stability based on the idea that our claims about which model is best should remain logically consistent when we drop a small amount of data.
We develop framework for stabilizing black-boz model selection procedures by combining bagged model selection with the inflated arymar (Soloff et. al., 2024b), which was formulated in the context of multiclass classification. The key insight of this paper is that if we formulate model selection as a multi- class classification problem, where each class corresponds to one of the candidate models, then we can directly use the inflated argmax and inherit its corresponding theoretical stability guarantees for model selection, which are agnostic the model class and data distribution. Additionally, the stability guarantee holds even if the "true" model is not considered in the class of candidate models. We emphasize that our approach is highly generic and can be applied to many base model selection tasks, such variable selection, constructing decision trees, selecting the order of an autoregressive model, and more.
We show empirically that the inflated argmax can improve black-box model ion stability with examples in clustering, decision trees, and variants of variable selection tasks. For variable selection (a special case of model sclection), we show this result for data with correlated covariates, sctting known to yield instabilitics in standard model selection algorit hmss. We compare against and outperform numerous conventional model selection approaches, described in $2.2. These conventional methods, unlike our approach, are not adaptive to the underlying level of uncertainty of the model selection procedure.
2 Model selection
More formally, model selection focuses on selecting one or more models from among a set of candidate models, denoted M A model selection procedure is a map M from data D to set of selected models M C M The output M = MM(D) is a nonemply subset of M+. I always outputs a singleton (i.e., |M = 1), M is called simple.
For example, in variable selection (a special case of model selection), the data consists of n variable-response pairs D = {1 v= where x; E RRd. We seek to identily relevant variables (e.g., variables of x most predictive of y.), so the model class M = 2id] is the power set of [d] (l....,d) d). However, in more general settings, the set of all candidate models M+ can contain any class of models. For instance, in graphical model selection, M+ is the set of all graphs G on d vertices, where an edge kj,k) in G encodes conditional dependence between variables j and l (Drton and Perlman, 2004; Friedman et al., 2008). In dynamical system identification, + is a set of possible differential equations governing the dynamics underlying noisy data (Brunton et al., 2016). In econometrics, M could be a collection of competing time series forecasts (Masini et al., 2023). To keep our discussion general, the model class is an abstract, countable set M+.
I Note that this general formulation of the model selection problem, and does not assume that there is true data generating model in M+.
2

---

Published in Transactions on Machine Learning Research (12/2025)
Algorithm Bagged model selection (Breiman, 1996a;b)
input Model weighting algorithm A; data D; ensemble size B; bag size for b I...., do
Construct D6 by sampling uniformly (with or without replacement) K times from D Compute model weights u = A(D*)
end for
output Ensemble weights v = B Σx= 17
We decompose the model selection procedure M into two stages
M=SoA,
(1)
where A maps data D to weights a = A(D) M where for each m M+ tm is the weight assigned to model m E M This weight may be loosely interpreted the probability of the candidate model being the best model given the data among the choices in M+ The selection rule S takes these weights é as input and returns a non-empty set of selected model(s) M S M+ that provide the best performance as a function of the training data, where performance may reflect some combination of fit to data and structural properties of the model.
2.1 Stage 1: Assigning weights to candidate models with A
Base algorithm A. A base algorithm A is any off-the-shelf model weighting algorithm. Many model selection methods lead to weight vectors á that place all weight on a single model (i.c., the algorithm A returns only one model, so w is a vector of zeros with an entry of 1 corresponding to the selected model). Examples of such algorithms are LASSO for variable selection (Tibshirani, 1996), SINDy for model discovery for dynamical systems (Brunton et al., 2016), and graphical LASSO for graph subset selection (Fricdman et al., 2008), all of which are experimentally explored in this work. A known issue of these procedures is that these algorithms produce unstable model selections. One mechanism that has been shown to help its stability is bagging.
Bagged algorithm ÂK.B. We denote ÁK,B as the algorithm resulting from applying a bagging techmique on a base algorithm A, where K is the number of samples per bag and B is the number of bags. The output of the weighting algorithm AK, must be bounded. For simplicity, we assumne that the weights can be normalized such that t belongs to the probability simplex A A1 2 Bagging averages an ensemble of weight vectors fit on different randomly sampled bags; see Algorithm Sampling with replacement is known bagging, and sampling without replacement is known as subbagging.
For base algorithm A that places all weight on a single model, bagging A to compute o.. simply counts the fraction of bags where m was selected. Even if the set of candidate models is infinite (i.e., |M+| co), the weights v will have at most B nonzero entrics, so bagging is still tractable. Generally, increasing the number of bags B used in ÂK, provides a more stability, as we empirically show in Figures 12 and 13.
2.2 Stage 2: Conventional selection rules S
We discuss conventional selection approaches S and their limitations.
The standard argmax. The default selection rule selects the best model(s) M E Eargmax,, Φn (M may contain more than one model in the case of exact ties for the top weight.) If the weights contain near- maximizers, the argmax can be sensitive to small perturbations of w, meaning that a different model may be selected with small perturbation to the training data.
Top-k. Instead of returning the very best model(s) using the argmax, an alternative is to return the top k models that receive the largest weight, where is user-specified. Similar to the argmax, top-k fails to
2he probability simplex Mh/-1 consists of nonnegative weights that sum one: Emexe Wmm
3

---

Published in Transactions on Machine Learning Research (12/2025)
automatically adapt to the inherent uncertainty in the model selection process. In particular, if v has lar more than k near-maximizers, top-k(1) will be unstable. On the other hand, if i has fewer than & mear-maximizers, top-k returns a set that is unnecessarily large.
Stability selection. In the context of variable selection, the stability selection method (Meinshausen and Bühhlannn 2010) selects a single model based on whether each covariate is included in the model with probability at least т. Concretely, consider a vector of weights á associated with candidate models M M+, where [M+ 20d since each of the d covariates is cither included or not for each model in M Define the marginal variable probabilities II,
Emglduem Wm, for j =I,.... L....,d d, where j E m corresponds to variable being included in model m. In other words, 11 corresponds to the frequency of variable being included in any selected model. Stability selection returns (1) == { [d] fi, 2 for some use-specified threshold т. Similar to the standard argmax (without ties), ip, (w) only returns a collection of variables, which logether constitute a single model. This selection procedure can be unstable if some marginal probabilities II, are near the threshold т. SA presents more detailed comparison between this method and our method, particularly in the setting of correlated variables.
2.3 Model selection stability
The stability of M corresponds to how much the output of M may change with small perturbations to D for any dataset D. Note that this cannot be estimated empirically because we cannot test all possible datasets D with finite computational resources. We adapt the definition of stability in Soloff et al. (2024b) to the context of model selection.
Definition 2.1. A procedure M has model selection slabilily 1 - 5 at sample size n if, for all datasets D with n samples,
n 1 MPV =2 56
(2)
n
where M = M(D) and M = M(D))) for each i E [n].
In Definition 2.1, D'Vi refers to dataset D with the i-th sample removed. In plain language, the inequality in (2) quantifies stability in terms of whether the sets of selected models overlap when we drop a single observation at random. Note that this notion of model selection stability holds for any "black box" model selection method of the form (1), in contrast to prior works that focus on particular scttings, such as variable selection (e-g., Meinshausen and Bühlmann (2010)) or k-means clustering (e-g., Ben-David et al. (2007)). While this definition of stability may seem weak (i.c., in practice, much larger perturbations to the data could be made), we will see empirically in our results that this minimal change to the data can have a substantial impact on model selection stability. Moreover, stability according to this criterion is necessary for most stronger notions of stability to be satisfied.
3 Related work
Model selection has a rich history in the statistics literature. Ding et al. (2018) provides a recent overview of model selection techniques. Various aspects of model selection have been studied in previous work, including hyperparameter selection via cross-validation (Raschka, 2020), preventing overfitting (Cawley and Talbot, 2007), sclection critcria (Rao et al., 2001), among others. In this work, we are particularly interested in the aspect of model selection stability, and bagging has been key tool to address instability in model selection.
Returning model selection set. Returning set of models to acknowledge inherent uncertainty in model sclection is an old idca. Early literature on subset selection for lincar regression (Spjøtvoll, 1972) focused on selecting the true model with high probability. Breiman (2001) coined the terms "Rashomon set" to refer to a complete set of nearly equally performing models (e.g., Xin et al. (2022) finds the full Rashomon set for decision trees), and the "Rashomon effect" to refer to the phenomenon that multiple, vastly distinct models may all perform equally well on out-of-sample data for particular problem. More recent works extend the idea of model confidence set to more general modeling settings (Hansen et al., 2011; Zhang
4

---

Published in Transactions on Machine Learning Research (12/2025)
et al., 2024). These works have stochastic model for the data and construct confidence set containing the best-fitting model with high probability. By contrast, our framework places no stochastic assumptions on the data, instead targeting a more modest goal of stability for the selected set.
Stability in variable selection. There is extensive literature on stability in variable selection; see, e.g., Khaire and Dhanalakshmi (2022) for a recent overview. Work in this area has overwhelmingly focused on "simple" procedures, returning a single subset of variables, and stability is measured based on the similarity of the selected subsets of variables. By contrast, Definition 2.1 is not satisfied by simple procedures unless they return the same model most of the time.
Bagging to improve stability. Bagging is an important tool for model selection and variable selection, and has been widely applied lo stabilize model selection procedures (Breiman, 1996a;b). For example, stability selection thresholds the marginal bootstrap inclusion probabilities of each variable (Meinshausen and Bühlmann, 2010; Shah and Samworth, 2013). Additionally, BayesBag (Bühhmann, 2014) directly averages posterior distributions over different resamplings of the data. Huggins and Miller (2023) analyze the accuracy and stability properties of the BayesBag approach specifically for Bayesian model selection. Soloff et al. (2024a;c) shows that bagged weights are stable, meaning that dropping single sample cannot change the bagged weights d by much in the Euclidean norm. Bagging can thus express model uncertainty, but, as we show empirically in $5, selecting the most, frequent model among multiple close contenders can still be quite unstable. In order to achieve any stability guarantees of the selected model(s), a new selection rule is needed, as discussed in $4.1.
Stabilizing classification. Our work builds on ideas from stable classification. Soloff et al. (2024b) introduces framework for sct-valued classification, combining the inflated argmax with bagging to guarantee stability for classification.
4
Our approach
We now provide a brief overview of our pipeline for stabilizing a model selection procedure M = S A.
First, we bag the weighting algorithm A, i.c., average the weights fit on different random samples from D. We fully define bagging in Algorithm
Next, we pass the bagged weights d through the inflated argmaz, which selects set of near- maximizers of u that is robust to small perturbations of the weights. We define the inflated argmax for model selection in Definition 4.1.
In $4.2, we provide a stability guarantec for this pipeline that can be applicd to any algorithm A and any data set D. We first discuss our pipeline in detail.
4.1 The inflated argmax for model selection
In this section, we propose an alternative selection procedure S based on the inflated argmax. Soloff et al. (2024b) introduced the inflated argmax to guarantee stability in the context of multiclass classification. A key insight of this paper is that if we formulate model selection as multi-class classification problem, where each class corresponds to one of the candidate models, then can directly use the inflated argmax and inherit its corresponding theoretical stability guarantees for model selection. With this insight, the inflated argmax, an alternative model selection criteria S that we denote by argmaxf, will allow us to decide which subset M of M to return.
Definition 4.1 (Inflated argmax). For m E M let
R {
E E A Wm 2 max Won' + #m
(3)
/2 S
where 1|M+1|1 is a probability simplex. For any w E 1 and E > 0, the inflated argmax is defined as
argmax (w) {m E M+ dist(w, R < <s},
(4)
5

---

Published in Transactions on Machine Learning Research (12/2025)
Figure 1: Visualization of argmax" when |A/+|=
[0.0.1]
3. The weights assigned to each model A can be mapped onto a 3- 3-dimensional probability sim- plex shown. The example weight vector _ 5+ 5 lands in an area that is within e dis- tance of regions RS and Rž, which are the re- gions corresponding to returning models my and Ma, respectively. Because this weight vector is an "uncertain" region based on & for distin- guishing whether model m2 or my provides better fit to the data, argmax will return both models: {m2, ma)- In contrast, weight vector = lg 9 3] would land in the area labeled Rg, and argmax will return my only.
R2
of w with dist(w, RD) < e and dist(w, Rg) <e
4
9
Rf
R2
[1.0,0]
w = [0,1,0]
where dist(w, Rm.) = inf Re, 11g wll, and |I- is the Euclidean norm.
A visualization of a simple example of argmaxx where more than one model is returned is shown in Figure 1. An intuitive explanation of Definition 4.1 that m E argmna (us) indicates that m would have the largest weight in w by some margin if the weights w were slightly perturbed, where both the size of the margin and the weight perturbations scale with This definition naturally allows for the number of returned models to adapt to the uncertainty in the model selection: for example, if the top two competing models are close in probability, then the inflated argmax will return both models for some choice of small On the other hand, if the top two competing models have large enough gap in their two weights, the inflated argmax will only return the top model for a choice of small e. The paramcter & is not something to be tuned. It is chosen to satisly a user-specified stability lolerance; see $4.2.
4.2 Stability guarantee
In this section, we state a theorem guaranteeing the stability of our approach based on our definition of stability in Definition 2.1.
Theorem 4.2. (Adapted from Soloff al. (2024b, Theorem 17).) For any model selection procedure S A, our method argmax o AK,B satisfies model selection stability instability level
1 ( E²
)
1
P
16c
8
+
[M+
12
11
B
where p
K. for subbagging and 1 (1 - 2 for bagging.
The proof immediately follows from the proof presented in Soloff et al. (2024b). The main use for this theorem is that the user can choose tolerable worst-case instability 8 for their model selection, input the known parameters of the problem (i.e., |M+| candidate models, K samples per bag, total samples, and B bags), and directly solve for the that should be used for argmax'.
Recall that e determines how much we inflate the argmax, so when it's smaller, the inflated argmax is closer to the standard argmax. Accordingly, when è is smaller, Theorem 4.2 gives larger d, reflecting a weaker stability guarantee. Further, note that the instability level 8 decreases (ie., stability increases) by decreasing the bag size K and increasing the number B of bags. The size of our model class |M+1 has only a small influence on the level of instability of our procedure. For example, when there are two models, |M+ 2, the factor 1 M- 1 giving constant-factor boost in stability over the case where M+ is countably infinite, where 1 1
J becomes 1. This theoretical result is the first to provide theoretical guarantees for model selection stability that is agnostic to the underlying models, base algorithm, and data distribution.

---

Published in Transactions on Machine Learning Research (12/2025)
NOTATION
SELECTION CRITERION
argmax
Model(s) with maximum weight assigned by A
argmax AK,B
Model(s) with maximum weight assigned by Ak,B
top-k o AK,B
models with & largest weights assigned by AK.B
ip, AK.B
Model with variables with inchusion probability 2 T in AK.B
argmax" AK,B
e-inflated argmax Ak,B
Table 1: Table of notation for each model selection algorithm. A refers to base algorithm mapping a dataset D to empirical predicted probabilities v for each possible subset, and Ax,n is the subbagged version of A for bags with K randomly selected (without replacement) samples.
5 Experiments
Our experiments assess and compare the stability of five different. model selection approaches, described in Table 1. We perform experiments in three settings: identifying important variables in a synthetic dataset that contains a small collection highly correlated variables (see Appendix C), identifying the latent Lotka- Volterra governing differential equations from noisy trajectory data generated from these equations, and constructing graph to represent cell-signaling pathways of 11 proteins. In the first two settings, generate synthetic datascts for N trials, each contaminated with independent Gaussian noise, allowing us to examine the stability across collection of datasets where we know the ground truth m².
In each experiment, we compute weights & across the set candidate models M using base algorithm A and a subbagged version of A denoted Ax,B. Based on (2), we compute the stability for trial j E [N] using
6, n n HM, Ma =0).
(5)
This empirical measure of stability measures, for trial j, the proportion of LOO selected models M) for E [n] that have overlap with the set of selected models M, returned by training with access to the entire dataset. We also compute the empirical cumulative distribution function (CDF)
Σ 1 56,
(6)
N
as a function of E [0, 1] to highlight variation across trials with different random seeds (shown in Figure 3 and Figure 8). Curves higher on the plot are better in the sense that the corresponding method achieves model selection stability 8 (Definition 2.1) for larger fraction of trials, and our primary intorest is in small 5 (i.e., high stability).
We additionally compute weighty- accuracy, inspired by Zaffalon et al. (2012), using
1 1m e M,
(7)
M1
where m is the "true" data generating model in the synthetic experiments. This measure down-weights the accuracy by the number of models returned |M,]. For example, if M 1, and the true model contained in M (i.e., m E M), then both the accuracy and the weility- accuracy are 100%. However, if [M = 1, 000 and m e M, the accuracy 100%, but the weighted accuracy 0.1%. The utility- weighted accuracy mcasure aims to identify returned model sets that are both accurate and interpretable, and we argue that smaller number of returned models more interpretable. This metric used to compare methods in the center column of Figure 3.

---

Published in Transactions on Machine Learning Research (12/2025)
kolterr Dynamica P. 00 1.31
Emplrical Comelation Matro a Deginee 2 Payrumial Library
Figure 2: (Left) Lotka-Volterra trajectories for one simulated dataset. The time span ranges from 0 to 16, and there are 21 time points with predator-prey information contaminated with (0,, 0.05) distributed noise. (Right) Empirical correlation matrix for the variable library input into SINDy and E-SINDy. We use the library (c, x(), (u(()))),, (29), u(2), (u(2))9), where is a constant.
4.5--
1.35
Covariste
0:5
w/t..
Tirro
Observer
Coveratte
5.1 Model discovery of Voltera dynamics
This experiment focuses on discovering Volter- dynamics from data using SINDy. The Voltera equations were first derived to model chemical reactions in the seminal work (Lotka, 1910) and later studied in the now classic use case of predator-prey dynamics (Totka, 1925). The Lotka- equations are
(1)()) = (1)- (u) w/a) (1) = (1) + Cu (1))
(8a)
(8b)
where u1)) is the population density of prey, (²)) is the population density of prodators, a()) and a(2) are the population growth rales, and corresponds to time. The parameters a, 3, 7, C > 0 respectively refer lo the prey's growth rate, the predator's 's impact on the prey's death rate, the predator's death rate, and the prey's impact on the predator's growth rate (Lotka, 1925). We now imagine we do not know our system is governed by the Lotka- Volterra equations and wish to discover the governing model from data. Our concern is in recovering the correct terms rather than estimating the parameters.
5.1.1
Base algorithm A: Sparse Identification of Nonlinear Dynamical Systems (SINDy) with Sequentially Thresholded Ridge Regression (STRidge)
One application of sparse regression is the discovery of governing equations of system from observations. Specifically, we seek to learn the governing equations of a system of the form
d dt u(t)= f(((t),
(9)
where u(t) (1), ((),,((((()))))) (1)) (() is d-dimensional state vector containing quantities of interest time t, and f is the function describing how the derivative of u(t) evolves with time. The Sparse Identification of Nonlinear Dynamical Systems (SINDy, Brunton et al. (2016)) algorithm represents a differential equation of the form (9) as a weighted sum of terms in a library of functional forms,
(t) (w)).
(10)
where @(u) RTXp (n samples of p functionals) is the library of functionals of u(t) (e.g., polynomial transformations), and E Roxd are the fitted coefficients. =\t is assumed to be sparse, so our goal is to find which elements of E should be nonzero, which will immediately give the estimated structure of the governing equations. SINDy finds a sparse solution to E by solving a sparse lincar regression problem; see $B for more details on sparse regression.
As a bagged extension of SINDy, Ensemble SINDy (E-SINDy, Fasel et al. (2022)) creates an ensemble of model fits, and with this ensemble, the final model is chosen via thresholding inclusion probabilities of each library term based on tolerance This approach corresponds to applying stability selection (Meinshausen and Bühlmann, 2010; Shah and Samworth, 2013) to SINDy, using sequentially thresholded ridge regression (STRidge, Rudy et al. (2017)) instead of the LASSO for A. E-SINDy aggregates the coefficients by selecting covariates whose inclusion probability is greater than some prespecified value T E (0, 1). Of the variables with an inclusion probability greater than 7, the corresponding coefficients are set to the average fitted values across bags. The collection of variables with a nonzero coefficient correspond to the selected model. We refer to this ensemble approach as ip, Ax AR where the base algorithm A is STNDy solved with STRidge.
8

---

Published in Transactions on Machine Learning Research (12/2025)
5.1.2 Data generation
We generate synthetic datasets containing n 21 time points of abundance data for predator and prey populations, solve from the VolterVo equations and contarninate with Gaussian noise. This data is generated match key characteristies of real world predator- prey dataset from the Hudson Bay Company (Hewitt, 1921) in terms of the periodicity of observations in relation to the dynamics. N 100 trials (i.e., independent datasets) are independently generated according to the same data generation process. We provide further details of this data generation process in $D.1. We visualize one set of trajectories in Figure 2.
To construct our library of functions, we selected degree-two polynomials as our covariates, and for one dataset, we visualize the corresponding covariate correlations Figure 2. As seen in the correlation matrix, the covariate structure contains highly correlated variables (e.g., x(1) and (u()))),, which impacts the model selection stability of SINDy. In our experiments, we solve the SINDy linear regression problem with STRidge penalization (Rudy et al., 2017) our base algorithm A, as in Fasel et al. (2022).
5.1.3 Results
We see that across various choices of E, the inflated argmax provides enhanced stability while the same time yielding small and accurate sets of selected model(s). This boost in utility-weighted accuracy is due to the adaptivily of the inflated argmax: this method will output single model as often as possible, whereas for example, top-k will always output k model sclections, leading to comparatively lower utility in many cases. This point is further confirmed in the plot of the median number of models returned vs. worst-case instability in Figure 3. For a fixed level of worst-case instability empirically seen, we can visualize the empirical CDF assess the distribution of instability values 6, across trials j E [N]. We see Figure 3 that the inflated argmax provides a distribution of instability values much closer to (i.c., higher stability) compared to all other conaentional approaches investigated.
In summary, Figure 3 shows that (a) selection methods that return single model (argmax and ipr, even with subagging the base algorithm) are less stable than methods that return one or more selected models (i.e., the inflated argmax and lop-k), (b) the inflated argmax provides higher utility-weighed accuracy scores for similar worst-case stability levels as top-k, meaning that the inflated argmax is returning smaller sets containig the correct model more often and (c) the inflated argmax provides smaller sets of returned models compared to top-k for a wide range of worst-case stability values. We additionally plot the stability across the parameter B in SE
Lotka-Volterra
10-
05
5000
0.6
Urikty pegghee accurnt
Modion Medioni
>WA 0.4
02
0.
0.-
3.0
02
0.8
0.6
5.8
0.6 Woust Woust-Cose Instoliii By
0.0
Werst-Case rotatelity (197a 6
A (Brunton
ar.
(Mainshausen Bünimann)
argmaxi (Curs)
Figure 3: Lotka- Volterra results. (Left column) Empirical CDF for the stability measures 8y, computed in (5), across j E [N] for each S. We chose the parameters 7 = 0.63, k 2, and = 0.09 since these values yield a weigtty-weihhted accuracy approximately 0.3 for each S. (Center column) weility accuracy, computed in (7); versus the worst-case instability across 5y. (Right column) Median number of models returned across versus the maximum d, across j E [N] for each S. (Center and right columns) We plot the range of values & E [1,...,6] 6), and (0, 1) and (0, 1) with approximately evenly spread values across their supports.
9

---

Published in Transactions on Machine Learning Research (12/2025)
Top Graph
Top Graph
Figure Visualization of the top two graph struc- tures selected via subbagged gLASSO with 10,000 bags. The red connection shown in the right graph highlights that this connection is the only difference between the top and top 2 seleeted graph struc- tures. The top graph was selected for 9.31% of bags, and the top graph was selected for 8.54% bags.
plcp
plog
pattser2
poitueiig
PMC PS0
PRC P30
5.2 Graph selection on flow cytometry data
We consider a model selection problem where the model to select the sparse connection structure within a graph. The base algorithm A we consider for this problem is Graphical LASSO (gLASSO), which is developed in Friedman al. (2008). The gLASSO estimates the structure of graph by performing a maximum likelihood estimation of the precision matrix @ subject to an penalty of the estimated precision matrix. More preciscly, the gLASSO assumes the data ~ NO θ-*), where Rd across samples E [n].
Empirical Correlation Matrix
reaf
1.00
prek
0.75
p.c
0.5
and finds an estimate è such that
HF2
MPS
0.00
argmax log det θ tr(EO) X||O||1.
(11)
palbv173
-035
RCI
-0.50
where ||- ||a is the entry-wise I norm, tr(-) is the trace, Σ is the sample covariance matrix, and À is the penalization hyperparameter. In our ex- periments, we select À via cross-validation as described in $F, and hold this quantity fixed across our experiments.
-0.75
ROK
-1.00
pod
pmok
pas
Figure 4: Empirical correlation matrix of the 11 proteins. Data from Sachs et al. (2005).
We compute LOO stability results from a flow cytometry datasel in Sachs et (2005) (labeled "od3ee22icam2++u1125 (nodes) and n 759 cells (samples). Figure shows the empirical corre- lation structure of these 11 proteins in this dataset, which clearly shows grouped relationships among proteins.
containing =1 proteins
SELECTION METHOD LOO INSTABILITY AVG. LOO SET SIZE
argmax
0.453
1.00
ipo.5 Âx,B
0.013
1.00
argmax AK,
0.112
1.00
top-2 AK.B
0.008
2.00
(Ours) argmax
AK,
0.008
1.58
Table 2: Table of stability and the average number of model structures returned across LOO trials for selecting the structure of the inverse precision matrix for a flow cytometry dataset taken from Sachs et al. (2005). In this experiment, is the graphical LASSO (Friedman et al., 2008), ÅK, the bootstrapped graphical LASSO with B 10, 000 and K 700.
5.2.1 Results
Since we have one dataset to compute the LOO stability, N = 1. We compute stability results for the algorithms gLASSO (A) and subbagged gLASSO (Ax,n), where K = 700 and B = 10, 000, for the selection methods argmax, ipo.s, top-2, and 0.00 The value for for the inflated argmax was chosen so that the LOO instability closely matched the LO0 instability of top-2, and r in ip, was chosen so that the number of selected connections was approximately the the same as the number of connections given by the top selected subbagged graph. These results shown in Table 2.
10

---

Published in Transactions on Machine Learning Research (12/2025)
As seen in Table 2, our method, argmax0002 provides the lowest LOO instability while minimizing the average number of graphs selected across the 759 TOO model fits. We also note that the percentage of model fits on bags that resulted in the top 2 models shown in 5 are relatively close: 9.31% for the top graph structure versus 8.54% for the sccond most selected graph structure, which suggests a reasonable degree of uncertainty in the "best" model. Due to this uncertainty the selected model, selection methods that. relurn only one model (i.c., argmax and ip-) cannot sufficiently account for this uncertainty in the top model. Since the top two graph structures differ by one connection, this kind of result can provide evidence for follow up experiments, for example, to elucidate the relationship between the pakts473 and pmek proteins.
Though there is no agroed-upon "ground truth" graph connection structure for this real world example, we can compare the identified connections between proteins with those established in the literature. For example, Sachs et al. (2005) visualizes connections between proteins based on scientific evidence, which include PTP3-pleg, praf-pmek, PKA-praf, PIP2-PIP3, PKA-p38, and others, where some of these listed connections are mediated via intermodiate proteins. These validatod connections are seen Figure We note, however, that, hyperparameter À in (11) largely impacts the accuracy by controlling the number of connections retained in the graph. In experiment, this parameter was chosen in purely data-driven way in $F, but different value may have been chosen with domain expertisc.
5.3 Selecting the number of clusters for k-means
In this section, we present results on simulated data of an unsupervised k-means clustering (MacQueen, 1967) example where the task is to identify the number of clusters x in the data. This examples ralustrates our algorithm's use non-variable selection types of problems. Another example on decision trees can be found in SH. We generate N = 100 independent datasets, cach with sample size n = 30 with a true number of clusters x =3 3. Further details of our synthetic data generation can be found in SG.
5.3.1 Results
In this setting, each model corresponds to different number of clusters K, so the model selection problem is to choose In our experiments, our base algorithm A is deterministic computation of the "elbow" of the sum of squared distances versus number of clusters K, which is heuristic approach introduced in Thorndike (1953) and commonly used in unsupervised clustering tasks. For each k, clusters are partitioned with the Lloyd-Forgy algorithı (Forgy, 1965), the standard algorithm for k-means. Further details of our data generation and implemented methods can be found in &G.2. We construct AK,B B by bagging this base algorithm by computing B 10, 000 bags, each with K 25 samples per bag. Figure 6 shows the main results of our experiment. We notably do not include results for the model selection method ip, AK,B since that method is only applicable to variable selection problems.
The leftmost panel shows the empirical distribution of stabilities across methods, with fixed values of e = 0.3 and h = 2 for argmax" and top-k. Choosing k- = 2 provides 100% stability, meaning there are seemingly only two reasonable choices for the number of clusters k. The argmax0.3 showed an empirical worst-case instability of 0.067 and an average returned set size 1.18 across LOO trials, which provides high stability with minimal returned models. The middle plot of Figure 6 shows that our approach provides the best weighte accuracy and stability compared to competing methods. In this setting, accuracy is measured based on how often the truc = 3 is returned in M, across all j = .,.... N trials). This result shows that the inflated argmaxx returns sets that are both accurate and small as possible at each fixed è value. This point is further confirmed by looking at the rightmost plot of the median number of returned models versus worst- case instability. Across wide range of choices of worst-case instabilities, the inflated argmax consistently returns interpretable sets.
3In this paper, we use K to represent the number of samples bag, correspond the number of models returned by t selection procedure, and x to represent the number of clusters clustering problem.
11

---

Published in Transactions on Machine Learning Research (12/2025)
6.5
2.7
0.00
2.00
0.52
Vadaa Pocas
2.23
(
0.1
0.42
0.1
0.0 24
1.0
600
Marst basterily
Mayi-Cooe bolibiityy
argmax-.
argmezn Are
bop-ke Ale
angmae! Ae (Ours)
Figure 6: Unsupervised k-means clustering results. (Left column) Empirical CDF for the stability measures 8,, computed in (5), across [N] for each S. We chose the parameters k = 2, and 0.3. (Center column) Ueility- accuracy, computed in (7), versus the worst-case insta- bility across 65. (Right column) Median number of models returned across j versus the maximum 95 across j E [N] for each S. (Center and right columns) We plot the range of values E [1, and e [0.005, 0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.95,
6 Discussion
Our framework provides multiple opportumities for further investigation, including exploring the interplay of stable model selection and identifiability, as well as stabilizing model selection without bagging, therefore reducing computation. Moreover, this work primarily focuses on providing theoretical stability guarantees for model selection, but in practice, practitioner would also be interested in maximizing model performance. Additional future work could explore maximizing performance subject to stability guarantees.
We believe another interesting direction is developing similar theoretical stability results for other measures of stability aside from Definition 2.1. Our definition of stability is relatively weak, reflecting a bare minimum criterion necessary (though not sufficient) for many other notions of stability. Moreover, a limitation of our stability measure is that it may favor larger sets since two returned sets need overlap by only one model to achieve a score of for that LOO trial. For example, if |M,] and |M]] are large but [M, N MA the LOO stability for trial i would be 1, which may not reflect what one would colloquially deem as "stable." An interesting alternative stability measure to theorctically analyze would be one that computes the percent overlap between elements of M, and M)."
A natural follow up question to our approach is what a practitioner should multiple models are re- turned. As noted in $5.2.1, top models that substantially overlap can prompt directed follow-up experiments. Moreover, the size of the set of returned models may be a useful indicator of model uncertainty for a prac- titioner. The decision for how to handle multiple models should be carefully considered and align with the practitioner's 's ultimate goals (e.g., model interpretation or predictive accuracy).
The combination of bagging and the inflated argmax described in this paper can be applied to any black-box model selection procedure and offer stability guarantees described in Theorem 4.2. We define a notion of stability that distinct from other works, one that computes the proportion of LOO selected models that are disjoint from the selection made by training with access to the full data. This in contrast to, for example, covariale-level stability in linear models.
Our experiments illustrate that, in addition to providing theoretical stability guarantees, our method out- performs ad hoc procodures. Specifically, choosing the most frequently selected model across bags (i.c., argmax o AK,n) provides limited stability gains relative to an unstable base algorithm. Selecting the top-k models across bags enhances stability but not adaptive to the underlying empirical uncertainty of the model selection process; k models will always be returned to the user. Choosing the returned model basod on marginal inchusion probabilities (ip=) of individual covariates limited covariate stability tasks and presents hyperparameter selection challenges (i.e., how to select r). In contrast, the inflated argmax adapts
12

---

Ben-David, S., Pál, D., and Simon, H. U. (2007). Stability of k-means clustering. In Bshouty, N. H. and Published in Transactions on Machine Learning Research (12/2025)
to the number of models returned based the underlying uncertainty, provides approach for empirically protecting against, a worst-case instability, and provides competitive weility-weighted accuracies (i.e., returns the correct model in the smallest possible M).
Acknowledgments
We gratcfully acknowledge the support of the National Science Foundation via grant DMS-2023109 and the NSF-Simons National Institute for Theory and Mathematics in Biology via grants NSF DMS-223545) and Simons Foundation MP-TMPS-00005320. MA gratefully acknowledges the support of the NSF Graduate Research Fellowship Program under Grant No. DGE-17460455. JS gratefully acknowledges the support of the Office of Naval Research via grant N00014-20-1-2337,, and the Margot and Tom Pritzker Foundation.
References
Gentile, C., editors, Learning Theory, pages 20-34, Berlin, Heidelberg. Springer Berlin Heidelberg.
Breiman, L. (1996a). Bagging predictors. Machine Learning, 24:123-140
Breiman, L. (1996b). Hcuristics of instability and stabilization in model selection. The Annals of Statistics, 24(6):2350 2383.
Breiman, L. (2001). Statistical modeling: the two cultures (with comments and a rejoinder by the author). Statistical Science, 16(3):199- 231.
Breiman, L., Friedman, J., Olshen, R. A., and Stone, C. J. (1984). Classification and Regression Trees. Chapman and Hal//CRC, edition.
Brunton, S. L., Proctor, J. L., and Kutz, J. N. (2016). Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences of the United States of America, 113(15):3932-3037
Büihlmann, P. (2014). Discussion of big Bayes stories and BayesBag. Statisticnl Science, 29(1):91-9) -44
Cawley, G. C. and Talbot, N. L. C. (2007). Preventing over-fitting during model selection via bayesian regularisation of the typer-parameters. Journal of Machine Learning Research, ((31):8418861
de Silva, B., Champion, K., Quade, M., Loiseau, J.-C., Kulz, J., and Brunton, S. (2020). Pysindy: A python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5149):2104.
Ding, J., Tarokh, V., and Yang, Y. (2018). Model selection techniques: An overview. IEEE Signal Processing Magazine, 35(6):16 34.
Drton, M. and Perlman, M. D. (2004). Model selection for Gaussian concentration graphs. Biometrika, 9133)5591-602.
Effrosynidis, D. and Arampatzis, A. (2021). An evaluation of feature selection methods for environmental data. Ecological Informatics, 61:101224.
Fasel, U., Kutz, J. N., Brunton, B. W., and Brunton, S. L. (2022). Ensemble-sindy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control. Proceedings of the Royal Sociely A: Mathematical, Physical and Engineering Sciences, (78(2260)20021990)
Forgy, E. W. (1965). Cluster analysis of multivariate data: Efficiency versus interpretability of classifications. Biometrics, 21:768-766
Friedman, J., Hastie, T., and Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatisties, 9(3):432 441.
13

---

Published in Transactions on Machine Learning Research (12/2025)
Hansen, P. R., Lunde, A., and Nason, J. M. (2011). The model confidence set. Econometrica, 79(2):453-49) 497. Hewitt, C. G. (1921). The conservation of the wild life of Canada. C. Scribner's Sons, New York.
Hoerl, A. E. and Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1):55-67.
Huggins, J. Н. and Miller, J. W. (2023). Reproducible model selection using bagged posteriors. Bayesion Analysis, 1811)799.104
Kaptanoglu, A. A., de Silva, B. M., Fasel, U., Kaheman, K., Goldschmidt, A. J., Callaham, J., Delahunt, C. B., Nicolaou, 7. G., Champion, K., Loiseau, J.-C., Kutz, J. N., and Brunton, S. L. (2022). Pysindy: A comprehensive python package for robust sparse system identification. Journal of Open Source Software, 7(699:3994.
Khaire, U. M. and Dhanalakshmi, R. (2022). Stability of feature selection algorithm: A review. Journal of King Saud Universily Compuler and Information Sciences, 34(4):1060-14 1073.
Lotka, A. J. (1910). Contribution to the theory of periodic reactions. The Journal of Physical Chemistry, 14(3):271-3 2744.
Lotka, A. J. (1925). Elements of Physiccul Biology. Williams and Wilkins Company, Baltimore, Maryland.
MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. Berkeley Symp. on Math. Slatist. and Prob., 5.1:281-299
Masini, R. P., Medeiros, M. C., and Mendes, E. F. (2023). Machine learning advances for time series forecasting, Journal of Economic Surveys, 37(1):76-11) 111.
Meinshausen, N. and Bühlmann, P. (2010). Stability selection. Journal of the Royal Statistical Society, Series B (Statistical Methodology), 7444)4417-473.
Petzold, L. (1983). Automatic selection of methods for solving stiff and nonstiff systems of ordinary differ- ential equations. SIAM Journal on Scientific and Statistical Computing, 4(1):136-148.
Rao, C. R., Wu, Y., Konishi, S., and Mukerjee, R. (2001). On model selection. Lecture Notes-Monogruph Series, 38:1 64.
Raschka, S. (2020). Model evaluation, model selection, and algorithm selection in machine learning. httts://arxiv..ee www..orggabbba aaaaaass/11.// 1881..288.
Rudy, S. H., Brunton, L., Proctor, J. L., and Kutz, J. N. (2017). Data-driven discovery of partial differential equations. Science Advances, 3(4):1-6.
Sachs, K., Perez, O., Pe'er, D., Lauffenburger, D. A., and Nolan, G. P. (2005). Causal protein-signalin networks derived from multiparameter single-cell data. Science, 308(5721) (570555:553--333 -:529-599.
Saeys, Y., Inza, I., and Larrañaga, P. (2007). A review of feature selection techniques in bioinformatics. Bioinformatics, 23(19):2507-2517,
Shah, R. D. and Samworth, R. J. (2013). Variable selection with error control: another look at stability
selection. Journal of the Royal Statistical Society, Series B (Statistical Methodology), 75(1):55 80.
Soloff, J. A., Barber, R. F., and Willett, R. (2024a). Bagging provides assumption--fre stability. Journal of Machine Learning Research, (513)):
Soloff, J. A., Barber, R. F., and Willett, R. (2024b). Building stable classifier with the inflated argmax. In The Thirty-eighth Annual Conference on Neural Information Processing Systems.
Soloff, J. A., Barber, R. F., and Willett, R. (2024c). Stability via resampling: statistical problems beyond the real line. https:/ httss://arxiv.ooo httpss//arxiv.org/aae www.ii.orgg/ab//445..og
14

---

Published in Transactions on Machine Learning Research (12/2025)
Spjetvoll, E. (1972). Multiple comparison of regression functions. The Annals of Mathematical Statistics, 43(4):1076-1088
Thorndike, R. L. (1953). Who belongs in the family? Psychomelrika, 18(4):267 276.-776.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statisticnl Society, Series B (Methodologieal), 58(1):267-22)
Veleslavov, C. and Stumpf, M. P. (2020). Decision tree models and cell fate choice. 12.00......1...2.2 423622.
https:/ htsss//wwwwp http///wwwww/ww/oo/////////t///p www.xiv.org/coooeeeeeoe
Vrieze, S. I. (2012). Model selection and psychological theory: discussion of the differences between the A Akake information criterion (AIC) and the Bayesian information criterion (BIC). Physchological Methods, 17(2):228 243.
Xin, R., Zhong, C., Takagi, T., Seltzer, M., and Rudin, C. (2022). Exploring the whole rashomon set of sparse decision trees. Advances in Neural Information Processing Systems, 35:14071 14084.
Yu, B. (2013). Stability, Bernoulli, 19(4):1484 1500.
Zaffalon, M., Corani, G., and Mauá, D. (2012). Evaluating credal classifiers by uillity-discounted predictive accuracy. International Journal of Approzimate Reasoning, 53(8):1282-1331.
Zhang, T., Lee, H., and Lei, J. (2024). Winners with confidence: Discrete argmin inference with an applica- 02060.
tion to model selection.
A Example: Linear regression with correlated covariates
Consider the setting in which we have set of covariates [x(1),, xx()) .... ,...,,.......................... )[]] where
E = 2] = (1) + (3).
(12)
and we seek to learn this relationship using sparse variable selection. (1) is highly correlated with x;12) and (3) is highly correlated with (44) and 25)) (see covariance matrix in Figure 7), then subtle changes the dataset could result. in selecting any of the following models mx(x).
mx(x) =1x(1), x0)) m²(x) =(x("),, x0())
mx(x)= (1) (1))),,())), ma(x) =(x(²), ,,())))
(13)
mx(x) (2x(2), (0)) ma(x) =(x(2), (0))
If each of the six models above were selected with approximately equal frequency across bags, the output of Algorithm 1, would be approximately
à (1/6) 1,,
(14)
A selection method that finds the argmax of u would be unstable despite the stability of à derived from bagging. We soek method that outputs models corresponding to selected variable sets
{ (100), (2)) x (x) (4) 1
(15)
accurately indicating that our method suggests either In or 22 should be included in the selected model, but not necessarily both. In other words, this output reflects the fact that our model selection method cannot reliably distinguish among the six models in (13).
Limitations of the argmax. The argmax applied to in (14) will be highly sensitive to small perturbations in w. Furthermore, the single best model does not reflect the inherent uncertainty in the model selection and that several models almost equally likely.
15

---

Published in Transactions on Machine Learning Research (12/2025)
Limitations of top-k. If the user chooses k 2, only 2 models will be returned to the user, which does not reflect that any of the 6 models are approximately equally likely. This failure to adapt to uncertainty also applies to the other direction: if the user chooses k > 6, top-k will still give back & models, despite the fact that only 6 are highly likely.
Limitations stability selection The stability selection method (Meinshausen and Bühlmann, 2010) will select the model corresponding selected variables {x(1) x(("), 200), 204) زققه if 7 is sufficiently small. This result can be misleading, in that it ultimately selects model with five selected covariates opposed to indicating that we lack certainty about which of multiple two-covariate models is correct.
Benefits of the inflated argmax. For a large enough choice of e, the argmax" would return all 6 models in (13). This output can be reported as a Boolean expression as in (15), which accurately reflects the uncertainty in the model selection.
B Variable selection in linear models with LASSO
In the linear case, the data is modeled as the form
yi (2 B) + Thi ' N(0,
(16)
where o² IR a vector containing the response error variance, and 3 g is the vector of regression coefficients. In many practical scttings, the true population is hypothesized to be sparse, and the goal is to accuralely estimate which elements of 8 are nonzero.
Consider the following common model selection approach: first solve the least absolute shrinkage and selection operator (LASSO, Tibshirani (1996)) problem
| = argming
Σ U3 (z,, B))² AR(3)
(17)
n
where À is a hyperparameter controlling the sparsity of B, and the regularization term R(3) = |11811. Then select variables via
m(D) = 1118 > 0}
(18)
where 3 implicitly depends on the dataset D. In this context, the LASSO provides unstable model selections, particularly in settings with nontrivial correlations among covariates (Meinshausen and Bühlmann, 2010). Moreover, since many variable selection algorithms, including the LASSO, only return a single set of selected covariates, these methods do not directly provide an approach for assessing uncertainty in the selected model.
C Experiment: Penalized linear regression experiment with highly correlated covariates
This experiment is based on the motivating example in (12), where we have two important covariates and a small collection of highly correlated covariates.
C.1 Data generation
In our standard regression example, we simulate datasets for j E [N] trials, where N 200, according the following process:
D; N0,
(19)
where x; e Rd, N denotes a Normal distribution, and C E gaxd d is visualized in Figure Figure 7 visualizes the covariance structure of the 20 covariates in the design matrix X, where the nonzero covariance between distinct covariates is 0.99, corresponding to extremely correlated covariates. For single trial simulate dataset with d = 20 covariates and E [n] samples. We take the first and third columns as the covariates used to create the response y,
Di x + xr
(1)
(3) Uis is ~ N(0,0.0."
(20)
16

---

Published in Transactions on Machine Learning Research (12/2025)
Covariance Matrix for Covariates
Covariate Index
Covariance
0.
-02
11
Covariate Index
Figure 7: Covariance matrix for each simulated dataset. Covariate indices and 3 were used to construct the response
Regression
0.2
0.0
Worth-Case Instabllty (20
Worti-Caee ietabiityy (ress
argmaxr.4 (Tibshirani)
argmext
ip, Ac. (Meinshsusen and Bühhmann)
top-kr.Ax.
argmax" Acc (Ours)
Figure Regression results. In the left plot, we plot the stability curves for the parameters
0.3,
k = 2, and 0.8 since these values yield ntility-weightee accuracy of approximately 0.18 for each S. (Left column) Empirical CDF for the stability measures &,, computed in (5), across E [N] for each S. (Center column) Utility-weighted accuracy, computed in (7), versus the worst-case instability across 6j. (Right column) Median number of models returned across j versus the maximum 6; across j E [N] for each S. (Center and right columns) We plot the values [l....,6], and e E (0,1) and T (0, 1) with approximately evenly spread values across their supports.
where == [1, x and x; (3) respectively the first and third covariates for sample i, and represents independent Gaussian distributed obscrvation noise with mcan 0 and standard deviation 0.3 for sample í.
We solve each optimization problem (i.e., fits on the full datasets D, LOO datasets DV for E [n]] and bags D' for E [B] across trials j [N]) for å using the optimization function in (17) for R(3) ||8]|1, which is equivalent to the LASSO.
C.1.1 Optimization hyperparameter
Selecting À via cross-validation lcads to much denser models than the ground truth model, so we increase À to a level that induced more sparsity. This approach implicitly assumes we have prior knowledge of approximately how many coefficients are important the regression task, but not necessarily which are important. We select = 0.5 since this penalization level led highly sparse variable selections. We note, however, that the theory behind the inflated argmax is model agnostic, meaning that stability guarantees apply for any choice of X.
17

---

Published in Transactions on Machine Learning Research (12/2025)
Selected Covariates
0.218
0.092
0.09
Top Model.
[
0.084
0.057
Weights
0.056
0.051
0.033
Figure 9: Top-10 selected models with their corresponding weights + Σ j==1 across trials. Selected covariates for cach model are shown with a black square.
( which averages the weights
C.2 Small sample size (n= 30n=30
We generate independent datasets (i.e., (rials) for small sample size, n 30. For each of the N 200 trials, generate dataset with d 20 covariates and 30 training samples.
C.2.1
Base algorithm AA: LASSO
Our learning algorithm A is the LASSO, which corresponds to solving (17) with R(3) [| and AK,B corresponds to subbagged LASSO with B 10, 000 bags and K = 25 training samples per bag.
C.2.2 Results
Figure 8 provides an empirical CDF of ô (6), the case instability, and the median number of sets returned across N trials versus the worse-case instability, The empirical CDF across our chosen selection methods S show that the selection methods that return single model, the argmax and ipr, are the least stable. Subagging, which helps in stabilizing the estimated weights à assigned to each model m, yields marginal gains in stability when using the argmax. Top-k and the inflated argmax, on the other hand, provide enhanced stability due to the fact that these mothods are able to return more than one model. The ip, selection method fails to control the worst-case instability across a variety values of t, but able to achieve good utility-weighted accuracy. Top-k is able achieve good weility-eighhted accuracy with stable estimates. However, the inflated argmax is able attain similar stability as top-k, while achieving higher utility- accuracy across range of stability levels.
accuracy (7) versus the empirical worst-
We confirm that the inflated argmax frequently outputs single model by plotting the median number of models returned compared to the worst-case instability across various choices of e compared to various choices of
To illustrate that returning multiple models, when necessary, can be represented compactly, we borrow notation from logic to express an example of 6 models, reported in Figure Figure 9 shows that the model weights averaged across trials ( ( The top 6 models can be reported as
{ (ox()x(0)) x fz(3),((1) (1)) x (0))
which succinctly describes the model options.
C.3 Larger sample size (n 300n=300)
To illustrate the utility of the inflated argmax in a setting where this method is particularly useful, we con- struct an example where the model weights á estimated via bootstrapping suggest a high level of uncertainty in the model selection. We focus on one trial, so N 1. In constructing this example, we follow similar procedure described in (19) and (20). We set d 200, n 300, B 10, 000, K 25, 0.5, and
18

---

Published in Transactions on Machine Learning Research (12/2025)
U; ~ 0.5.. where 0.5 is the response error standard deviation for each sample i. The correlation struc- ture of covariates is extended to 200 covariates, where similarly the first 5 covariates are highly correlated with a correlation of 0.99 as in Figure 7, and z(1) and x3() are used to construct with 3 = [1,1]'
Figure 10 shows the selections of among the first 5 covariates out the 200 total for the top-10 models (covariates through 200 are not selected in any of the top-10 models) using the full dataset versus a 1.00 dataset. As an illustrative example, we choose a LOO dataset that yielded different top model compared to the top model using the full dataset.
The model weights clearly show uncertainty in the selection since the weights estimated for each of the top-10 models are very close together, suggesting that small perturbations of the dataset could impact the top model selected. The empirical stability for argmaxoAx,n using (5) yields 6, = 0.14 for this one trial j. This result means that 14% of the time in our 300 LOO trials, when a single sample was removed from the full datasct, model other than (a(²), .211) (the model selected using all 300 samples) was selected.
The inflated argmax, on the other hand, will select more than one model in this instance for relatively small due to the large uncertainty in the model weights. To concretely show this result, compute the inflated argmax's s model selection(s) and choose using the result in Theorem 4.2. Setting 5 to a fixed value, we can compute via the following expression for subbagging
K/n Kfn
E
V
(21)
n
Since |M+ = 2200, we set |ar+] to 0. Additionally, since the Monte- Carlo error term 16e² B is known to be
1
overly conservative, and since we use a large number of bags B in this experiment, we set this term equal to 0 as well. With 300, K = 25, and ő = 0.05, we compute that = 0.078.
The instability of arrmaxx 0.078 o Ax, is 0, meaning that for 0% of the 300 leave-one-out datasets, the inflated argmax returned a completely disjoint model selection. Additionally, the inflated argmax returned models, which correspond to
{ (100), (x)
(1) (3)
This set of returned models is very reasonable given the underlying correlation structure among the covariates.
Moreover, the lype of stability in this work only measures changes the stability due to LOO changes to the dataset, which is a relatively small perturbation to the data. In reality, a method would ideally provide robustness to more substantial changes in the dataset (e.g., providing the same, or approximately the same, model selection(s) with a different dataset drawn from the same underlying distribution). The inflated argmax is better equipped than the argmax in this scenario by returning six models, reflecting the large uncertainty in selecting model.
Main experimental details
In our experiments, we utilize a cluster computing system to distribute parallel jobs across CPU nodes.
D.1 VolterVolterra simulations
This section provides further details of the experiment in $5.1. We utilize the pysindy Python package (Kaptanoglu et al., 2022; de Silva et al., 2020) for their implementations of these methods to perform out experiments with SINDy.
D.1.1 Data generation
We generate data from the Lotka-Voltene differential equations described in (8) for (a, 7, (2/3, 4/3, 1, 1) for 21 evenly spread time points t in the range [0, 16) and initial condition [1, IT. Solu- tions of this ordinary differential equation were obtained with the LSODA (Petzold, 1983) solver. For each
19

---

Published in Transactions on Machine Learning Research (12/2025)
Selected Covariates
Selected Covariates (LOO)
(Full Data)
X11
x
x44 x
X41
x2
x
x4
yS
0.062
0.056
0.055
0.049
Top-10 Model Weights (W)
0.044
0.048
0.036
0.043
0.032
0.032
0.026
0.029
0.02
0.019
0.014
0.014
0.014
0.013
0.012
0.012
Figure 10: Models with the top 10 largest weights assigned via subbagging with B 10, 000 bags for one trial. The y-axis labels correspond to the model weights estimated with 10,000 bags. Covariates within each model (row) that were selected to have a nonzero estimated coefficient by LASSO are shown in black, and covariates that were not selected are shown in white. Only the first. covariates are shown since the top models did not select covariates any of the covariates ((()) ,x(7), ,..., (2200), (Left) Top-10 models using the full dataset where n 300. (Right) Top-10 models using LOO dataset where
(200))
299.
(0, 16), non-inclusive of the initial condition, we contaminate each (t) and (t) with No 0.053) Gaussian noise, where 0.05 is the standard deviation of the measurement noise level.
To obtain sparse solutions the regression problem and therefore estimate the governing equations, we solve the optimization problem with STRidge (Rudy al., 2017), which itcratively (a) solves the equation in (17) with R(3) || corresponding solving Ridge regression (Hoerl and Kennard, 1970), and (b) thresholds 3 based on w, where w is the minimum magnitude needed for a coefficient in the weight vector. Any coefficient with an estimated magnitude below is set to This process repeated until convergence criterion is met.
D.1.2 Cross-validating optimization hyperparameters
We perform grid search across two paramcters find combination that leads to low validation MSE. The 5-fold validation MSE measured
5
val-MSE(A, 5) 5 5()) o (EV(1)
Σ
[ - -u(1)
(22)
where V(v) the collection of time points included in the validation evaluation for fold v E 1,2,.... 5;
a(t; A, w) a vector of solutions, containing (t;) X, and (t: X,w), the model selected by SINDY using optimization parameters À and w at time t; u(t) is a vector of the observed values for (t)) and u (22) (t)
at time The trajectories were temporally split into 5 disjoint sets, where each set contains temporally adjacent time points.
We provide a table of validation MSEs averaged across the 5 folds in Table 3. Given our discretization of A and w values, we choose w = 0.18 and = 0.01 and keep these values fixed in our experiments presented in Appendix C. Figure 11 visualizes the corresponding model selected by 5-fold cross validation. The best estimated model is
(1) = 0.636G (1) + 1.255 (1)
(23)
(()) 0.839 (2) 0.833 (1)
(24)
20

---

Published in Transactions on Machine Learning Research (12/2025)
w
X
0.16
0.17
0.18
0.19
0.20
0.25
0.0075
0. 0.0538 0.0538 0.0539
0.0982
0.0963 0.1156
0.01
0.0538 0.0538 0.0538
0.0981
0.0963 0. 0912
0.02
0. 0.1285 0.1285 0.1287
0.1275 0.1256 0.4403
0.03
0. 1285 0.1285 0.3160
0.3166 0.3202 0.3451
0.15
0.0788
1
10.65
14.94
Table 3: Table of hyperrarameter combinations and corresponding validation MSEs for choosing the STRidge and w hyperparameters. Auy combination of hyperparameter values that led to a mull model (selecting none of the covariates to include in the model) for any of the validation folds is reported the hyperparameter combination that should give sparser models (larger À and larger w) in the case of tied validation MSE, which leads to choose À 0.01 and 0.18 (bolded and highlighted).
4_". We choose
Prey
Obs Prey
Obs. Predator
1.50
Abundance
1.25
0.50
0.25
Time
Figure 11: Lotka- Lotka-Volterra solutions (solid lines) for the equations selected using SINDy with hyperparameters 0.01 and w 0.18 using one observational dataset (points).
which selects the correct terms in each differential equation with low error on the fitted å coefficients.
E
Stability with varying numbers of bags: Lotka- example
The computational complexity of estimating the weights associated with each possible model m E M+ is not determined by the cardinality of the set M+, but rather the number of bags B used to estimated the weights associated with each model. Intuitively, we can generally expect that a larger number of bags (i.c., larger B) will correspond to a better estimate of the weights for each model m across selection methods. Figures 12 and 13 visualize the empirical CDF of instabilitios 6, across trials j E [N], where N 200 in the experiment. A separate line is plotted for each
regression experiments and N 100 in the choice of B.
In both plots, we see that B = 10 provides the least stability, evidenced by the larger proportion of 8ys near = The regression example provides evidence that generally, increasing the number of bags improves stability. However, this plot does not appear to clearly evidence this hypothesis as well as the Lotka- Volterra experiment. We hypothesize that since the covariates are highly correlated in the regression example with many reasonable models able accurately predict the response, this setting may yield more variable stability results across the number of bags compared to the Lotka- Volterra example, which has clear, unique model.
21

---

Published in Transactions on Machine Learning Research (12/2025)
Regression
argmaax-A
top-2=Ak.0
argmax An,
iPoo + Ar.
10-
0.4
0.0
B 10
B 100
0 1000
5000
B 10000
Figure 12: Empirical CDF of instability mcasures 6, across trials j E [N] in the regression experiment for varying number of bags across four different selection methods: the argmax, top-2, argmax008 ipo.3- The values of k e, and chosen to be consistent with the CDF plot in Figure 8.
Lotka Volterra
top-2+.Ak
argmax 0.02 A
Arc
00
0.2
0.0
0.0
08
1.0
0.0
0.2
0.8
1.0
0.0
0.2
0.4
1.0
0.0
0.2
0.6
0.
Figure 13: Empirical CDF of instability measures 83 across trials E [N] in the Lotka- Volterra experiment for varying number of bags B across four different selection methods: the argmax, lop-2, argmax0.0g ipo.63. The values of k E, and t were chosen to be consistent with the CDF plot in Figure 3.
F
Graph subset selection
F.1
Cross-validation to select
We select the penalization hyperparameter À in (11) via 5-fold validatii,.. The validation log- lig-liill averaged across the 5 folds computed as
wal-lo--ikeliihoo((:3333
1 Σ
(og θ,, -r() ~11@,11.).
(25)
where Ô, is computed by solving (11) using the training data in fold index v, and Σ, is the sample covariance matrix computed from the validation data fold index u. Computing (25) on the flow cytometry data across various choices of E [1,..., 500] yields the curve shown in Figure 14. Based on this figure, the best choice of is X 77, which we keep constant throughout experiments in $5.2.1.
22

---

Published in Transactions on Machine Learning Research (12/2025)
Log Log-Liken Likeiihoood (5 folds)
Average Validation
-66-
-68
-69
100
200
300
400
500
A
Figure 14: Average 5-fold validation log-likelihood across different choices of A. The that maximizes the log- likelihhod averaged over 5 held-out validation folds is À 77.
G K-means clustering
G.1 Data generation
We randomly generate N = 100 independent datasets in Rax2 with = 30 total data points per datasel and
k = 3 distinct clusters. To generate the three clusters, we generate n2 = 5 samples from N ([1.5,
5 samples from ([3.5, 3.5] 12), and n3 20 samples from (M2.5,, 2.5] 0.3² where N denotes a Gaussian distribution, and In x 2 identity matrix. An example one of our generated datasets is shown in Figure
Example Simulated Clustering Data
Cluster
Cluster
Cluster
1
2
5
X1
Figure 15: Example clustering dataset where the sample size is n = 30 and the number of clusters is x: = 3.
G.2 Base algorithm A: the "elbow" plot from A- x-means
A common technique in unsupervised clustering for selecting the number of clusters is via the "elbow" of the plot of the sum of squared distances (SSD) vs. к. We deterministically select as described in Algorithm 2, which selects the number of clusters by some hyperparameter for the x value that leads to the minimum observable slope the SSD vs. к plot. In equation 26, S, is the set of points in cluster i, and 14 is the cluster center for cluster selected via k-means. In our generated example, the maximum number of clusters
M 29, and the slope tolerance 5. Additionally, for each run of k-means, the clusters are randomly initialized.
23

---

Published in Transactions on Machine Learning Research (12/2025)
Algorithm 2 Computing the number of clusters input Maximum number of clusters M, dataset X, slope tolerance w for == I...., M do Compute the sum of squared distances (SSD):
SSD = Σ
Σ
||, - -0:11-
(26)
i 2yE8,
Compute the slope of SSD curve using the current SSD, and the previous SSD
SSD, SSD
slope,
(K
if slope, <w then output Number of clusters K end if end for output Number of clusters M
H Decision tree classification on single cell transcriptomics data
We investigate a model selection problem involving deciding the structure of a decision tree for multiclass classification. In this setting, we only measure two decision trees as distinct from each other if the trees split on the same features for given split. Therefore, we do not consider the exact split value for continuous features in deciding whether two selected models are the same.
We explore the decision tree stability with data from single cell transcriptomics to classify stem cell fate based gene expression. Mouse embryonic stem cells were sequenced for their gene expression and subsequently labeled as three different cell types: embryonic stem cells (ESC), epiblast-lib Mike intermediate cells (EPI), and neural progenitor cells (NPC) (Veleslavov and Stumpf, 2020). The cell fate label noted to be assigned by clustering, so there uncertainty in the true cell type label.
The dataset contains 547 cells and 96 gene expression values per cell. To explore stability in setting where the sample size is much lower than the feature dimensionality, which is common problem in genetics, we randomly subsample the number of cells to 50 and compute decision tree stabilities. As shown in Figure 16, gene expression is high correlation among genes, which may lead to instabilities in the fitted decision tree for small perturbations in the training data.
The base algorithm A corresponds to a multiclass classification tree. Broadly, decision trees create a sequence of splits on features that sort the training data samples into nodes, where a particular split is chosen minimize the impurity of the samples in the two split nodes. Greedy algorithms are typically implemented to sequentially choose the locally best split at each node. More information decision trees can be found in Breiman et al. (1984). Additionally, Xin et al. (2022) formulates algorithm for finding the full Rashomon set for decision trees, which can be used as an alternative base algorithm to rank decision tree structures.
For this classification problem, we choose to measure dissimilarities among samples in a node via Gini impurity, which can written as
Σ pra(1 -Pia)
(27)
for d E [D] classes at node I, where the proportion of class d in node I is
1 Σ T(y= d)
Pid
(28)
ni yEQ:
24

---

Published in Transactions on Machine Learning Research (12/2025)
Empirical Correlation Matrix
11
Correlation
-4.2
Figure 16: Empirical correlation matrix of gene expression across 50 cells from data in Veleslavov and Stumpf (2020). Rows and columns that show null value correspond to genes that showed variability in expression across the 50 cells.
25

---

Published in Transactions on Machine Learning Research (12/2025)
and Qi denotes the split of node 1 based on feature index and split value ty:
Qt {(xx y)1 24 t Q = Qi \ Q/
(29)
(30)
Our experiments use the default hyper-parameters in the sklearn. tree. www.ssionrreeeCassiiiee... function.
H.1 Results
Since we have one dataset of gene expression across cells, N = in this example. Table 4 shows the stability results and average number of models returned across TOO trials. Notably, ip, is absent from this table because this selection method is specific to variable selection, which a special class of model selection.
Based on Table 4, significant instabilities are seen for both the argmax A and argmax ÁK, model selection methods, while top-2 ÁK, provides enhanced stability, For the inflated argmax, we chose 0.03 that the stability matched that of the stability of top-2 in order for to directly compare the number of returned models. For fixed level of instability 8= 0.12, the inflated argmax returns a smaller number of models on average across LOO trials, showing that our method adaptively returns the smallest set of models based on the underlying uncertainty in the selection.
For = 0.04, we obtain 0.04 instability, corresponding to a high degree stability with a LOO average of 2.88 returned models. We visualize the top selected models across bootstrap samples Figure 17. From these trees, we can see that the first and second most common selected trees across bootstrap samples are completely disjoint, potentially indicating strong correlations among the features for both trees. However, the third most commonly selected tree differs from the top tree by only one split, suggesting there may be high correlation between the gene expressions of Lin28a and Nanog. Future investigations can focus on disentangling the associations among the genes selected in the lop tree versus the second most selected tree to elucidate their associations with the cell fates.
SELECTION METHOD LOO INSTABILITY AVG. LOO SET SIZE
argmax A argmax o AK,B top-2 AK.B (Ours) argmax0000
0.40
1.00
0.24
1.00
0.12
2.02
AK.B
0.12
1.74
Table 4: Table of stability and the average number of models returned across LOO trials for selecting the structure of the multi-class decision tree for a gene expression dataset from Veleslavov and Stumpf (2020). In this experiment, AA is a classification decision tree, and Áx, B bootstrapped classification decision tree with B 10, 000 and K 45. In computing the returned models, if there is an exact tie in the top models, all ties returned, which explains why top-2 2 returns slightly more than 2 models on average across the LOO fits.
26

---

Published in Transactions on Machine Learning Research (12/2025)
Top Decision Tree
Top 2 Decision Tree
Top 3 Decision Tree
Gdf3
Pou5f1
Gdf3
s
Lin28a
Zfp42
NPC
Klf4
Nanog
Zfp42
NPC
EPI
lgf2
ESC
Fgfr2
ESC
NPC
EPI
Igf2
ESC
$
$
EPI
NPC
EPI
NPC
EPI
NPC
Figure 17: Visualization of the decision tree structures for the top 2 trees chosen by bootstrap sampling the data for K = 45 samples per bag across B = 10, 000 bags. Differences between the top and second tree are indicated in blue in the second tree, and the differences between the first and third trees are indicated red. The names of the genes that determine the split at each node are shown in ovals, and the class labels NPC, ESC, and EPI are listed for each leaf. The top 1 decision
was selected by 4.74% of bags, the top 2 decision was selected by 2.80% of bags, and the third decision tree was selected by 2.24% of bags. The exact values used for the splits cach node were not considered in our stability results.
27


Reviews:
[{"review_id": "TFTJleOzQi", "reviewer": "Reviewer_TwV9", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this work, the authors introduce a new technique for model selection. The paper is clear and well written. I found the proposed technique new and the multi-class classification perspective is fresh and interesting. From there, using the newly introduced inflated argmax is an interesting choice. The technique only requires one interpretable hyperparameter. The experimental results provide reasonable evidence that the technique has practical merits.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The new perspective proposed in the paper is well justified and supported by experimental results. The idea is intuitive enough that its introduction and discussion would be of interested to the community.", "audience_yn": "Yes", "audience_explanation": "Although there is vast literature on the topic of model selection, this paper brings a new perspective. I think that TMLR readers and practitioners will find this paper interesting.", "requested_changes": "I have a few questions/comments (in no particular order):\n\nA- In section 2, model selection is described as choosing \"the model(s) [...] that provide the best fit to the data.\" I believe that this description is slightly misleading as fitting is not the only quantity involved in model selection. For example, in variable selection, selecting all variables will (in general) fit the data better than selecting a subset but this complete model is often not desirable. It is clear to me that the authors understand this difference, but a discussion around this topic might help practitioners reading the paper and trying to use their technique.\n\nB- Although one of the claims in the paper is that the proposed approach extends beyond variable selection, all experiments in the main manuscript only involve variable selection. There is one experiment involving decision trees in Appendix F (not mentioned in the main text). I would have appreciated one example (and preferably two) in the main text to support this claim. As a suggestion, selecting k in k-means is a classical example that the authors mention in the introduction. There are many techniques/measures for this algorithm (the silhouette score, the adjusted Rand index, adjusted mutual information, and information-theoretic approaches, for example). Although I understand that algorithmic-specific techniques may yield better results than a general technique like the one proposed here, it would be interesting to see the differences. Moreover, could these standard measures be used as algorithm A in the proposed framework?\n\nC- The clustering literature has analyzed the problem of the relationship between the selected number of clusters and the algorithmic stability of the optimization techniques under subsampling. See [Shai Ben-David, Dávid Pál & Hans Ulrich Simon. Stability of k-Means Clustering. COLT 2007], for example, where a negative association between both concepts if formally proven. The authors should discuss how these results relate to their approach.\n\nD- There is also ample theory on model selection for the LASSO itself. Although I understand that LASSO can be regarded as a form of variable selection technique, it is a model in itself. It commonly comes with a hyperparameter controlling the strength of the L1 regularization. In the paper, these type of parameters are described as \"optimization hyperparameters\" (in Appendix D.1.2, for example). However, their values affect directly the cardinality of the solution in non-trivial ways. Hence, their values cannot be optimized independently, especially in the conditions analyzed in the paper (correlated variables) as the ensuing choice will be (as recognized by the authors) unstable. Would it be possible to discretize the range of this parameter and fold it into the proposed framework? Evaluating this experiment would be very interesting.\n\nE- I understand the statistical merits of returning multiple models instead of just one when this selection is unstable. However, I would like to see a discussion of what it means for the practitioner. How to deal with multiple models? In certain areas, like decision trees, one can just treat them as a forest, for example. In other cases, the interpretation is more uncertain. For variable selection, what would happen if the variables selected by two returned models do not overlap (or overlap only slightly)?\n\nF- The main stability result links the fraction of samples to the stability to the selection. Has the set of \"instability inducing\" samples any interpretable meaning? Are these samples just outliers? What should the practitioner do with these samples (if anything at all)?\n\nG- What happens if there are no models at the required confidence level? Can this even occur? Would lowering the confidence level just solve the issue? For problems where the model category (e.g., LASSO) is a crude approximation of the underlying data is delta stable? In the sense that the proposed technique should not transition sharply from returning no models to returning \"all\" models as delta is relaxed.\n\nH- The proposed technique comes with a stability guarantee but Theorem 4.2 has no proof. If its proof follows immediately that of Theorem 17 by Soloff et al., please state this in the manuscript. If it does not follow immediately, a proof is needed. Without this clarification, the theory feels vacuous.\n\nI- What about uncountable set of models (with a continuous parameter, for example)? If this is out of scope, the authors should clearly identify this problem constraint where they state the main claims of the paper (abstract, Section 1.1, conclusions). If it is not out of scope, they should describe how to deal with such a scenario.\n\nJ- The statement in Section 6 \"We define as notion of stability [...] that computes the proportion of LOO selected models that are disjoint from the selection made from by training on the full data\" seems slightly misleading. Since bagging is used, no model is trained on the full data. Did I miss a critical step or should the wording be adjusted?\n\nMinor points:\n- Appendices D, E, and F are not mentioned in the main text.\n- There is extensive literature pertaining model selection. The overview of the topic in the paper seems to be constraint to a specific type of techniques (i.e., bagging). I believe that the authors should devote a couple of sentences to the bigger picture.\n- Paragraph after Equation (1) in the sentence \"this weight can be thought of as a probability of each candidate model,\" the authors should clarify which event this probability is related to.\n- In Section 2.1, I found the sentence \"We denote bagging a base algorithm A as A_{K,B}\" very unclear. Maybe something like \"We denote by A_{K,B} the algorithm resulting from applying a bagging technique on algorithm A\" would be clearer. There may be better formulations that the one I propose here.\n- In section 2.1, the sentence \"the ensemble weight w_m simply counts the fraction of bags wherein model m was selected\" is unclear and needs more elaboration. Maybe clarify that for simple model selection each vector w^b is one-hot would help.\n- In section 2.1, \"A_{K,B} provides a more stability\" -> \"A_{K,B} provides more stability\".\nAt the start of Section 2.2, the notation \"\\hat{M} \\in \\argmax_m \\hat{w}_m\" implies that \\hat{M} is a single element and not a set. Replacing \"\\in\" by \"\\subseteq\" may be a better choice.\n- I believe that should Definition 2.1 involve \"a given dataset D\" instead of \"all datasets\"? The measure s dataset specific.\n- In Section 4.1, \"denotd\" -> \"denote\".\n- In Definition 4.1, which norm is used in \"dist\"?", "broader_impact_concerns": "Model selection, although of critical importance in machine learning, can be clearly misused to infer all sorts of misleading results if applied without care. As such, the authors should probably add a paragraph explaining that a model that exhibits high stability does not necessarily imply that the \"ground truth\" model was found with high probability or anything remotely similar. Of course, this does not subtract value form the contribution, but such a general statement would be highly welcomed.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "pyduBgXPAV", "reviewer": "Reviewer_TSM4", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper is motivated by the observation that standard frequentist model selection procedures can be very unstable, potentially returning different results even when a single observation is dropped. The authors propose to stabilize black-box model selection algorithms by (1) using bagging to obtain weights for each model and (2) applying an \"inflated argmax\" procedure to select a set of top models that is stable across dropping observations. The authors show that this procedure provides to a precise form of \"model selection stability.\" Their theoretical result relates the degree of stability to algorithm's parameters, which aids users in selecting these parameters. Experimental result show that the proposed approach provides equal or greater stability compared to alternatives. \n\nStrengths:\n\n+ easy-to-use, black-box method\n+ works well empirically \n+ simple but precise theory that guides algorithm tuning \n\nWeaknesses:\n\n- leave-one-out stability is a fairly weak guarantee \n- unclear why the goal of selecting a stable \"top\" model is practically relevant", "claims_and_evidence_yn": "Yes", "claims_explanation": "**Claim 1:** *Our method selects a small collection of models that all fit the data*\n\nSupported by experimental results – 1 simulation study, 2 real data examples, all using quite different models \n\n**Claim 2:** *...and it is stable in that, with high probability, the removal of any training point will result in a collection of selected models that overlaps with the original collection.*\n\nSupported Theorem 4.2\n\n**Claim 3:** * the proposed method yields stable, compact, and accurate collections of selected models, outperforming a variety of benchmarks.*\n\nSupported by experimental results.", "audience_yn": "Yes", "audience_explanation": "Black box model selection is an important and common problem in ML. Users of ML methods want their procedures to be \"robust,\" including being stable to small changes in the data. So, I could see this method being further developed and/or used in practice.", "requested_changes": "I have two requested changes, both of which aim at clarifying the motivation and better justifying the specific goals of the work: \n\n1. **Justification of stability criterion:**  The leave-one-out stability criterion seems quite weak. For example, the theory from Huggins \\& Miller (2023) shows that BayesBag can be stable even when the full dataset is resampled. So, I think this choice needs better justification.\n\n2. **Goal justification:** I would like to see some justification for why choosing a single (or small set) of \"top\" models is a sensible goal. Indeed, once we effectively have a distribution over models, why not select the smallest set of models with total probability at least $1 - \\delta$? This seems like a more interpretable goal, as it aims to guarantee the \"true\" top model is getting returned. For example, if the models have weights (.35, .35, .25, .1), the inflated argmax might select models 1 and 2 but leave out model 3, even though the third model still has considerable weight. The 90% \"credible set,\" on the other hand, would return models 1–3.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "UYpdLfMOVJ", "reviewer": "Reviewer_X5xN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "# Summary\n\nThe authors consider the problem of stable model selection, i.e., model selection which gives consistent results even in the presence of small changes to the dataset. They define a novel notion of model selection stability which applies to methods which return a set of possible models, given a dataset. The notion measures the frequency with which the set of models selected by applying the selection algorithm to the entire dataset is disjoint from the set of models selected when one point is removed from the dataset. Inspired by methods from stable classification, they propose the inflated argmax for stable model selection. They prove that models selected by the inflated argmax are guaranteed to be stable with respect to their proposed notion of model stability. Finally, they validate their method with three experiments, two on synthetic data and one on real flow cytometry data.\n\n# Strengths\n\nThe fundamental problem considered (stable model selection) is interesting and useful. The shortcomings of existing selection procedures (e.g., when several models are near optimal, selecting the argmax will be unstable; if the selection probability of several models is close to the selection threshold, this will also be unstable) were clearly explained. Framing the problem of model selection in a multiclass classification framework in order to employ techniques from that literature is also an interesting connection.\n\n# Weaknesses\n\nFor major weaknesses, see the discussions below.\n\nTypos/notational clarifications:\n- In Definition 4.1, what is the norm being used in $||q-w||$? Is it the L1 norm of the probability distributions considered as vectors?\n- What is the value of $\\varepsilon$ in Fig. 1?\n- At the end of Section 5.1.1 it is stated that the base algorithm $\\mathcal{A}$ for the experiment is SINDy, while then end of Section 5.1.2 it is stated that the problem is solved using STRIdge penalization as the base algorithm.\n- The formatting of the headings of Table 2 should be fixed. There is no spacing between the headings, making it difficult to read.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The evaluation criteria used in the experiments are insufficient for showing the benefit of the proposed method, specifically on the real data.\n\nIn particular, considering only stability without a notion of correctness is not convincing, because any algorithm which outputs a single model regardless of the input dataset is perfectly stable while also outputting a prediction set of minimal size. So in the real dataset, such an algorithm would have obtained the best results if it was only evaluated according to the metrics reported in Table 2.\n\nFor the Lotka-Volterra model experiment, the results include a notion of correctness (utility-weighted accuracy, equation (7)) which also takes into account the number of models returned by a given method. This is an improvement over the real data experiment. Unfortunately, the results are not favorable for the proposed method. For the two strictest worst-case instability levels, the inflated argmax obtains 0% utility weighted accuracy (Fig. 3, middle plot). This seems especially problematic as it is stated that \"our primary interest is in small $\\delta$ (i.e., high stability)\" (pg. 7, after equation (6)).\n\n---\nPost rebuttal: The proposed method produces several model sets on the Pareto frontier in the Lotka-Volterra experiment. The discussion of utility-weighted accuracy was also helpful. In order for the flow cytometry results to be meaningful, however, some notion of correctness (even a qualitative discussion based on domain knowledge) is necessary.", "audience_yn": "Yes", "audience_explanation": "The problem of stable model selection is of interest to TMLR's audience, but I don't believe the paper has adequately motivated why the notion of stability that it proposes is a reasonable one. In particular, the definition of model selection stability (Definition 2.1) seems unfairly biased against simple model selection procedures. For instance, if the model selection procedure selects $\\hat{M}$ which contains 1000 models, and after removing a single datapoint, the resulting $\\hat{M}^{\\setminus i}$ differs from the original in 999 of the 1000 models but has a single model as overlap, this will not be penalized at all in the stability calculation, even though intuitively it seems that this is a very unstable selection procedure. A natural fix to this is to use the IoU of the two selection sets, i.e.,\n\n$$\\frac{1}{n}\\sum_{i=1}^n \\left(1 - \\frac{|\\hat{M}\\cap \\hat{M}^{\\setminus i}|}{|\\hat{M} \\cup \\hat{M}^{\\setminus i}|}\\right).$$\n\nWhen $\\mathcal{M}$ is a simple model selection procedure, this reduces to the formula in equation (2). In the case of the example above, such an example would still be highly penalized (it would add 1998/1999 to the sum).\n\nIn order to meet the TMLR criteria, I believe some justification should be provided for why the proposed stability notion is still reasonable in light of this example.\n\n---\nPost rebuttal: In light of the discussion and the other reviewers' interest, I have updated my score. I still believe a more thorough discussion of the limitations of the stability notion (implicitly favoring larger decision sets; the need for also considering utility) should be included to give readers a complete picture.", "requested_changes": "1. Please explain how the proposed stability notion can be considered useful in light of the example provided above.\n\n2. Please provide some objective measures of correctness on the real data experiment.\n\n3. How should the results showing 0% accuracy at high levels of stability (small $\\delta$) in Fig. 3 be interpreted?", "broader_impact_concerns": "N/A", "additional_comments": "There are some other recent works on model set selection which would be interesting to discuss. Specifically, there has been some recent interest in constructing Rashomon sets of models, which are defined as model sets which have performance close to that of the single best model on a given dataset. See, e.g.:\n\n>Xin, R., Zhong, C., Chen, Z., Takagi, T., Seltzer, M., & Rudin, C. (2022). Exploring the whole rashomon set of sparse decision trees. Advances in neural information processing systems, 35, 14071-14084.\n\nThis particular work focuses on sparse decision trees in particular, so the methods proposed in the submitted manuscript are more general. Nevertheless, it seems like there would be some strong connections between model selection stability considered in the submitted manuscript and sets of near-optimal models considered by these previous works. This connection should at least be mentioned and I believe it would strengthen the paper to discuss the connection (or dissimilarity) in some detail.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)