Exploring and Interacting with the Set of Good Sparse Generalized Additive Models
Authors: Chudi Zhong
Zhi Chen
Jiachang Liu
Margo Seltzer
Cynthia Rudin
Venue: NeurIPS 2023 poster (2023)
Keywords: Interpretability
human-model interaction
generalized additive model
Rashomon set

Abstract: In real applications, interaction between machine learning models and domain experts is critical; however, the classical machine learning paradigm that usually produces only a single model does not facilitate such interaction. Approximating and exploring the Rashomon set, i.e., the set of all near-optimal models, addresses this practical challenge by providing the user with a searchable space containing a diverse set of models from which domain experts can choose. We present algorithms to efficiently and accurately approximate the Rashomon set of sparse, generalized additive models with ellipsoids for fixed support sets and use these ellipsoids to approximate Rashomon sets for many different support sets. The approximated Rashomon set serves as a cornerstone to solve practical challenges such as (1) studying the variable importance for the model class; (2) finding models under user-specified constraints (monotonicity, direct editing); and (3) investigating sudden changes in the shape functions. Experiments demonstrate the fidelity of the approximated Rashomon set and its effectiveness in solving practical challenges.

Exploring and Interacting with the Set of Good Sparse Generalized Additive Models
Chudi Zhong
Zhi Chen'
Jiachang Liu'
Margo Seltzer?
Cynthia Rudin'
Duke University 2 The University of British Columbia {chudi zhong, zhi chen1, jiachang liueeduke edu mseltzerècs ubc.ua, ca, cynthia@cs duke. edu
Abstract
In real applications, interaction between machine learning models and domain experts is critical; however, the classical machine learning paradigm that usually produces only single model does not facilitate such interaction. Approximating and exploring the Rashomon set, i.c., the set of all near-optimal models, addresses this practical challenge by providing the user with searchable space containing a diverse set of models from which domain experts can choose. We present algorithms to efficiently and accurately approximate the Rashomon set of sparse, generalized additive models with ellipsoids for fixed support sets and use these ellipsoids to approximate Rashomon sets for many different support sets. The approximated Rashomon set serves as a cornerstone to solve practical challenges such as (1) studying the variable importance for the model class; (2) finding models under user-specified constraints (monotonicity, direct editing); and (3) investigating sudden changes in the shape functions. Experiments demonstrate the fidelity of the approximated Rashomon set and its effectiveness in solving practical challenges.
Introduction
A key ingredient for trust in machine learning models is interpretability; it is more difficult trust. a model whose computations we do not understand. However, building interpretable models is not easy; even if one creates model that sparse and monotonic in the right variables, it likely still not what a domain expert is looking for. In fact, domain experts often cannot fully articulate the constraints their problem requires. This leads to an often painful iterative process where a machine learning algorithm, which is typically designed to produce only one model, is now asked to produce another, and possibly many more after that. Hence, the classical machine learning paradigm broken, because it was not designed to facilitate choice among models. In other words, the classical machine learning paradigm was not designed for the real-world situation which domain experts see and understand a model, identify shortcomings in it, and change it interactively in specific way. There is instead an interaction bottleneck.
We work in a new paradigm where an algorithm produces many models from which to choose, instead of just one [.3). Algorithms in this new paradigm produce or approximate the Rashomon set of a given function class. The Rashomon set is the set of models that are approximately as good as the best model in the class. That is, it is the set of all good models. Rashomon sets for many real problems are surprisingly large [4 and there are theoretical reasons why we expect many simple, good models exist (3]. The question we ask here is how to explicitly find the Rashomon for the class of sparse generalized additive models (sparse GAMs).
'Equal Contribution
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

GAMs are one of the most widely used forms of interpretable predictive models [7- and have been applied to complex problems such as prediction of health outcomes from medical records 110 [1 where they can provide the same accuracy as the best black box models. GAMs can use sparsity regularization to generalize and to avoid constructing overly complicated models. Different regularization strategies have been used [8] 12 [3, and sparse GAMs have performance that is essentially identical to black box models on most tabular datasets [12]. Thus, an important previously-
missing piece of the Rashomon set paradigm is how to obtain the Rashomon set of sparse GAMs.
There are two major challenges: (1) Different from regression, the Rashomon set of GAMs in the classification setting has no analytical form; (2) The Rashomon set should lend itself easily to human- model interaction, mitigating the interaction bottleneck. To tackle these challenges, we find the maximum volume ellipsoid inscribed in the true Rashomon set. Our gradient-based optimization algorithm can find such an ellipsoid for a fixed support set. Leveraging the properties of ellipsoids, we can efficiently approximate Rashomon sets of numerous support sets that are subsets of the original support set. We show in Section 5 that this approximation captures most of the Rashomon set while including few models outside of it.
More importantly, using the maximum volume inscribed ellipsoid enables us to solve many practical challenges of GAMs through efficient convex optimization. We first use it to study the importance of variables among a set of performing models, called variable importance range. This is essentially the Model Class Reliance, which measures the variable importance of models in the Rashomon set [5]: : previously, ranges of variable importance have been studied only for linear models [4] and tree models (2 16. We also show how to efficiently use the ellipsoid approximation to search for near-optimal models with monotonicity constraints. The Rashomon set empowers users to be able to arbitrarily manipulate models. As users edit models, their edits are either already in the Rashomon set (which we can check easily) or they can easily be projected back into the Rashomon set using the technique in Section 4.4 Thus, being able to approximate the full Rashomon set for GAMs brings users entirely new functionality, going way beyond what one can do with just a set of diverse models (7-23]. Also, the ability to sample GAMs efficiently from the approximated Rashomon set allows us to investigate whether the variations or sudden changes observed in a single model are indicative of true patterns in the dataset or simply random fluctuations.
Since our work provides the first method for constructing Rashomon sets for sparse GAMs, there is no direct prior work; this is novel problem. There has been prior work on related problems, such as constructing Rashomon sets for decision trees [2] and the Rashomon sets for linear regression are simply described by ellipsoids [3; Rashomon sets for linear regression models have been used for decision making [24], robustness of estimation (25), and holistic understanding of variable importance
14 26
Our work allows these types of studies to generalize to GAMs.
2 Background
A sample is denoted as (x, y). where x is the p dimensional feature vector and y is the target. The th dimension of the feature vector is x... A generalized additive model (GAM) 7 has the form
g(Elv)) = wo (x1) + f2(22) +... fz(pp)
where wo is the intercept, ʃÿ's are the shape functions and is the link function, e.g., the identity function for regression the logistic function for classification. Each shape function Jj operates only on one feature Zj, and thus the shape function can directly be plotted. This makes GAMs interpretable since the entire model can be visualized via 2D graphs. In practice, continuous feature is usually divided into bins [8], thereby its shape function can be viewed as step function, i.c.,
B,-1
[5(,)) - Σ wyA 1.ba < xz S b3.0+1];
where (b are the bin edges of feature j. leading to By total bins. This is linear function on the binned dataset whose features are one-hot vectors. w = [W Wj11 , is the weight vector. For simplicity, we use this formulation in the rest of this paper.
Given a dataset D = ((x, y:) }/ we use the logistic loss L.(w, wn, D) as the classification loss. To regularize the shape function, also consider weighted &g loss (2( on the cocfficients and
2

---

total number of steps in the shape functions L,(w) as penalties. Specifically,
Bj-1
(w) =
Σ
K=
1by. Ei= Tay S bj,+1]] is the proportion of samples in bin k of feature j. The weighted & term not only penalizes the magnitude of the shape function, but also implicitly centers the shape function by setting the population mean to O: any other offsets would lead to suboptimal C2 (w). This approach is inspired by [8] which explicitly sets the population mean to O. To make the shape functions fluctuate less, penalize the total number of steps in all shape functions, i.e.,
where -= -
B
(o(w)=
Σ
1wzA + w3,+1]s
3=1
k=
which is similar to the lo penalty of [12]. Combining classification loss and penalty yields total loss:
[lw.wo, D)= Le(w, Wu. D) + XzC2(w) + (oC,(w).
(3)
Following the definition of 3 we define the Rashomon set of sparse GAMs as follows:
Definition 2.1. (0-Rashomon set) For binned dataset D with samples and m binary features,
where w E I defines a generalized additive model. The 8-Rashomon set is set of all w R with ((w, wo, D) at most 0:
R(0, D) == (w E Ih wo E R: wo, D) S 0}.
We use ((w,wo) represent C(w,wa, wo. D) and R(0) represent R(0, D) when D clearly defined.
3
Approximating the Rashomon Set
3.1 GAM with Fixed Support Set (Method 1)
Suppose merge all adjacent bins with the same wij. The set of bins after merging called the support set. For a GAM with fixed support set, the L. term is fixed and the problem is equivalent to logistic regression with weighted & penalty. For simplicity, in what follows, w also includes wg. In this case, the loss [(w) is a convex function, and the Rashomon set R(0) is convex set. Hence, we find the maximum volume inscribed ellipsoid centered at we approximate R(0). Specifically, the approximated Rashomon set is:
R = (w E RT
(w w.)
TO( - w.) s 1}.
(5)
where Q and w. are the parameters defining the ellipsoid that we optimize given different As.
Gradient-based optimization: Our goal is maximize the volume of î while guaranteeing that almost all points in it are also in R(0). Mathematically, we define our optimization objective as
min det Q) + C' -Ew (mxx([(w) - 0, O)].
(6)
Q.w.
The volume of R proportional to det Q) 11/2, therefore we minimize det Q) (22m)) in the objective. m is included in the exponent to normalize by the number of dimensions. In the second term of Eq (6), we draw samples from R( Q, w.) and add penalty if the sample is outside R(0); C' is large to ensure this happens rarely. We use gradient descent to optimize the proposed objective function.
Note that the "reparameterization trick" can be used to convert the sample w drawn from ft into purely stochastic part and deterministic part related to and wa, and calculate the gradients accordingly. Details are shown in Appendix To facilitate the optimization, we need a good starting point.
Since the support set is fixed, we treat L, as a constant. We initialize we using the empirical risk minimizer w" for L(w) and initialize Q by H/ (2(( L(w*))), ))), where H is the Hessian. The initialization is motivated by the Taylor expansion on L(w) the empirical risk minimizer w":
Initialization:
L(w) [(w") + VL (w w*) + 1 1/( - w H(w w²) O(w w*))).
3

---

Since w" is the empirical risk minimizer, the gradient VC = 0. If we ignore the higher order terms, we get quadratic function, i.c., [(w) L(w")= 5(4 ([w--))"H(>>>>>> ww)). Plugging this quadratic function into the Rashomon set formula, we get exactly an ellipsoid, i.e.,
ww).
"(w-""''(W--- S [(w").
Therefore, given this approximation, we initialize we with w" and Q with H/(2(0 - [(w)))) and then optimize through gradient descent.
3.2 GAMs with Different Support Sets (Method 2)
Method can accurately approximate the Rashomon set given a fixed support set. However, want the Rashomon set for many different support sets. applying Method repeatedly can be time-consuming. We next present a blocking method, Method 2, to efficiently approximate the Rashomon sets with many smaller support sets after getting approximated Rashomon set for large support set by Method 1. Approximating the Rashomon set with an ellipsoid for large support can be helpful in approximating the Rashomon sets with smaller support sets. Specifically, if the bins of the smaller support set are generated by merging bins in the original support set, is equivalent to forcing the coefficients of these merged bins to be the same. Formally, this is hyperplane P in the solution space defined by set of linear constraints
== (w E R wk = wh+1 = wk].
where the htth to the th bins in the original support set are merged to create a single bin in the new support set. The hyperplane's intersection with any ellipsoid, e.g., fi, is still an ellipsoid, which can be calculated analytically. Specifically, if we rewrite Q in the form of block matrices, i.e.,
Q 1.1--1
Q Jaks
Ql:k 1os+
Q
Q 1.
Q.+1m
Q
likk 1,1- 1am
Qx+1sm
we obtain the quadratic matrix Q for the new support set by summing up the 4 to t rows and columns, i.e.,
Q1:k- 1,1k-1
qisk-1
Q
Q
9i
Qmid
9k+1:m
11:k- 11:k-1xx+10 1,
qk+ 1ct1mm,
where gmid = Dimk L5=x Qiji 91:k- and qu+ Iim
column vectors obtained by summing all columns of Q Q Q1:k-1.6 let and Qe+I:mmkea The linear term (coefficients the term) of the original support is = Qwe. Similarly, we can obtain the linear term of the new support set by summing up the 1t to th elements the vector, i.e., e [fI,... ,......,,,, - Σx x li, fett,.... em, With this, we can calculate the new ellipsoid's center, We = Q and the upper bound of the new ellipsoid u 1 w.. Qw: + W5. Qãe.. This new ellipsoid Rp == P n R, which satisfies both Equations 5 and 7 is defined as
Rp = (w E R" (w (w) -. w.) S u),
where u O, Rr is empty. If want to merge multiple sets of bins, the calculations can be slightly modified to summing up sets of rows and columns separately corresponding to the sets of bins to be merged. Ry serves as a simple approximation of the Rashomon set for the smaller support set. Since essentially all of the original ellipsoid ft the Rashomon set, and Rp is subset of ît, every in Rp is also in the Rashomon set. In fact, since the support set gets smaller, L, becomes smaller too, and the loss even lower, i.c., for w in Rp, [(w) < 0 - Ao - k).
Explore different support sets with size K: Using the blocking method mentioned above, we can approximate the Rashomon set for many different support sets that are subsets of the original support set. The approximation of the ellipsoid for smaller support set more accurate when the size difference between the original support set and the new support set is small. Therefore, instead of using a very large support set at the beginning, we use a sparse GAM algorithm such as FastSparse [2] with relatively weak sparsity penalty to get support set S whose size K is moderately larger

---

than the size of support sets we want to explore, but covers all bins we want to potentially merge to form new support set. For simplicity of defining the loss threshold of the Rashomon set, here we consider only support sets with size K. If we want to explore different K, we can repeat this process. Suppose our goal is to explore the Rashomon set with loss threshold 5. The first step is to use the methods in Section 3.1 to obtain the ellipsoid ft approximating the 0-Rashomon set on the original support set, where I 5 + Ao K.) Then we can enumerate all possible ways ( (K=))
kc) ways
in total) to merge bins in the original support set to get any subset S whose size is K, and calculate the intersection (i.e., ellipsoid) Rp; Pg n R between R and the hyperplane P3 corresponding to merging S into S. All nonempty Rp are valid approximations of the Rashomon set.
4
Applications of the Rashomon Set
We now show four practical challenges related to GAMs that practitioners can now solve easily using the approximated Rashomon sets.
4.1 Variable Importance within the Model Class
Variable importance can be easily measured for GAMs. For example, Explainable Boosting Machines 27 measures the importance of each variable by the weighted mean absolute coefficient value
of bins corresponding this variable, i.e., Vx(x.,) = Σk=o TJ T This measure is based on one model, but given the Rashomon set, we can now provide more holistic view of how important a variable is by calculating the range of variable importance among many well-performing models; this is called the Model Class Reliance (14). Specifically, VI_ and VI, are the lower and upper bounds of this range, respectively. For example, the range of importance of feature j is defined as
[VI- (x.,). g)] min VI(x.j), max VI(x WER(
R()(
A feature with large VI_ is important in all well-performing models; a feature with small VI+ is unimportant every well-performing model. We use R to approximate R(0) and estimate VT_ by solving a linear programming problem with a quadratic constraint, and we estimate V VI+ by solving mixed-integer programming problem. Solving these problems gives a comprehensive view of variable importance for the problem, not just a single model.
Lower bound of variable importance: VI_ of feature j can be obtained by solving the following linear programming problem.
B, Σ
min [otg 0...0..........
ky s.t. (w w.) we) 1
(9)
Since Ty.k 2 O. Problem can be solved using LP solver. The objective minimizes variable importance, and the constraint ensures the solution is in the Rashomon set.
Upper bound of variable importance: The maximization problem cannot be solved through linear programming since [uJj.0: ...: -1] can be arbitrarily large. Instead, we formulate as a mixed- integer program. Let MM be relatively large number (e.g.. 200 is large enough for real applications (see Appendix Eb) and let J be binary variable.
By-1
max
Wyo y.
k=
s.t. (w w.) TQ( wo) S
(10)
Wjik + M
w9, Wiik + M (1 1x) wh Vk E (o,..., B,. 1}
W5,k S w/y,: -Wj, S w/S 4k E (o,...,B--11 13.
The last two constraints ensure that w/, defined as the absolute value of Wj,k- However, solving -ixeedintteger programs is usually cine-consuming. We propose another way to solve the maxi- mization problem. Since wj.A can be either positive or negative, we enumerate all positive-negative combinations for wj.k, & E {0,.... B3 1}. solve the IP problem with the sign constraint enforced, and choose the maximum value (see Algorithm 2
in the Appendix E
5

---

4.2 Monotonic Constraints
use-case for the GAM Rashomon set finding an accurate model that satisfies monotonicity constraints. This can be obtained by solving a quadratic programming problem. For example, if user wants fj(x;) to be monotonically increasing, the optimization problem can be formalized as:
min(w W - wee) s.t. Wj.k S Wj,++1, (..... By 1/,
(11)
where we check that the solution I 1; i.e., that solution in the approximated Rashomon set.
4.3 Find robustly occurring sudden jumps in the shape functions
Sudden changes (e.g., a jump spike) in the shape functions of GAMs known to be useful signals for knowledge discovery [T0] [1] and debugging the dataset (28]. However, due to the existence of many almost equally good models that have different shape functions, it is hard to identify if change is true pattern in the dataset or just random artifact of model fitting. However, with the approximated Rashomon set, we can calculate the proportion of models that have such change among all near-optimal models and visualize the shape functions (see Section 5.4).
4.4 Find the model within the Rashomon set closest to the shape function requested by users
As discussed, the main benefit the Rashomon set is to provide users a choice between equally-gowe models. After seeing GAM's shape functions, the user might want to make edits directly, but their edits can produce a model outside of the Rashomon set. We thus provide a formulation that projects back into the Rashomon set, producing model within the Rashomon set that follows the users' preferences as closely as possible. We find this model by solving a quadratic programming problem with quadratic constraint. Let Wrey be the coefficient vector that the requests. The problem can be formulated as follows:
s.t. (w w.)² Q(w w.) S 1. min ||ا Wrog
(12)
Problem 12 can be easily solved through quadratic programming solver.
5 Experiments
Our evaluation answers the following: 1. How well does method approximate the Rashomon set? How does the approximated Rashomon set help us understand the range of variable importance in real datasets? 3. How well does the approximated Rashomon set deliver model that satisfies users' requirements? 4. How does the Rashomon set help us investigate changes in the shape function?
We use four datasets: recidivism dataset (COMPAS) [29], the Fair Isaac (FICO) credit risk dataset 30 a Diabetes dataset (31, and an ICU dataset MIMIC-II
5.1
Precision and volume of the approximated Rashomon set
As we mentioned in Section 3 our goal is create an approximation of the Rashomon whose volume is as large as possible, while ensuring that most points in the approximated Rashomon set are within the true Rashomon set. To measure precision of our approximation. models within and outside the true Rashomon set R(0) are considered as "positive" and "negative" respectively; models within and outside the approximated Rashomon set /t are considered as "predicted positiveeneegatiee" Accordingly, the precision defined as the proportion of points (w, wo). equally weighted and within îg, having L(w,wn) S 0, i.e., they are in R(0). The recall is defined as the proportion of points in Voll (
R(0) that are also in R. Note that recall = (8()) = precision x Vai[ Since Vol(R(()) is
Vol(
a constant, given the same volume Vol( RR), recall is proportional to precision.
We will develop some natural baselines to compare to. For fair comparison, we rescale the Rashomon set approximated by baselines to have the same volume our ît. Then, we sample 10,000 points from our R and the rescaled baseline Rashomon sets, calculate their loss, and estimate precision.

---

diabetes
mimic2
compas
fico
preciaion
precision
precision
04
0.-
0.0
0.005
0.01
0.05
0.1
0.005
0,05
0.1
0.005
0.01
0.05
8.01
(0- C"M
(0- WC
(0 6" VV2
ebm
sphere
hessian
ours
(a) Precision of the approximated Rashomon sets as function of the Rashomon threshold 0. The larger the
value (0 - L*)/ L", the larger the Rashomon threshold. The precision obtained by our method dominates the baselines, indicating that, when the sizes of the approximated Rashomon sets are similar between methods, the Rashomon set approximated by our method is better inscribed in the true Rashomon set. (A, = 0.001, X2 0.001). Hessian is the starting point for optimization.
diabetes
mimic2
compas
fico
1.0
trecsionn
precision
00
precision
0.4
0.2
0.2
0.0
0.0
0.0
0.90
0.99
10
1.01
1.02
0.90
0.99
100
101
1.02
1.01
1.02
0.99
1.0
1.03
3.52
scaling ratio
scoling ratio
scaling ratio
scaling ratio
ebm
sphere
hessian
ours
(b) Tradeoff' between precision and the size of the approximated Rashomon sets (scaling factor to the power m is volume, which is proportional to recall). Our method dominates baselines on all four datasets. As the scaling ratio increases, the precision starts to decrease. (A, 0.001, x2 = 0.001, 0 1")1 Hessian is our starting point.
Figure 1: Precision versus the size of the approximated Rashomon set.
denoted as "hessian" One baseline is using our Hessian initialization of Q i.c., 1 in comparison. Another baseline is to use the identity matrix to approximate Q. i.c., the approximated set a sphere. Another type baseline comes from fitting GAMs on bootstrap sampled subsets of data. For fixed support set, calculate coefficient vectors from two different methods (logistic regression and EBMs [8) on many subsets of data. This approach samples only a finite set of coefficient vectors rather than producing closed infinite set such therefore its volume is not measurable. Thus, we use optimization to find the minimum volume ellipsoid that covers most coefficient vectors. This is the MVEE problem 33;; more details provided in Appendix [ Sampling from the posterior distribution of a Bayesian model is also a way get many different GAMs. We tried Bayesian logistic regression (using Hamiltonian Monte Carlo), but it is too slow (26 hours) to converge and produce enough samples to construct the minimum volume ellipsoid. So we did not include this baseline.
ellipsoid, and
Figure I compares the precision and volume of our Rashomon set with baselines. Figure la shows that when the volume of the approximated Rashomon sets the same between methods, our Rashomon set has better precision than the baselines. In other words, the Rashomon set approximated by our methods has the largest intersection with the true Rashomon set. The results for the hessian is worse than ours but better than other baselines, which means our proposed initialization is already better than other baselines. Also, as 0 becomes larger, the hessian method becomes better and sometimes close to the result after optimization. Logistic regression and EBM with bootstrapping do not achieve good precision, because coefficient vectors trained by these methods on many subsets of data can be too similar to each other, leading to ellipsoids with defective shapes that do not match the true Rashomon set boundary even after rescaling.
In general, shrinking the approximated Rashomon set leads to higher precision and expanding leads to better recall. Figure [b] shows the performance of different methods under this tradeoff. Our method still dominates others even when we rescaled the volume with different factors.
The previous experiments evaluate the performance when the support set is fixed. To explore many different support sets, we use the blocking method (Method 2) from Section 3.2 Method 2 is much faster than Method 1, we need to check that it performs similarly, according to
7

---

compas
compas
fico
fico
1.6
volume ratio
precision ratio
ratio
1.4
1.4
1.2
precision
1.2
1.1
1.2
1.0
8
og
of000
1.0
1.0
0.9
0.8
15
13
26
24
K
Figure Box plots of the volume ratio and precision ratio of Method over Method across different K, for two datasets. shows the precision and volume ratios are close to 1, indicating that Method performs similarly to Method 1.
both the volume and precision of the approximated Rashomon set. The volume ratio, defined by V (ometaett))a)))) VVl Rmethoaz), should be if the two methods perform similarly. The precision ratio, the precision of Method over the precision Method 2, again should close to the methods perform similarly.
Figure[showe the precision ratio and volume ratio 100 different new support sets where u i.c., the approximated Rashomon set is non-empty. Clearly, Method 2 will break down if we choose the support size K to be too small. The results from Figure indicate what we expected: precision and volume ratios are close to 1, even for relatively small K. Importantly, Method 2 is much more efficient than Method For example, the running time <0.001 seconds for Method but ranges from 300 to 100 seconds for Method on all three datasets (COMPAS, FICO, MIMIC-II). More results are Appendix
diabetes
shape function
shape function
DiabeesPPdigreeuunction
dw. dix)
(hot flx)
Insulin
predrcted po
SkinThickness
BloodPressure
Glucose
Pregnancies
00
3.0
2.0
variable importance range
not fix other coefs
fix other coefs
zio
(a) Variable importance range of the Diahetes dataset.
(b) Shape functions of "Glucose."
Figure 3: (a) Variable importance range of the Diabetes dataset and (b) shape functions of "Glucose." "
"Not fix": the shape function is obtained by solving Problem (IO). "Fix" coefficients of all other
features except "Glucose" are set to the same values as in w.. ()a 0.001,,s = 0.001, 1.012*)
5.2 Variable importance range
We solve Problem 9 and 10) to bound the variable importance; results on the Diabetes dataset shown by dark blue segments in Figure The light blue segments show the range with an additional constraint that coefficients of features other than the one we are interested in are fixed. Adding this extra constraint leads to lower VI+ but has negligible impacts on VI_. The "Glucose" feature has dominant VI_ and VI+ compared with other features, indicating that for all well- performing GAMs, this feature is the most important. In fact, the large VI_ indicates that Glucose remains important Rashomon set. This makes sense, because an oral
even for the model that relies least on it in glucose tolerance test indicates diabetes if the plasma glucose concentration is greater than two hours into the test. Features "BMI" and "Age" also have high VT_, indicating both of them are more important than other features.

---

1.-
the user ueg requires
1.5
requires
after adding monotonic constraint
predicted logit
1.-
predicted logit
1.0-
predicted logit
after editing
0.5
0.5
before editing
0.0
0.0
after cditing
before editing
-0.5
editing
-0.5
-05
100
120
20
100
120
BloodPressure
Pressuree
BloodPressure
(a) Monotonically decreasing
(b) User request case 1
(c) User request case 2
Figure 4: Example shape functions of "BloodPressure" that satisfy users requirement. (a) The shape function of "BloodPressure" is desired to be monotonically decreasing. But has a jump at BloodPressure ~60. The optimization time to find the monotonically decreasing shape function is 0.04 seconds. (b) Removes the jump by connecting to the right step. The optimization problem solved in 0.0024 seconds. (c) Removes the jump by connecting to the left step. The time used to get the model is 0.0022 seconds. (Ao = 0.001, 22 = 0.001, 0 C*)
Shape Function (Pirato)
Shape Function (PPretio)
Shope Funct on iPiretio)
Shape Function (Piratio)
Me
Me Aammmple
u.somple
A.somple
preteted ligh
j_samppl
100
FFraboo
PPratic
Frooto
PRrttic
Figure 5: Different shape functions (in red) of "PFratio" in MIMIC-II sampled from R.
Figure 3 shows shape functions of we and the coefficient vectors obtained when the lowest and highest variable importance are achieved for "Glucose." " All these shape functions follow the same monost-toootooicc trend, which means the higher "Glucose" level, the more likely it is that a patient has diabetes. The variable tends to have higher variable importance when it has a larger range in the deviate more
predicted logit. Shape functions based on WV1, without fixing other coefficients from wii when we do not fix other coefficients, there is more flexibility in the shape functions.
5.3 User preferred shape functions
A significant advantage of the Rashomon set is that it provides the opportunity for users to contribute their domain knowledge flexibly without retraining the model. Note that all models within the approximated Rashomon set predict equally well out of sample (see Table Appendix), can choose any of them. We use the "BloodPressure" feature in the Diabetes dataset as an example. Figure 4 shows the shape function of "BloodPressurre" (in gray), in which jump occurs when blood pressure is around 60. A potential user request might be to make the shape function monotonically decreasing. By solving Problem (TT) with the constraint
coefficients related "BloodPressure," we get the shape function colored in yellow in Figure The jump is reduced by dragging up the left step. Suppose user prefers to remove the jump by elevating the left step (see inset plot in Figure (4)). By solving Problem (2). we find that the specified shape function is not within the Rashomon set, and we find the closest solution in green, which still has a small jump 60. Another user might prefer to remove the jump by forcing both left and right steps at the same coefficients (see inset plot in Figure 4t) This specified shape function is also outside the Rashomon set, and the closest solution keeps the most steps as required with a small step up. Different user-specified shape functions lead to different solutions. The approximated Rashomon set can thus provide a computationally efficient way for users to find models that agree with both their expectations and the data. Our work could be integrated with tools such as GAMChanger [3 which allows users to directly edit a GAM, but there is no guarantee that the resulting model still performs well.
9

---

5.4 Changes in shape functions
When we see a jump the shape function, we might want to know whether this jump often exists in other well-performing models. Using the approximated Rashomon set, we can give answer. For example, the shape function of "PFratio" in the MIMIC-II dataset has a sudden dip around 330 (see the gray curve in Figure 55. The jump is deep that might wonder if such a jump is present in other well perfooeeron GAMs. We sample 10,000 points from our R and find that 7012 samples have a downward jump at the same position. We are surprised to find that in some of the remaining 2988 samples, the jump is instead upward. Figure Shows four different shape functions of "PFratio" sampled from R. The dramatic magnitude change in either direction could be caused by the small number of data points falling in this bin (Tin = 1.2c-4). Since the weight is so small, even a tremendous change leads to only small impact on the loss. In practice, using our proposed optimization methods in Section 4 users can reduce or amplify the jump in the shape function or even create a different curve, as in Figure 4 Note that our ability to find diverse set of shape functions illustrates an advantage of the Rashomon set that cannot be easily seen with other approaches.
6 Conclusion and Limitations
Our work approximates the Rashomon set of GAMs in user-friendly way. It enables users to explore, visualize, modify, and gain insight from GAM shape functions. Our work represents a paradigm shift in machine learning, enabling unprecedented flexibility for domain experts to interact with models without compromising performance. Our contributions open a door for enhancing human-model interaction through using the new toolkit provided by the Rashomon set.
One limitation that could be creatively considered in future work is that there are numerous possible support sets, leading to the question of how human might comprehend them. A carefully-designed interactive visual display might be considered for this task.
Acknowledgements
We thank the anonymous reviewers for their suggestions and insightful questions. We gratefully acknowledge support from grants DOE DE-SC0023194, NIH/NIDA R0I DA054994, NSF IIS- 2130250, NSF DGE-2022040,, and DOE DE-SC0021358 We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). Nous remercions le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG) de son soutien.
Code Availability
Implementations of our methods is available at http:///tthcccccoccccccconoo/oae
References
[1] Lco Breiman. Statistical modeling: The two cultures (with comments and a rejoinder by the author). Statistical Science, 16(3):199-231) 2001.
[2] Rui Xin, Chudi Zhong, Zhi Chen, Takuya Takagi. Margo Seltzer, and Cynthia Rudin. Exploring the whole rashomon set of sparse decision trees. Advances in Neural Information Processing Systems, 35:14071-14004 2022.
[3] Lesia Semenova, Cynthia Rudin, and Ronald Parr. On the existence of simpler machine learning models. In ACM Conference Fairness, Accountabiliry, and Transparency (ACM FAccT), 2022.
[4] Alexander Amour, Katherine Heller, Dan Moldovan, Ben Adlam, Babak Alipanahi, Alex Beutel, Christina Chen, Jonathan Deaton, Jacob Eisenstein, Matthew D Hoffman, et al. Un- derspecification presents challenges for credibility in modern machine learning. Journal of Machine Learning Research, 2020.
10

---

[5] Charles Marx, Flavio Calmon, and Berk Ustun. Predictive multiplicity in classification. In International Conference on Machine Learning. pages 6765-6774, 2020.
[6] Hsiang Hsu and Flavio Calmon. Rashomon capacity: A metric for predictive multiplicity in classification. Advances in Neural Information Processing Systems,
2022.
[7] TJ Hastic and RJ Tibshirani. Generalized Additive Models, volume 43. CRC Press, 1990.
[8] Yin Lou, Rich Caruana, Johannes Gehrke, and Giles Hooker. Accurate intelligible models with pairwise interactions. Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 123-631,22 2013.
[9] Rishabh Agarwal, Levi Melnick, Nicholas Frosst, Xuezhou Zhang, Ben Lengerich, Rich Caruana, and Geoffrey Hinton. Neural additive models: Interpretable machine learning with neural nets. Advances in Neural Information Processing Systems, 34:4699-47 2021.
[10] Rich Caruana, Yin Lou, Johannes Gehrke, Paul Koch, Mare Sturm, and Noemie Elhadad. Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission. In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1721-1730, 2015.
[11] Benjamin Lengerich, Mark Nunnally, Yin Aphinyanaphongs. Caleb Ellington, and Rich Caruana. Automated interpretable discovery of heterogeneous treatment effectiveness: A covid-19 case study. Journal of Biomedical Informatics, 130: 104086, 2022.
[12] Jiachang Liu, Chudi Zhong, Margo Seltzer, and Cynthia Rudin. Fast sparse classification for generalized linear and additive models. In Proceedings of Artificial Intelligence and Statistics (AISTATS), 2022.
[13] Asad Haris, Noah Simon, and Ali Shojaie. Generalized sparse additive models. Journal of machine learning research, 23(70), 2022.
[14] Aaron Fisher, Cynthia Rudin, and Francesca Dominici. All models are wrong, but many are useful: Learning a variable's importance by studying an entire class of prediction models simultaneously. Journal of Machine Learning Research, 20(177): 2019.
[15] Gabriel Laberge, Yann Pequignot, Foutse Khomh, Mario Marchand, and Alexandre Math- ieu. Partial order: Finding consensus among uncertain feature attributions. arXiv preprint arXiv:2110.135 13369, 2021.
[16) Gavin Smith, Roberto Mansilla, and James Goulding. Model class reliance for random forests. Advances in Neural Information Processing Systems, 33:22305- 2020.
[17] Emilie Danna, Mary Fenelon, Zonghao Gu, and Roland Wunderling. Generating multiple solutions for mixed integer programming problems. In International Conference Integer Programming and Combinatorial Optimizution, pages 280-294. Springer, 2007.
[18] Izuwa Ahanor, Hugh Medal, and Andrew C Trapp. Diversitree: Computing diverse sets near- optimal solutions to mixed-integer optimization problems. arXiv preprint arXiv:2204.0382. 2022.
[19] Satoshi Hara and Takanori Machara. Enumerate lasso solutions for feature selection. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 31, 2017.
[20] Satoshi Hara and Masakazu Ishihata. Approximate and exact enumeration of rule models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018
[21] Salvatore Ruggieri. Enumerating distinct decision trees. In International Conference on Machine Learning, pages 2960-2968. PMLR, 2017.
[22] Salvatore Ruggieri. Complete search for feature selection in decision trees. Journal of Machine Learning Research, 20: 104-1, 2019.
[23] Nicholas Kissel and Lueas Mentch. Forward stability and model path selection. arXiv preprint axX1v:2103033462, 2021.

---

[24] Theja Tulabandhula and Cynthia Rudin. Robust optimization using machine learning for uncertainty sets. In International Symposium on Artificial Intelligence and Mathematics, ISAIM 2014, Fort Lauderdale, FI, USA, January 6-8, 2014, 2014.
[25] Beau Coker, Cynthia Rudin, and Gary King. A theory of statistical inference for ensuring the robustness of scientific results. Management Science, 67(10))6174-1193 6197, 2021.
[26] Jiayun Dong and Cynthia Rudin. Exploring the cloud of variable importance for the set of all good models. Nature Machine Intelligence, 2(12):810 122-84442 2020.
[27] Harsha Nori, Samuel Jenkins, Paul Koch, and Rich Caruana. Interpretml: A unified framework for machine learning interpretability. arXiv preprint arXiv:I909.00
2019.
[28] Zhi Chen, Sarah Tan, Harsha Nori, Kori Inkpen, Yin Lou, and Rich Caruana. Using explainable boosting machines (ebms) to detect common flaws data. Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 534-551. Springer, 2021.
[29] J. Larson, S. Mattu, L. Kirchner, and J. Angwin, How analyzed the COMPAS recidivism algorithm. ProPublica, 2016.
[30] FICO, Google, Imperial College London, MIT, University of Oxford, UC Irvine, and UC Berkeley. Explainable Machine Learning Challenge. https:// https:://comunitty htt////ommum.tcc.mmmm wft.initiyfffo....ooom.ooo fico. com/s/ // 2018.
explainable machime- hww.siiabbmmmmciie-eeeee
[31] Jack W Smith, James E Everhart, WC Dickson, William C Knowler, and Robert Scott Johannes. Using the adap learning algorithm forecast the onset of diabetes mellitus. In Proceedings of the annual symposium on computer application in medical care, page 261. American Medical Informatics Association, 1988.
[32] Mohammed Saeed, Christine Lieu, Greg Raber, and Roger G Mark. Mimic ii: a massive temporal icu patient database to support research in intelligent patient monitoring. In Computers in Cardiology, pages 641 641-644. IEEE, 2002.
[33] Michael J Todd and E Alper Yıldırı. On khachiyan's algorithm for the computation of minimum-volume enclosing ellipsoids. Discrete Applied Mathematics, 155(13)11311-174, 2007.
[34] Zijie J Wang, Alex Kale, Harsha Nori, Peter Stella, Mark Nunnally, Duen Horng Chau, Mi- haela Vorvoreanu, Jennifer Wortman Vaughan, and Rich Caruana. Gam changer: Editing generalized additive models with interactive visualization. Advances in Neural Information Processing Systems, Bridging the Gap: From Machine Learning Research Clinical Practice (Research2Clinics) Workshop. 2021.
[35] Steven Diamond and Stephen Boyd CVXPY: A Python-embeddded modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1-5, 2016.
[36] Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd. A rewriting system for convex optimization problems. Journal of Control and Decision, 5(1):42-60 60, 2018.
12

---

A Summary of datasets
We show experimental results four datasets: recidivism dataset (COMPAS) 29 the Fair Isaac (FICO) credit risk dataset (30 used for the Explainable MI Challenge, the Diabetes dataset [31] and an ICU dataset MIMIC-II
32 [ summarizes all the datasets. We note that these real worldw
datasets using features that are known to be used in systems of this type. Specifically, the COMPAS dataset is criminal recidivism dataset and we use features that are based on criminal history, age, etc. For FICO, the data were provided by FICO itself based on how credit scores are usually constructed. MIMIC-II is based on features observed in ICU patients that are predictive of death, and the Diabetes dataset also uses typical measurements taken during pregnancy that are indicative of diabetes. The features for each dataset are generally available for the decision and they are interpretable.
Dataset
Samples
Features
Classification task
COMPAS
6907
Predict if someone will be arrested I years of release
FICO
10459
23
Predict if someone will default on loan
Diabetes
768
9
Predict whether a pregnant woman has diabetes
MIMIC-II
24508
17
Predict whether a patient dies in the ICU
Table I: Summary of datasets
B
Sampling Uniformly from Ellipsoid
Algorithm 1 SampleFromElliisooid Q we)
Input: parameters of the ellipsoid Q E m2 we 1
Output: a point inside the ellipsoid w E 12m
I: u (0, I) // sample an m dimensional vector from standard multivariate Gaussian distribu-
tion
2: u u/ u/||u||z // normalize
get a unit-vector
3: r U), 1) get sample from uniform distribution
r: < /1/m // rescale to get the radius of sample in a unit sphere
5: y f- ru / is a random point inside a unit sphere
6: A, V = Eig( Q) // digen-decomposition, diagonal of A are the eigenvalues, and columns of V
are eigenvectors
x 4- VT // get point the ellipsoid x" Qx $
8: w x + w. // shift it so that the center is We
return w
Algorithm D describes the algorithm uniformly sample from the ellipsoid {w R" (w " w.) < 1). The algorithm first samples a random point inside a high dimensional unit sphere (line 1-5), and applies linear transformation (calculated from and w.) to get the point in the target ellipsoid (line 6-8). The whole process can be decomposed into a purely stochastic part, ie., sample in the unit sphere, and a deterministic part, which differentiable. Using this sampling algorithm, we can get samples for objective (6) and use gradient-basee methods to optimize Q and w. In addition, the algorithm is also used sample data for the problem in Section 4.3 and to calculate precisions in Section 5.I
C Precision and volume of the approximated Rashomon set
To run our method, set the learning rate to 0.0001 and run 1000 iterations. C' is set to 500. For the logistic regression and EBM baselines, we sample 2000 coefficient vectors by fitting GAMs on the bootstrap sampled subsets of data. We run logistic regression with lz penalty and EBM with interaction terms. We then find the minimum volume ellipsoid that can cover most coefficient vectors. Given set of coefficient vectors Wkamples, we solve the following problem:
2000
min Q,
det Q
+
max(II 1/2 ( - -w
10,-1000 1,0).
(13)
2000
13

---

We solve this problem via gradient descent. Ç is set to 1000. The number of iterations and learning rate in GD are set to 1000 and 0.01, respectively. We initialize Q by the ZCA whitening matrix and we by the average of Wamples-
After rescaling the Rashomon set approximated by baselines, we sample 10,000 points from our A and rescaled baseline Rashomon sets, calculate the loss, and get the precision. We include more figures in this appendix that are similar to Figure to compare with basclines using different values of 0 and 0.
2
0.0005
diabetes
mimic2
compas
fico
1.0
1.0
grecision
precision
preccion
0.8
0.2
0.0
0.0
0.0
0.005
0.1
0.05
0.1
0.005
0.01
0.05
0:005
0.01
0.05
0.1
0.003
0.02
(0-c'
10-
C'Ne"
10- C'Ne"
(8 19-6 2' Me2
2g3 = 0.001
diabetes
mimic2
compas
fico
1.0
3.0
precsion
0.8
grechloon
precision
0.6
0.4
3.0
00--
0.0
0.005
0.01
0.05
0.1
0.005
0.01
0.05
0.015
0.01
0.65
0.1
VVV'"
(0-
e'Me²
(0-
2"
1/2
(8 18-68 c" Me"
A, = 0,0015
diabetes
mimic2
compas
fico
1.0
1.0
oraccion
0.8
precision
precision
0.8
0.0
0.
0.6
0.6
0.
0.0
0.055
0.01
0.05
0.1
00005
0.01
0.05
0.1
0.005
0.05
0.025
0.02
0.00
01
(0-
e'Me'
(0
VC
(0-
VC
ebm
sphere
hessian
ours
Figure 6: Precision of the approximated Rashomon sets as function of 0. Our method always dominates other baselines. hessian is our starting point.
Figure&/compare coompares the precision of our method and baselines when the volume is fixed. The Rashomon set approximated by our method has the largest intersection with the true Rashomon set. This pattern consistent across all datasets and values of X, The Rashomon set approximated by the Hessian has lower precision but is always better than the other baselines. As f becomes larger, the Hessian method becomes better and sometimes comes close to the result after optimization.
shows the tradeoff' between the size and precision of the approximated Rashomon set for each method. The Rashomon set approximated by our optimization method is better than baselines given different values of 0.
Figure 7
We also report the optimization time of our method and baselines in Table 2 In most cases, proposed method has a run time slightly longer than logistic regression with bootstrapping but shorter than the EBM baseline. Baselines "hessian" and "sphere" do not require the optimization step, so they finish instantaneously. For this table, we ran the gradient descent on a CPU, whereas had we used GPU, it would be at least 10x faster.
Though we show that our method can find the ellipsoid with high precision, one may still wonder how well the ellipsoid approximates the true Rashomon set using other metrics. To answer this question, we show the scaling ratio that is needed to ensure points sampled from the surface of the ellipsoid are outside the true Rashomon set in Table 3 On average, each dimension needs to scale by only
14

---

1.005.C*
diabetes
mimic2
compas
fico
grecsion
0.6
preciiion
0.6
0.2
proctsion
0.6
0.8
0.-
0.0
100
0.99
1.02
0.8
0.99
3.0
1.1
102
scaling ratio
scaling ratio
Scaling ratio
0 = 1.01C"
diabetes
mimic2
compes
fico
1.0
3.0
preccon
procision
0.4
0.4
0.4
0.2
0.2
0.0
0.90
2.01
1.02
0.50
0.30
1.0
1.01
1.02
scalling
ing
Scaling retto
ratio
1.055*
diabetes
1.0
mimic2
1.0
compas
3.5
fico
predsion
0.5
prectsion
0.0
grecision
precssion
0.5
0.5
0.0
0.8
scaling ratio
0.50
0.50
scaling ratio
1.0
1.01
0.0
0.00
0.59
scaling ratio
1.0
1
1.02
0.0
0.98
6.50
scaling ratlo
100
1.01
1.52
0 1.1C
diabetes
mimic2
compas
fico
1.0
1.0
grecision
0.6-
precision
0.6
grectsion
precition
0.8
0.4
04
0.0
0.0
0.00
0.50
1.0
3.03
1.02
0.92
0.90
1.0
1.00
3.02
scaling ratio
ratio
Scaling racio
icaling ratio
Ir
ebm
sphere
hessian
ours
Figure 7: Tradeoff between precision and the size of the approximated Rashomon set given different 0. Note that the scaling factor to the power m is volume, which is proportional to recall. Our method dominates baselines on all datasets. As the scaling ratio increases, the precision starts decrease.
small ratio to cover the true Rashomon set. For example, ~55 is needed on FICO cover the true Rashomon set. A small scaling ratio means our ellipsoid captures almost the whole Rashomon set, so that scaling by a small amount can cover the whole set. (Note that practice we would not want to do this, it would be better to use our optimized approximation to the Rashomon set to avoid false positives inside the approximation.)
GAMs with different Support Sets: In these experiments, we keep 90%, 80%, and 70% of bins trained with >, 0.0005, X2 0.001. For the baseline method, C is set to 3000, and the learning rate and the number of iterations are set to 0.001 and 2500, respectively. Since (K2)) could be very large, we first sample 10,000 different merging strategies and compare at most 100 Rp Table shows more detailed results. Merging 30% of bins for the MIMIC-II II dataset leads to empty Rpg for 10,000 merging strategies and merging bins for the Diabetes dataset also leads to empty Rps-
D brradenn-based optimization for log determinant
we find the maximum volume inscribed ellipsoid by optimizing Eq 6 Eq 6 is not guaranteed to be convex, but using the log determinant, construct
As we discussed in Section 3.1
15

---

Hessian
Dataset
x
Ours
LR+sampling
EBM+sanpling
Sphere
(our initialization)
5e-4
17.11
10.25
458.61
instant
instant
Diabetes
le-3
13.27
7.01
288.87
instant
instant
3.5e-3
12.17
6.31
229.76
instant
instant
5e-4
437.87
786.85
3142.36
instant
instant
MIMIC-II
Te-3
390.49
576.31
1862.03
instant
instant
1.5e-3
383.91
572.8
1556.97
instant
instant
5e-4
94.92
22.65
246.14
instant
instant
COMPAS
le-3
84.61
18.92
155.77
instant
instant
1.5e-3
81.05
18.38
170.49
instant
instant
5e-4
131.71
89.61
1169.73
instant
instant
FICO
le-3
119.08
73.02
975.37
instant
instant
I.5e-3
117.14
68.28
901.0
instant
instant
Table running time seconds of our method compared to baselines using logistic regression and explainable boosting machine. Baselines "hessian" and "sphere" do not require the optimization step,
so they finish instantaneously.
scaling ratio (normalized)
Diabetes
MIMIC-II
COMPAS
FICO
1.0052*
1,081
1.145
1.096
1.014
LOII
1.079
1.111
1.087
1.020
1.05C*
1.108
1.076
1.089
1.041
1.1L*
1.119
1.078
1.095
1.048
Table 3: Scaling ratio needed to ensure points sampled from the surface of the ellipsoid are outside
the true Rashomon (Ag 0.001, A2 0.001).
convex objective function. In this appendix, we explore how performance changes we optimize convex objective.
Let us first define the new objective function:
min Qw.
log(det(Q-+)) + C. [mw. 0, 0)]-
(14)
2m
Similar to Eq 66. the first term is used to maximize the volume of the ellipsoid. The volume of an ellipsoid is proportional to det(Q- 3). We normalize it by i.c. (det(Q-a)) ((((Q-=))+ = Maximizing this term is equivalent minimizing (+))
Q is positive definite since it is the quadratic form for ellipsoid. Then Q! is also positive definite and det(Q-1) (det(Q))-1 -1 We also know that det( (Q²) (det(Q)) Therefore,
(det(Q-+)) - (det (22QQ²2
(det (Q)) 2mn
(det(Q-1))) -
We can take the log on the hight-hand side term, and the objective is to minimize is thus 2m log(det(Q-1)). It is well known that the log determinant is concave. After multiplying by 2m this term is convex.
The second term penalizes the points sampled from the ellipsoid they are outside R(0). L(w) is convex w.r.t Then L(w) is also convex. Finding the maximum between convex function and a constant convex and the expectation is known to be convex, so the second term is also convex. Therefore, the objective function convex and then gradient descent to find the minimizer.
The results obtained by minimizing Eq (4) are shown Figure They are almost the same as the shows similar results compared to Table 4 These
results shown in Figures T and And Table 5
16

---

Dataset
K
precision ratio
volume ratio
time Method 2 (sec)
time Method 1 (sec)
COMPAS
15
1.01 + 0.03
1.05 = 0.07
4.84e-4 + 6.25e-5
324,85 = 8.70
COMPAS
14
1.01 E 0.05
1.12 0.12
4.76e-4 + 9.33e-5
296.31 41.46
COMPAS
13
0.99
0.01
1.20 = 0.12
4.77e-4 + 1.13e-4
272.64 6.64
FICO
28
1.01
+ 0.04
0.90 + 0.06
9.99e-4 4.54e-
377.10 + 14.95
FICO
26
0.99
+ 0.06
+ 1.10+0.13
+ 1.12e-4
374.28 10.24
FICO
24
0.99
0.02
1.29 - 0.13
5,55e-4
+ 4.04e-5
297.15 = 8.26
MIMIC-I1
35
1.00 + 0.20
0.93 = 0.10
5.89e-4 119e-411. 1.72e-4
1127.10 + 141.42
MIMIC-II
32
0.99 - 0.01
1.10 = 0.06
5.29e-4 2.76e-5
1067.50 112.9
Table 4: Precision, volume, and time comparison between blocking method (Method 2) and optimiza-
tion (Method 1). This table shows that Method 2 is faster for the same task while performing equally well.
results indicate that optimizing the convex function doesn't bring us better results and during the experiments, we find hyperparameter tuning is even harder. Therefore, we use the previous results (optimizing the determinant, not the log determinant) in the remaining experiments.
Dataset
K
precision ratio
volume ratio
time method 2 (sec)
time method (sec)
COMPAS
15
1.00 + 0.02
1.09 0.06
4.33e-4 + 7.14e-5
158.12 12 + 4.16
COMPAS
14
1.00 + 0.03
1.12 + 0.08
3.23e-4 + 6.79e-5
136.90 + 3.73
COMPAS
13
0.99 1 0.01
1.20 + 0.12
5.33e-4 + 1.28e-4
324.19 1 8.71
FICO
28
1.00 = 0.01
0.95 = 0.04
5.04e = 8.99e-4
467.49 + 5.02
FICO
26
1.00 + 0.04
1.09 + 0.08
3.44e-4 + 1.04e-4
465.45 + 5.05
FICO
24
0.99 + 0.01
1.25 + 0.10
3.02e + 2.67e-5
449.33 15.93
MIMIC-II
35
0.99 + 0.01
0.98 - 0.06
6.42e-4 = 1.91e-4
1683.11 + 17.64
MIMIC-II
32
0.99 1 0.02
1.05 0.03
5.77e 4.2le-5
1666.22 =6.57
Table 5: Precision, volume, and time comparison between blocking method (Method 2) and optimiza-
tion (Method 1) trained by optimizing Eq (14).
The results are almost the same those shown in
Table
E Variable importance range
We first show more results on variable importance range and then compare the time taken compute VI for MIP-based formulation and LP formulation.
We show the shape functions of "Glucose" when the lowest and highest variable importance achieved in Figure 3 When the importance of "Glucose" is minimized or maximized, one might be interested in how the shape function changes for other features. We show such variations in Figure g Most features keep the trend as We with some variation in magnitude.
Figure TO shows the variable importance range on the MIMIC-II FICO, and COMPAS datasets. For the MIMIC-II dataset (left subfigure), features "PFratio", "GCS". and "Age" have relatively higher VI_, which means these features are, in general, important for GAMs in the Rashomon set. For the FICO dataset (mid subfigure), features "ExternalRiskkstiin qexcl7days" have higher VT_ either fixing or not fixing other coefficients, indicating that these two features are important. Feature aprior_coumte in the COMPAS dataset has slightly higher VI_ than feature "age."
and "MSinneMostReeenttnn
As we mentioned in Scction 4.1 the lower bound of variable importance V V__ obtained by solving a linear programming problem with a quadratic constraint, while to get the upper bound of variable importance we need to solve mixed integer programming problem. We use Python CVXPY package [35] 36) for solving LP problems and Cplex 12.8 for MIP problems. Note that as long as M > 2 w. the MIP formulation is valid. We set M to 200. This is usually large enough to bound the absolute value of each coefficient given that we have the &a penalty. Since sol a MIP problem usually time-consuming we also propose another way to get VT4. Note that it's the absolute value in the objective that restricts us use LP solver. Therefore, enumerate all positive-negative
17

---

A,
0.0005
disbetes
mimic2
compas
fico
10
2.0
brecsoo
preciiion
Drecsion
000
0.0
0.003
0.03
0.355
0.01
0.33
0.000
0.01
0.55
0.066
0.31
(9 C'Ve"
VC"
18
10 6'V
Ag
0.001
diabetes
mimic2
compas
0.3
0.8
0.2
0.01
0005
0.33
0.5
0.33
0.05
0.05
10- C""C"
HC
(0- C'
(9- E'
A₂ 0.0015
diabetes
mimic2
compas
fico
grection
proctsion
0.055
0.53
0.5
03
(9-c')0
(0-C
y'v
(6-
ebm
sphere
hessian
ours
(a) Precision of the approximated Rashomon set as a function of 0.
0
1.005CC*
diabetes
mimic2
compas
fico
3.0
3.5
2.0
orecsion
0.8
08-
0.58
2.09
scaling
1.0
1.00
0.0
0.90
Scaling ratlo
1.0
1.33
1.02
0.93
0.99
Scaling
1.0
101
1.02
scaling ratio
0
1.01."
diabetes
mimic2
compas
1.0
3.0
prccsion
5.0
0.0
1.01
3,02
0.60
009
0.00
Sceling
1.0
1.03
1.02
scaling ratlo
0 1.05.
diabetes
mimic2
compas
fice
1.0
3.8
0.
precsion
0.6
precsion
02
1.0
1.3
2,02
6.6
0.60
14
140
5c5 (ng ratio
scaling rallo
scaling ratio
scaling (2030
1.1.C*
diabetes
mimic2
compas
2.0
Srocsin
0.8
0.0
0.50
3.50
3.0
1.03
1.02
0.94
0.90
1.5
1.01
1
0.51
0.99
1.0
1.01
1.52
scaling
scaling rabio
scaling ratio
ebm
sphere
hessian
ours
(b) Tradeoff between precision and the size of the approximated Rashomon given different 0. The scaling factor to the power m volume, which is proportional to recall.
Figure& 8: Precision and volume of the approximated Rashomon set obtained by optimizing Eq 14) They are almost the same as the results shown in Figure and
18

---

shape function (VI-)
shape function (Vi-)
shape (VI-)
shape function (VI,)
da. ttst day, fret
dt. stcl
win, Jrst Bel
de tustii
togit
groett begy
predicted togit
predicted tegit
Mo. [94]
(ntt 14
-03
G Guuose
Giucose
Pregnancies
Pregnanties
shape function (VI_)
shape function [VI.)
shape function (V)_.
function (VI.)
Mu (not (o)
presicted ogt
preaicted age
0.1
-1.5
0.2
-2.0
-0.5
-25
-25
-1.
-1.0
302
125
130
125
20
100
B 1eodPressur
Sosuerassure
ScinThiccmsss
SkinThlkkeess
function (VI-)
shape function (VI,)
shape function
shape function (Vi..
May fnst
ein, (nst
predi loght
credd toght
poddetee legi
403
80c
201
45o
800
600
Insulin
Insulin
BNI
BHI
shape function (VI_)
shape function (V)..
shape function (W_.)
shape function (VI.)
3.0
peaderee loog
predicted logit
predictdd
(1ct
S111 Ui. 2not Bx
-18
-10
akp (nat
(X0C
60
20
35
60
10
1.
30
53
75
30
53
60
10
Diosetesled gresFunction
Figure 9: Shape functions of other features in the Diabetes dataset when the importance of "Glucose" minimized maximized.
mimic2
fico
compas
NetastaticCancer
silrubin
Netracctnnioeevnn@uuren
prios, couns.counc.
Renal
Percentrstal
Pfreto
Kumfatbdfectoryf Tradess
AverageNi-Fil
Age
Extema/RisEsttmatt
600 ols
range
variable importance range
variable importance range
fu other coets
Sa other conts
Figure 10: The variable importance range of the MIMIC-II, FICO, and COMPAS datasets, Each line connects VI_ and V44. (A, 0.001, X2 = 0.001, 0 = 1.01C*)
19

---

Algorithm 2 G- --rmm Q w.)
Input: proportion of samples each bin of feature j πy E RB,, parameters of the ellipsoid Q Rmxx ,W E IR
Output: point inside the ellipsoid w E R
1: Obj -00
2: for ( // try all positive-negative combinations for wjik, k (0....B
// solve the LP problem
3: Obj.w w & maxw, (TT @ I)" w, such that (w -w:) w.)" we) S
4: if Obj O0b":
5:
w" t w. Obj t Obj
6: return wt
predicted logit
10-
predicted llgit
predicted logit
predicted logit
-30
-1.
Pregnancies
BloodPressune
SkanTTiickesss
gredicted logit:
predicted logit
logit
predicted logit
2.5-
-25-
-10
-15
6
200
400
620
302
20
40
60
70
0.
05
2.
15
30
25
50
90
70
60
Insulin
BMI
RlabetesFadigreesunnti
Age
Figure 11: 100 different shape functions of the Diabetes dataset. The shape function at we is colored in gray. (A, 0.001, A2 0.001, 0 1.01C*)
combinations of wj {0, B; 11, and then solve LP problem with the sign constraint enforced (see Algorithm 2) Table compares the running time by solving both MIP and LP problems for V4+. The runtime required to solve LP problems with the sign constraint enforced is usually less than that required to solve MIP problem.
1}.
Fix other coefs
Not fix other coef's
Dataset
MIP
L.P
MIP
L.P
Diabetes
18,876 + 10.025
4.780-14.996
35.597 + 17.413
1.961 + 2.164
MIMIC-II
11.450
6.259
0.4974 0.166
14.977 1 12.921
0.265 1 0.139
COMPAS
23.193 + 0.176
4.138 = 0.033
46.782 + 5.709
3,043 = 0.017
FICO
11.749 + 7.841
0.816 + 0.766
30.133 + 7.11
0,425 + 0.457
Table Time comparison in seconds between solving MIP and I.P problem with the sign constraint enforced for VT+
F
Shape functions of GAMs in the Rashomon set
Next, show diverse set of coefficient vectors sampled from the approximated Rashomon set. Figure [1 and Figure 12 depict 100 coefficient vectors (in red) sampled from the ellipsoid and w. (in gray). the center the optimized ellipsoid of the Diabetes and MIMIC-II dataset, respectively. Various different red lines in each subfigure indicate that the approximated Rashomon set contains many different coefficient vectors. And these models are actually within the true Rashomon sets. Many of them may generally follow similar patterns as we can see from the trend of these shape functions, while some of them may have some variations (see Figure 5). In summary, using the approximated Rashomon set, easily get diverse set of shape functions for each feature.
20

---

predicted logic
1.5
predicted logit
1.0
losit.
predicted logit
1.0
0.5
800
0.2
-3.5
-1.0
-3.75
2:30
-03
23
10
350
10:15
23
4.
162
156
250
50
380
Age
GCS
SBP
HR
1.0
0.75
0.0
predicted logit
predicted logit
0.25
predicted
02
predicted
002
0.0
-6.30
-0.5
500
1000
1500
2000 2200
Temperature
Pfratio
Renal
Urea
1.0
presicted logie
logit.
predicted logit
1.80
logit:
2.5-
1.5
1.0-
1.5
1.0
05-
0.25
0.
0.0-
-05
200
430
660
903
10
85
30
60
-0.50
120
320
340
100
-1.0
20
8
WBC
CO2
Na
Bilirubin
predicted logi:.
prediitted logit:
1.0
0.0
52
04
0.6
5.0
1.5
0.0
02
0.1
0.6
50
1.5
MelastaticCancee
Lymphora
Figure 12: 100 different shape functions of the MIMIC-II dataset. The shape function at is colored
gray. (A, 0.0002, A2 0.001, 0 1.01
G User preferred shape functions
In real applications, users might have preference for shape functions that are consistent with their do- main knowledge. Our approximated Rashomon set makes it easy to incorporate such user preferences by finding a model within set that satisfies some constraints. We show two case studies here.
Figure 11 shows that jump occurs when blood pressure is around 60. One possible user request might be making this blood pressure shape function monotonically decreasing. By solving the quadratic programming problem with linear constraints, we can get the shape functions colored in yellow in Figure 13al The updated shape function for "BloodPressure" is monotonically decreasing. And the shape functions for other features are also updated. Almost all of them follow the trend in w. (in gray), the center of the optimized ellipsoid, with small changes in magnitude. Note that this optimization problem is solved in 0.04 seconds.
Diabetes:
One might also be interested in making the shape function of "BMI" monotonically increasing. By solving the optimization problem, we can get the shape functions shown in Figure T3b) The updated shape function for "BMI" is monotonically increasing (in yellow). The sudden decrease that occurs "BMI"=40 is removed by connecting the left step. Similar to Figure T3a the updated shape functions of other features follow the trend in w. (in gray), the center of the optimized ellipsoid, with small changes in magnitude. And this optimization problem is solved in 0.0004 seconds.
Sometimes a monotonic constraint is not what users want; they might have more specific preferences on the shape functions. Here we show some examples using hypothetical shape functions. Figure [4 extends the visualizations in Figure 4 It shows shape functions after imposing two different requests on "BloodPressure."
Figure 15 shows shape functions after imposing two different requests on "BMI." In Figure 15a the requested shape function of "BMI" (colored in red in the top-left subfigure) shifts below the original shape function but maintains monotonically increasing pattern. While the closest shape function
21

---

adding constraint
predictdd logit
-2.
preocceo
02
before editing
122
BloodFressure
SkinThickness
predictsa logr.
predictos
predicted
2 1
-10
Insulin
8M4
DistetsSPodggnooo Funcccnn
Age
(a) Shape functions of the Diabetes dataset with the monotonic constraint the "BloodPressuree" (in yellow). The optimization problem is solved in 0.04 seconds.
00
prod logt
85
lagil
predicted logi-
predicied fool:
predic
25
15
200
65
66
85
350
22
62
Pregnancios
Glacose
[RoodPhessure
Str1 ness
logit
predictos logi:
after adding monotonic
peedcooa togs
togs:
constraint
0:1
before editing
2.00
800
(b) Shape functions of the Diabetes dataset with the monotonic constraint on the "BMI" (in yellow). The optimization problem is solved in 0.0004 seconds.
Figure 13: Shape functions with monotonic constraints.
in the approximated Rashomon set is also below the original shape function, it is not necessarily monotonically increasing. Instead, it is more likely to follow the trend observed in the original shape function, as we aim to minimize the Euclidean distance between the requested shape function and the shape function in Å. Shape functions of other features change only slightly magnitude. Figure 15b shows a different case. The requested shape function of "BMI" forces certain steps to have the same coefficient. However, after solving the QP problem, the updated shape function shown in green in the top-middle subfigure a combination of theshape function before editing and the requested shape function.
MIMIC-II: In Figure 12 we can see jumps in several shape functions. For example, "PFratio" has sudden jump around 330, and "Bilirubin" has jump close to O. PFratio is measurement of lung function; measures how well patients convert oxygen from the air into their blood. And bilirubin is a red-orange compound that breaks down heme. The bilirubin level reflects the balance between production and excretion. The elevated levels may indicate certain discases. Missing values commonly exist in the real dataset and imputation is widely used. (28) shows these jumps are caused by mean imputation, and have no physical meaning. We can impose monotonic constraints on these two features simultaneously. We want the shape function of "PFratio" to be monotonically decreasing, while the shape function of "Bilirubin" is monotonically increasing. Figure 16 shows the shape functions after optimization. The shape functions of "PFratio" and "Bilirubin" satisfy the request as shown in the inset plots, and the shape functions of other features are only slightly changed.
Suppose user prefers to remove the jump in the shape function of "PFratio." One way is to remove the jump while keeping the last step unchanged (top-left subfigure Figure[T7p). Fortunately, by solving Problem (12). find that the specified shape function within the Rashomon set (shown by
22

---

0.1
user requires
predictnd logit
after editing
0.4
predicted
0.5
gredicted
02
0.0
-0.5
before editing
-03
before editing
-0.1
20
40
60
100
120
40
60
80
100
120
00
2.5
5.0
7.5
10.0
12.5
15.0
17.5
BloodPressure
BloodPressure
Pregnancies
logit
logit
-05
predicted logit
-1.0
-15
-2.
SkinThicknees
predicted logit
predicted logit
prodicted logit
-1.
10
20
10
50
00
70
0.0
0.5
1.0
1.5
2.0
25
20
30
40
50
60
BM
DiabeetesPedigreefuncccc
Age
(a) Case
before editing
before editing
predicted logitt
predicted logit
predicted logit
0.7
requires
after editing
60
50
100
150 175
BloodPressure
BloodProssure
Pregnancies
predicted logit
-0.5
predicted logit
presicted
-1.5
150 :35
200
80
50
100
200
400 Insulin
600
Glucose
SkinThickness
1.0
predicted logit
predicted logit
0.5
predicted
DiabottePodigreeFFnnttoon
0.5
2.0
2.3
2.0
255
30
Age
(b) Case
Figure 14: Shape functions the Diabetes dataset after hypothetical shape function on "Blood- Pressure" is requested. The red curve in the top-left subfigure the required shape function. The shape function colored in green in the top-middle subfigure is the closest shape function within Å..
23

---

before editing
before editing
predicted logit
predicted logit
predicted lagit
04
requires
after editing
-0.1
20
30
43
52
70
23
30
40
50
60
70
0.
5,0
7.5
10.0
12.5
15.0
175
BMI
BM1
Prognancies
predicted logit
predicted logit
1.0
predicted lagit
0.5
0.5
-55-
0.0-
-0.5
25
50
75
100
125
150
175
200
60
60
300
120
20
40
50
50
100
Glucose
BloodPressure
SkinThickness
predicted lagit
predicted logit
predicted logit
-2
-1.0
250
400
600
800
0.0
05
100
15
23
20
30
60
Insullin
Diabetes Pedigrrestcit
Age
(a) Case
before editing
before editing
0.4
predicted logit
predicted logit
predicted logit
0.3
0.2
user requires
af tee editing
10
30
40
60
10
20
40
50
70
0.
2.5
5.0
75
10.0
125
15.0
175
BMI
BMI
Pregnancies
1.5
predicted logit
1.0
predicted logit
-0.5
0.5
-1.
-1.
0.0
-05
20
40
60
80
100
120
20
40
100
BlodPPeessure
SkinThickness
1.0
predicted logit
0.5
1-
0.0
-4-
0.5
-5
-1.
DisbetesPedigreeeunnttoon
Age
(b) Case
Figure 15: Shape functions on the Diabetes dataset after hypothetical shape function on "BMI" is requested. The red curve top-left subfigure is the required shape function. The shape function colored in green the top-middle subfigure is the closest shape function within f
24

---

077
1.5
06
predicted logit
predicted logit
03
1.0
02
0.4
predicted
0.5
0.3
0.0
-37
0.1
02
-04
0.0
Age
0.8
1.0
after adding monotonic
0.6
predicted logit
0.4
predicted logit
0.5
constraint
0.3
0.0
0.
62
01
before editing
0.0
0.0
2500
HR
Temperature
Pfratio
0.
1.50
predicted looit
0.6
predicted logit
gredicted logit:
1.25
0.4
1.00
0.75
0.2
0.50
-21
0.25
-32
0.00
0.5
3.0
Renal
2.0
25
30
35
4.0
-0.25
200
400
WBC
600
180
0.6
0.8
predicted logit
0.4-
3.5
after adding monotonic constraint
predicted
2.0
0.4
0.5
0.2
0.1
before editing
CO2
50
100
110
730
130
Na
140
190
1/6
170
Bilirubin
1.0
predicted logit
0.5
predicted logit
0.8
0.6
0.4
6.44
0.0
0.0
MetastaticCancer
lymphoma
Figure 16: Shape functions of the MIMIC-II dataset with the monotonic constraints the "PFratio" and "Bilirubin" (in yellow).
25

---

3
w sampled from R
Dataset
9 (constant & ()
accuracy
auc
accuracy
auc
1.005
0.683 + 0.010
0.744 0.004
1.01
0.683 = 0.010
0.742 + 0.005
COMPAS
0.696
0.748
1.05
0.668 + 0.017
0.724 + 0.012
1.1
0.649 + 0.026
0.704 + 0.021
1.005
0.717 0.003
0.791
0.001
1.01
0.716110.00 0.004
0.790 + 0.002
FICO
0.720
0.792
1.05
0.708 + 0.008
0.780
0.006
1.1
0.700 0.010
0.770 - 0.009
1.005
0.761 + 0.004
0.818 0.002
1.01
+ ::000:00000 0.005
0,818 + 0.003
Diabetes
0.760
0.819
1.05
0.758 + 0.011
0.816 0.006
1.1
0.755 + 0.014
0.814 + 0.009
1.005
0.886 + 0.001
0.803 0.002
1.01
0.886 + 0.001
0.802 + 0.002
MIMIC-II
0.886
0.803
1.05
0.885 1 0.002
0.794 0.005
1.1
0,884 = 0.003
0.784 = 0.009
Table 7: Test accuracy and AUC comparison between w and w sampled from the approximated Rashomon set with respect to different Os.
the green curve in the top-mid subfigure in Figure 17a). Another user might not like this idea and prefers to remove the jump by connecting to the left step (Figure However, this specified shape function is not within the Rashomon set, and we find the closest solution in green, which still has a small jump at 330 (see the inset plot). Different user-specified shape functions lead to different solutions. The Rashomon set serve small but computationally efficient space in which users can find a model that is closest to their needs.
H Test performance
We now show the performance of models sampled from our approximated Rashomon set. We compare the accuracy and AUC between w" and ws sampled from it on the four datasets with different values of and results are shown in Table 7 We sample 1000 from R and show the average and one standard deviation. The larger value of 0 leads to larger Rashomon set, which means we allow models with higher loss. Therefore, as increases, both accuracy and AUC decrease. But when the constant is slightly larger than 1, such as 1.005 and 1.01, the test performance of models sampled from ft usually covers the performance achieved by w" in one standard deviation. This means general our approximated Rashomon set can return diverse set of models without compromising the performance.
26

---

user requires
editing
prodeemet leght
1.0
1.3
predeced lnggt
00
-1
before editing
before editing
-15
1248 Pfootio
2000
2083
500
2008
1200
2000
2500
Age
04
133
predectd logk
04
0.3
20
350
236
200
256
306
1.6
Temperature
loge
2460 cled loget
1.0
0.-
20-
006
02
-0.2
350
406
350
800
000-
150
250
300
133
120
330
180
pesteceedoop
predected
00
02
06
15
Mer CaberrCanee
(a) Case
user requires
ofter editing
0z
gred ttec Ing
paotic
before editing
before editing
553
1890
2503
2006
PProto
Age
055
0.0
predected
prodicted
0.
30
50P
HR
Temperature
Benal
1.50
0.-
2.5
loggs
1.00
predictdd loght
180
356
330
683
Usia
C02
14
1.0
predictd
1
1.4
Skimibin
Mettestaticaance
kamphoma
(b) Case
Figure 17: Shape functions the MIMIC-II dataset after hypothetical shape function "PFratio" requested. The red curve in the lep- subfigure is the requested shape function. The shape function colored in green in the top-middle subfigure is the closest shape function within f.
27


Reviews:
[{"review_id": "Kdl0bwydAN", "reviewer": "Reviewer_r4Ee", "summary": "This paper presents a novel algorithm (and an additional variant) to learn the Rashomon set of sparse generalized additive models (GAMs) approximated by an ellipsoid in the parameter space. The paper also presents 4 interesting use cases of the proposed algorithm. It includes a rich experiments section. In the experiment section, the authors compared the proposed algorithm with a few others and showed that the proposed algorithm has very favorable results in terms of both precision and volume. The experiment section also includes the demonstrations of the proposed algorithm in several use cases.", "questions": "1. line 89, p. 3: You mentioned \"$\\omega$ also includes $\\omega_0$\". I assume this is just to say that the notation $\\omega$ will include the bias term $\\omega_0$. However, we usually don't apply $L_2$-regularization to the bias term, so regularization weight $\\pi_0$ for the bias term still remains zero, right?", "limitations": "n/a", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Interesting problem, rigorous formulation, novel algorithm, clear presentation, and favorable experiment results.", "weaknesses": "The problem this papers studies is a niche in the general ML problem space.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review! We really appreciate it! See below for our response to your questions.\n\n> Q: The problem this paper studies is a niche in the general ML problem space.\n\nIn fact, GAMs are a widely used form of interpretable predictive model. The vast majority of healthcare scoring systems (e.g. the APGAR score, and CHADS2 scores - the most widely used models in medicine) are GAMS. In criminal justice, GAMs are used extensively for parole and bail predictive modeling. So the problem has very general applications and makes it easier to construct models in a huge number of domains.\n\nIn addition, our paper is closely connected to multiple important topics in machine learning, including interpretability, uncertainty quantification, and human-model interactions. And we believe our work will have a broad impact in the domain.  \n\n\n> Q: line 89, p. 3: You mentioned “$\\omega$ also includes $\\omega_0$”. I assume this is just to say that notation w will include the bias term $\\omega_0$. However, we usually don’t apply $L_2$-regularization to the bias term, so regularization weight $\\pi_0$ for the bias term still remains zero, right?\n\nRight. We include $\\omega_0$ in $\\omega$ just for notational simplicity. In our implementation, we don’t include $\\omega_0$ when calculating the $\\ell_2$ penalty. Thank you for pointing this out. We will make this point clearer.\n\n\n"}, {"review_id": "lmFwMvOBNi", "reviewer": "Reviewer_8bE1", "summary": "this paper designs a new editable gam model. They leverage the ellipsoid to approximate the optimal solution set. The model allows to edit the model according to different use-cases and requirements.\nThe edited parameter could be solved by constrained quadratic problem. They experiment the model and edited use-cases at four datasets.", "questions": "1, how many bins do you used in EBM? \n2, how to use the proposed method when the final model is the bagging of multiple GAMs？ ", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.providing a novel perspective of constructing editable GAM\n2.proving the upper and lower bound of the variances' importance\n3. showing the specific result of edited use-case of multi-datasets", "weaknesses": "1. the procedure of approximating the roshomon set has high complexity. The experiments don't evaluate the overhead   of this part.\n2. though the quadratic problem is easy to solve by programming, the constraint of ellipsoid and high dimensions of features may make the procedure to be high complexity. It needs to evalute the overhead in the experiment. Besides, as Q has large dimension, it is expensive to store it. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review! We really appreciate it! See below for our response to your questions.\n\n> Q: the procedure of approximating the Rashomon set has high complexity. The experiments don't evaluate the overhead of this part.\n\nWe do show the running time of approximating the Rashomon set in Table 2 of Appendix C. In most cases, our proposed method has a run time slightly longer than logistic regression with bootstrapping but shorter than the EBM baseline. Baselines “hessian” and “sphere” do not require the optimization step, so they finish instantaneously. For this table, we ran the gradient descent on a CPU, whereas if we had used GPU, it would be at least 10x faster.\n\n> Q: though the quadratic problem is easy to solve by programming, the constraint of ellipsoid and high dimensions of features may make the procedure to be high complexity. It needs to evaluate the overhead in the experiment. Besides, as Q has a large dimension, it is expensive to store it.\n\nIn terms of the linear programming and quadratic programming for the applications mentioned in our paper, Table 6 in Appendix E shows the time needed to calculate the upper bound of the variable importance. The time needed to obtain the lower bound of the variable importance is even faster than getting the upper bound, usually within 0.002 seconds. \n\nTo find a model with monotonic constraints (see Figure 4a), the optimization problem is solved in 0.04 seconds. Another example can be found in Figure 13 in Appendix G. The time used to get the model in Figure 4b is 0.0024 seconds and in Figure 4c is 0.0022 seconds. Thank you for pointing this out. We will include timing information in the caption. \n\nOur goal is to find the Rashomon set of sparse GAMs. Sparse GAMs are easy to interpret and have performance on par with more complicated models. With the sparsity penalty, Q usually doesn’t have extremely large dimensions. Even if we have (for instance) 1000 bins, Q is a 1000x1000 matrix, which is only 8MB if we store the matrix elements as double float numbers. It is much more memory efficient than many other machine learning models such as deep neural networks.\n\n> Q: how many bins do you use in EBM? \n\nFor a fair comparison, we use the same number of bins in EBM as the reference model, which is obtained by running FastSparse[1], a sparse GAM. To achieve this, we turn continuous features into categorical, where the categories are determined by which bin in the reference model the feature value belongs. We then fit an EBM on these categorical features by setting feature_types to “nomials”. \nYou might want to know how many bins are used in the reference model. In our experiments, the number of bins ranges from 8 to 61. \n\n[1] J. Liu, et al., Fast sparse classification for generalized linear and additive models. Proceedings of machine learning research,151:9304,345 2022.\n\n> Q: how to use the proposed method when the final model is the bagging of multiple GAMs？\n\nOur method can find the Rashomon set for a single GAM. If the base model is the bagging of multiple GAMs, one could simply aggregate them into a single GAM model (since the average of GAMs is still a GAM) and our algorithm can find the Rashomon set of this aggregated model. \n\n\n\n"}, {"review_id": "UUVBmRxect", "reviewer": "Reviewer_SBSW", "summary": "In contrast to the conventional machine learning paradigm, which typically yields a single model, Rashomon sets comprise a collection of near-optimal models. These sets allow users to select a model based on their specific preferences. This paper focuses on constructing Rashomon sets for sparse generalised additive models (GAMs). Furthermore, the paper demonstrates that by constructing Rashomon sets, users can estimate variable importance, enforce monotonic constraints, and interact with the shape functions of GAMs.", "questions": "I have difficulty understanding the distinction between fixed support sets and different support sets in Section 3.2. Could you please provide further clarification or explanation on this matter?", "limitations": "I appreciate the solidity of the paper, although I must admit that I personally find it challenging to comprehend its content. I want to emphasize that this difficulty arises from my own limitations as a reviewer, rather than any shortcomings of the paper itself. As I am not familiar with the surrounding literature in this field, I can’t provide a  judgment.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper addresses a novel problem of constructing Rashomon sets for sparse generalized additive models (GAMs). It highlights the scarcity of prior work in this specific area, emphasizing the novelty and significance of the research.", "weaknesses": "While the paper appears to be well-executed, I must admit that I lack familiarity with the concepts of Rashomon sets and generalized additive models. As a result, I am unable to provide an in-depth review or offer a confident assessment of the paper's claims. This limitation is due to my limited knowledge in these specific areas, despite my previous research experience in explainable AI.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review! We really appreciate it! See below for our response to your questions.\n\n> Q: I have difficulty understanding the distinction between fixed support sets and different support sets in Section 3.2. Could you please provide further clarification or explanation on this matter?\n\nWe are sorry that the paper is not understandable to you. We truly hope everyone can understand the paper even without much background knowledge. Would you mind pointing out which parts are unclear or difficult to understand? We will make these points clearer.\n\nWe can view continuous features as collections of bins. For example, blood pressure (bp) is a continuous feature used for prediction. It ranges from 0 to 125. A sparse GAM with fixed support set looks like this \n\n$g(bp) = b_0 + b_1 \\times \\mathbb{1}[ 0<=bp<60] + b_2 \\times \\mathbb{1}[60<=bp<65] + b_3\\times \\mathbb{1}[65<=bp<100 ] + b_4 \\times \\mathbb{1}[100<bp].$\n\nBinary features $0 \\leq bp < 60$, $60 \\leq bp<65$, $65 \\leq bp<100$, $100<bp$ consist of a fixed support set. And we can find the Rashomon set based on this fixed support set by optimization in Section 3.1. \n\nGiven a fixed support set, we can get different support sets by merging bins. For example, \n{$0 \\leq bp<65$, $65 \\leq bp<100$, $100<bp$} is a set. {$0 \\leq bp<60$, $60 \\leq bp<100$, $100<bp$} is another support set. {$0 \\leq bp<100$, $100<bp$} is also a support set. While we can use optimization in Section 3.1 to construct a Rashomon set for each of these three subsets, it seems not very efficient. Therefore, we propose the method in Section 3.2 to quickly find the Rashomon set for each of these support sets by using information from the approximated Rashomon set of the fixed support set. \n"}]

Author rebuttal:

Meta-review: This paper received only three reviews, despite my efforts to recruit more reviewers. This meta-review is based on my reading as well.

The reviewers unanimously accepted the paper. I agree with this decision. This well-written and clear paper addresses an important and burgeoning topic in machine learning: exploring the "Rashomon set" of good models. The authors focus on GAMs and provide a compelling algorithm for approximating the Rashomon set using ellipsoids. The algorithm is clear and intuitive.

This paper will make a very nice addition to NeurIPS and can impact how GAMs are used in practice.

In addition to comments noted in the reviews, my minor suggestions for the authors are (1) improve the clarity of their figures, both in terms of presentation (e.g., font size) and explanation, and (2) a disclaimer that there are several recent calls for refraining from using the COMPAS dataset, though I recognize its use as a benchmark.
Decision: Accept (poster)