# Obtaining Explainable Classification Models
using Distributionally Robust Optimization

Sanjeeb Dash, Soumyadip Ghosh, João Gonçalves, Mark S. Squillante
IBM Research, Yorktown Heights, NY 10598

###### Abstract

Model explainability is crucial for human users to be able to interpret how a
proposed classifier assigns labels to data based on its feature values. We
study generalized linear models constructed using sets of feature value rules,
which can capture nonlinear dependencies and interactions. An inherent trade-
off exists between rule set sparsity and its prediction accuracy. It is
computationally expensive to find the right choice of sparsity – e.g., via
cross-validation – with existing methods. We propose a new formulation to
learn an ensemble of rule sets that simultaneously addresses these competing
factors. Good generalization is ensured while keeping computational costs low
by utilizing distributionally robust optimization. The formulation utilizes
column generation to efficiently search the space of rule sets and constructs
a sparse ensemble of rule sets, in contrast with techniques like random
forests or boosting and their variants. We present theoretical results that
motivate and justify the use of our distributionally robust formulation.
Extensive numerical experiments establish that our method improves over
competing methods – on a large set of publicly available binary classification
problem instances – with respect to one or more of the following metrics:
generalization quality, computational cost, and explainability.

##  1 INTRODUCTION

The output of many popular machine learning (ML) methods for binary
classification such as random forests are hard to interpret and analyze. In
many applications that have societal impact, e.g., in areas such as criminal
justice and medicine, the use of ML models that are not interpretable is being
questioned (Rudin,, 2019; Rudin et al.,, 2022). In recent years, there has
been a lot of emphasis on explainability and interpretability in ML models.
This has led to methods such as LIME (Ribeiro et al.,, 2016) that explain the
predictions of a classifier by approximating it locally with an interpretable
model and also led to revisiting rule sets (Dash et al.,, 2018; Wang et al.,,
2017), rule lists (Lakkaraju et al.,, 2016), and decision trees (Hu et al.,,
2019) as interpretable classifiers. Ensembles of rules (i.e., linear
combinations of rule outputs) are considered by Friedman and Popescu, (2008)
instead of rule sets (i.e., disjunctions of rules), though there is less
emphasis on the sparsity of the ensemble. Dhurandhar et al., (2020) construct
pre-selected interpretable ensembles of models hand-picked for explainability,
and optimize the parameters of the models as well as the weights of each
ensemble member by comparing their output on a training dataset to that of a
good non-interpretable model. Lakkaraju et al., (2016) present a user-study
that compares the ability of humans to interpret rule sets and rule lists and
concludes that rule sets are more interpretable than rule lists.

In the above-mentioned papers on rules sets and decision trees, model
complexity is considered an important parameter for interpretability. For rule
sets, one measure of complexity (Dash et al.,, 2018) is the number of rules
plus the number of terms in the rules. Experimental results in (Schmid et
al.,, 2017; Muggleton et al.,, 2018) show that, for logic programs (which
include rule sets and decision trees), inspection time (the time taken by
humans to study and understand the program before applying it) is negatively
correlated with human predictive performance applying the logic program. It is
clear that higher model complexity will in general lead to more inspection
time. Thus compact rule sets are considered to be more interpretable than rule
sets with many rules, yet the latter usually have more predictive accuracy,
and so one needs to trade off model complexity and predictive accuracy.
However, more complex models can lead to the overfitting of training data and
subsequent poor generalization over unseen data assumed to be sampled from the
same underlying (unknown) distribution. Cross-validation is used by Dash et
al., (2018) to choose a good model complexity. This is computationally
expensive since each training with a given complexity bound requires solving a
combinatorial optimization problem.

In this paper, we focus on sparse ensembles of rule sets as binary classifiers
(for binary data). These generalize both rule set classifiers and ensembles of
rules (Friedman and Popescu,, 2008). We define the complexity of our ensembles
to be the sum of the complexities of the rule sets used in the ensemble (see
Section 2.1 for a definition of rule set complexity). We limit model
complexity to improve interpretability while trying to obtain more accurate
classifiers than rule sets of similar complexity levels. To avoid an expensive
cross-validation step to choose the appropriate model complexity and to obtain
good generalization, we use a distributionally robust optimization (DRO)
approach.

Recently, Ghosh et al., (2021) showed that DRO-based statistical learning
solutions generalize as well as those obtained by regularized learning
formulations that need computationally expensive fine-tuning via cross-
validation. They provided efficient algorithms to solve DRO problems over
models with continuous variables. However, their DRO formulation does not
translate directly to the cases of binary-valued data and rule set models (see
remarks at the end of Section 2.3). We extend this DRO approach in a couple of
important ways in order to apply it to binary-valued data and rule set models.

First, we consider the selection of our classifier from the class of all
_convex combinations_ of rule sets. Second, we set the complexity of such a
model to the sum of the complexities of the rule sets. Our formulation then
selects a convex combination that is minimax optimal for a DRO formulation
while satisfying a sparsity constraint. Convex combinations of binary-valued
models are an important class of models that have been considered in the
context of ensemble models (Sollich and Krogh,, 1995), random forests
(Breiman,, 2001), boosting (Freund and Schapire,, 1997) and their variants.
These approaches have been favored by practitioners for their superior
performance in out-of-training inference, and indeed they represent the state-
of-the-art for many of the datasets considered in our experiments (see Section
4). However, our approach fundamentally differs from previous work in giving
primary consideration to the sparsity of the chosen combination, while
existing methods produce dense and hard to interpret models.

Our Contributions. We construct our sparse ensemble of rules sets using two
formulations that are solved in coordination with each other. The first is a
distributionally robust combinatorial optimization problem that selects a
_dense_ collection of _sparse_ DNF classifiers that minimizes the worst loss
over a carefully chosen distributional ball. The chosen DNF models obey a
strict complexity bound. The solution approach utilizes column generation to
efficiently search through the space of all feasible DNF models. The second
smaller integer program (IP) selects a sparse _convex_ combination from the
dense set generated by the first formulation. This IP is easily solved with
standard tools, and can be used within fine-tuning procedures to select the
best complexity constraint for the chosen convex model. In Section 2, we
provide upper bounds on the generalization performance of individual rule sets
as well as their convex combinations in terms of their sparsity
characteristics. The relevant DRO literature is reviewed in Section 2.3 and
motivates the use of this approach in our algorithm, described in Section 3.
In Section 4, we provide the results of our extensive experiments on binary
classification models chosen from the public domain. We demonstrate that the
sparse convex ensembles of rule sets we construct (a) use less computational
effort than competing state-of-the-art methods; (b) while matching, and in
some cases improving, the generalization performance over both the sparse rule
sets of Dash et al., (2018) and the best state-of-the-art methods; and (c) are
of relatively low complexity leading to easy interpretability.

##  2 THEORETICAL ANALYSIS

In this section we present theoretical analyses and establish related results
that motivate and support our algorithm provided in Section 3. We start in
Section 2.1 with establishing the notation needed to describe the problem
setting and our algorithm. Section 2.2 explains why existing rule set
optimization methods require expensive computation to ensure good model
generalization. Section 2.3 provides an introduction of the recent
developments of DRO techniques to alleviate this concern. However, the
implementation of such techniques for rule set models is complicated by the
combinatorial nature of the model class. In Section 2.4, we propose the use of
convex combinations of rule sets both as a more flexible model that can
improve generalization and as a more suitable candidate for DRO training
formulations.

###  2.1 Definitions

Binary classification is the problem of finding a model $F(x)$ that correctly
classifies a data point $x$ as having label $y\in\\{0,1\\}$. We assume that
the data $x$ can be represented as an ordered set of $d$ binary valued
features, i.e., $x\in\\{0,1\\}^{d}$. Data with integral or categorical
features can be transformed into this setting by using, for instance, one-hot
encoding. Real-valued features can be encoded in this manner by binning the
range of all observed values. The true distribution $P^{t}$ that generates the
pair $(x,y)$ is seldom known, and model fitting is done based on a (finite)
training dataset ${\mathcal{S}}=\\{(x_{1},y_{1}),\ldots,(x_{N},y_{N})\\}$ of
size $N$. We refer to $P^{0}=(1/N,\ldots,1/N)$, the probability mass function
(pmf) that assigns equal weights to the points in the training dataset, as the
empirical pmf. The quality of a model $F$ is determined by the loss
$l(F(x),y)$ experienced at data $(x,y)$. In the sequel, we are interested in
minimizing expected loss with respect to pmfs $P$ over dataset
${\mathcal{S}}$: $\mathbb{E}_{P}\,l(h(x),y)\,=\,\sum_{i}l(h(x_{i}),y_{i})\cdot
P_{i}$.

Let $\mathbb{I}(\cdot)$ return 1 if its argument is true and 0 otherwise. A
_simple rule_ is a function of the form

$r(x)\;{:=}\;\mathbb{I}(x^{j}=0)\mbox{ or }r(x)\;{:=}\;\mathbb{I}(x^{j}=1)$

where $x$ is a datapoint and $j$ is an index between 1 and $d$. A _rule_ is a
conjunction $t(x)$ of simple rules and has the form

$t(x)\;{:=}\;r_{1}(x)\wedge r_{2}(x)\wedge\ldots\wedge r_{C}(x),$

or equivalently $t(x)=\prod_{v=1}^{C}r_{v}(x)$. The conjunction $t(x)$ checks
if a subset of components of $x$ have certain desired values. We also refer to
each conjunction as a rule. For a fixed $C$, there are $\binom{2d}{C}$ such
rules. A _disjunctive normal form_ (DNF) classifier has the form

$h(x)\;{:=}\;t_{1}(x)\vee t_{2}(x)\vee\ldots\vee t_{M}(x).$

The DNF term $h(x)$ is often associated with the set of rules
$\\{t_{1},\ldots,t_{M}\\}$, and thus is also called a rule set. The $h(x)$
assigns a label of $1$ to $x$ if it satisfies any of the conjunction terms
$t_{m}$. This form allows the classifier to capture non-linear relationships
between the features $x$ and the label $y$. As an example, consider the loan
application default risk classification problem from FICO, (2018), where the
label ‘high risk’ may be indicated by the disjunction of the two conjunction
terms

(#Loans $\geq 7$) $\vee$ ( (#Loans $\leq 5)\,\wedge\,$(Total Amount
$\geq\$10,000)\,)$.

To score the explanatory power of a rule set, we associate with a conjunction
$t(x)$ of $C$ simple rules a cost $c(t)=C+1$. The cost of a DNF term $h$ is
$c(h)=\sum_{m=1}^{M}c(t_{m})=M+\sum_{m=1}^{M}C_{m}$.

The DNF model predicts a label $\hat{y}=h(x)$ at data point $x$. Dash et al.,
(2018) construct _sparse_ rule sets with a constraint on their complexity
$c(h)\leq{\mathfrak{C}}$. Let ${\mathcal{H}}({\mathfrak{C}})$ denote the space
of all DNF rule sets that obey this constraint. Then, for a given
${\mathfrak{C}}$, the rule set is constituted from a possible
$\sum_{c=1}^{{\mathfrak{C}}-1}\binom{2d}{c}$ conjunction rules.

A convex combination of rule sets takes the form
$F(\cdot)=\sum_{k}v_{k}h_{k}(\cdot)$ where each $h_{k}(x)$ is a rule set and
the non-negative $v_{k}$ satisfy $\sum_{k}v_{k}=1$. It predicts labels by
applying a threshold: $\hat{y}=\mathbb{I}(F(x)\geq 1/2)$. The cost associated
with $F$ is the sum of the individual costs $c(F)=\sum_{k}c(h_{k})$.

###  2.2 Analyzing Rule Sets

The standard approach to finding a good learning model is to minimize the
expected loss under the empirical pmf $P^{0}$ over the training dataset
${\mathcal{S}}$. Sparse rule sets are constructed by the optimization
formulation:
$\min_{h\in{\mathcal{H}}({\mathfrak{C}})}\mathbb{E}_{P^{0}}\,l(h(x),y)$, where
the loss is $l(h(x),y)=|h(x)-y|$ to reflect the impact of $\hat{y}\neq y$. The
following result sheds light on the expected gap in performance of any model
on the (unknown) true data distribution $P^{t}$ and $P^{0}$.

######  Proposition 1.

Let $P^{t}$ be the underlying distribution and $P^{0}$ the empirical pmf over
a training set ${\mathcal{S}}$ of size $N$, and the loss is measured as
$l(h(x),y)=|h(x)-y|$. Then, for any $\delta\in(0,1)$, there exists constants
such that for all rule sets $h$ with max-complexity $c(h)\leq{\mathfrak{C}}$,
we have with probability at least $1-\delta$:

$\mathbb{E}_{P^{t}}[l(h(x),y)]\;\leq\;{\mathbb{E}_{P^{0}}[l(h(x),y)]}\;+\;\bigg{(}\frac{2}{N}\bigg{(}({\mathfrak{C}}-1)\log\frac{2d}{{\mathfrak{C}}-1}+\log\frac{1}{\delta}\bigg{)}\bigg{)}^{\frac{1}{2}}.$
(1)

Note that, for binary output $\hat{y}=h(x)$,
$\mathbb{E}_{P}l(h(x),y)=\mathbb{E}_{P}|\hat{y}-y|=\mathbb{P}_{P}[\hat{y}\neq
y]$. Detailed expressions on the dependence of the rate on $\delta$ are
provided in the proof presented in Section A of the supplement, which extends
the Vapnik-Chervonenkis bounds from Golea et al., (1997) and Schapire et al.,
(1998). The misclassification losses under the true distribution $P^{t}$ and
the empirical distribution $P^{0}$ represent the test and training performance
of a model, respectively. The typical learning approach chooses an optimal
model that minimizes the training error, the first term in the upper bound (1)
on the test error. However, the second term is independent of the particular
model chosen and grows as $O(\sqrt{{\mathfrak{C}}})$ for a given dataset,
i.e., with fixed $N$. This explains the phenomenon observed in the results of
Dash et al., (2018) that optimal rule sets satisfying a higher complexity
${\mathfrak{C}}$ constraint can reduce the training error to zero but do not
generalize as well as lower complexity solutions.

In learning settings where the models are parametrized with variables that
take values in a continuous space, the standard approach to alleviate this
generalization issue is to regularize the parameter value of the chosen model
in an $l_{p}$ norm to enforce some form of sparsity. Rule sets lack such a
regularizable parameter, but sparsity can be directly enforced by the right
choice of ${\mathfrak{C}}$. Dash et al., (2018) employ cross-validation to
achieve this, but at the expense of a large computational cost since the
procedure involves repeatedly solving a complex combinatorial minimization
problem.

###  2.3 DRO for Learning

DRO has emerged recently (Namkoong and Duchi,, 2017; Blanchet et al., 2019a, ;
Lam and Zhou,, 2017) as an alternative new statistical model training
technique with rigorous foundations to establish its ability to produce models
with improved generalization. Moreover, Ghosh et al., (2021) show that these
formulations can be solved to produce models of similar or higher
generalization quality on a significantly lower computational budget.

More specifically, DRO methods have been developed to handle the case where
models are parametrized with $\theta$, a continuous-valued variable. Let
$l(\theta,\xi)$ provide a loss value for the task at hand (classification,
regression, and so on) with model parameter $\theta$ at data $\xi$. As before,
we are given an $N$-sized training dataset ${\mathcal{S}}$ of $\xi$ values.
The DRO approach selects the model parameter that minimizes the _robust loss_
$R(\theta)$:

$\displaystyle\min_{\theta}$
$\displaystyle\;\quad\left\\{\;R(\theta)\,\,{:=}\,\,\max_{P\in{\mathcal{P}}}\;\mathbb{E}_{\xi\sim
P}l(\theta,\xi)\;\;\right\\}$ where
$\displaystyle\;\;{\mathcal{P}}{:=}\\{P\,\,|\,\,D_{\phi}(P,P^{0})\leq\rho\,,\,\,\sum_{i}P_{i}=1,\,\,P_{i}\geq
0\\}.$

The robust loss $R(\theta)$ is the worst loss observed for the model
parameterized with $\theta$ over a set ${\mathcal{P}}$ of probability vectors
representing weights over the training dataset ${\mathcal{S}}$. This set
${\mathcal{P}}$ is centered at the empirical pmf $P^{0}$ and allows all
weights that are within a $\phi$-divergence distance $D_{\phi}(\cdot,\cdot)$
of $\rho$ from $P^{0}$. The $\phi$-divergence distance function on the space
of weights is defined as
$D_{\phi}(P^{\prime},P){:=}\mathbb{E}_{P}\left[\phi\left(\frac{dP^{\prime}}{dP}\right)\right]$,
where $\phi$ is positive, convex and $\phi(1)=0$. Members of this class of
distance functions include the modified $\chi^{2}$ divergence, with
$\phi(s)=(s-1)^{2}$, and the Kullback-Leibler (KL) divergence, with
$\phi(s)=s\log s-s+1$.

Lam and Zhou, (2017) show that, for well chosen $\rho$, the set
$\\{\arg\min_{\theta}\,\mathbb{E}_{P}l(\theta,\xi),\,\,P\in{\mathcal{P}}\\}$
contains the minimizer of the expectation under the true data distribution
$P^{t}$ with high probability, even if $P^{t}$ itself does not belong to
${\mathcal{P}}$. Namkoong and Duchi, (2017) (Section 3.1) consider the DRO
formulation with a loss convex in $\theta$ and a $\chi^{2}$-divergence
uncertainty ball with radius $\rho=O({\mathfrak{V}}/N)$, where
${\mathfrak{V}}$ is the (finite) _Vapnik-Chervonenkis_ (VC) dimension of the
class of models parametrized by $\theta$. They show that, with high
probability, the robust loss $R(\theta)$ is bounded above by

$R(\theta)\;\;\leq\;\;\mathbb{E}_{P^{0}}l(\theta,\xi)\,+\,\left(1+o(1)\right)\sqrt{\frac{{\mathfrak{V}}}{N}\mbox{Var}_{P^{0}}l(\theta,\xi)},$

where $\mbox{Var}_{P^{0}}$ is the variance of the loss function with respect
to the empirical pmf $P^{0}$ and we follow the notation that a function
$g(v)=o(v)$ if $g(v)/v\rightarrow 0$ as $v\rightarrow\infty$. This upper bound
is a close first order approximation of the right hand side of (1). (The VC
dimension of the DNF rule set class of complexity up to ${\mathfrak{C}}$ in
(1) is ${\mathfrak{C}}$.) Thus, the model that minimizes the robust loss
stands to generalize better over unseen data, given the relation (1), and does
not require extensive fine tuning via heuristics like cross-validation. The
experimental evidence over a large set of (convex) classification examples
provided by Ghosh et al., (2021) bears this out in practice.

The development of this DRO theory motivates the consideration of an analogous
formulation for rules sets, that of minimizing the robust loss of a rule set:

$R(h)=\max_{P=(P_{i})\in{\mathcal{P}}}\bigg{\\{}\mathbb{E}_{P}l(h(x),y)\bigg{\\}},$

where, as before, the function $l(h(x),y)=|\hat{y}-y|$ measures the penalty
incurred when $\hat{y}=h(x)\neq y$. This formulation, however, can lead to
undesirable behavior because the misclassification loss $l$ only takes two
values. The solution
$P^{\ast}(h)=\arg\max_{P\in{\mathcal{P}}}\mathbb{E}_{P}l(h(x),y)$
correspondingly has a simple structure, partitioning the dataset
${\mathcal{S}}$ into two subsets of the correctly classified and the
misclassified, each with the same probability assigned within its members.
This is quite uninformative and easily leads to cycling when solving the
minimax formulation in alternating iterations over $h$ and $P$: subsequent
iterations of corresponding experiments where $P$ is updated were quickly
observed to cyclically over- and under-emphasize the same subsets. (Additional
details and results are provided in Section 4.1.)

###  2.4 Ensemble Rule Sets

The issue in not being able to directly adapt DRO formulations developed for
continuously valued parameter models to rule set optimization lies in the
discrete nature of the models. An extensive literature has emerged in the past
two decades on developing stronger models by taking linear combinations of
binary models; these include ensembling models (Sollich and Krogh,, 1995),
such as AdaBoost (Freund and Schapire,, 1997) and their variants, and an
extended treatment is available for collections of rule-set like decision
trees such as random forests (Breiman,, 2001).

We refine the generalization result in Proposition 1 for individual rule sets
to their convex combinations.

######  Proposition 2.

For any convex rule set combination $F=\sum_{k}v_{k}h_{k}$, define
${\mathfrak{C}}^{\prime}\;{:=}\;\max_{k}\\{c(h_{k})\\}$ to be the maximum
complexity of its constituent conjunctive terms, and let the predicted label
be $\hat{y}=\mathbb{I}(F(x)\geq 1/2)$. Then, for any given $\delta\in(0,1)$,
with probability at least $1-\delta$, we have

$\displaystyle\mathbb{P}_{t}(L(F,x,y)\leq 0)\quad$
$\displaystyle\leq\quad\mathbb{P}_{0}(L(F,x,y)\leq\Delta)\quad+\quad
O\bigg{(}\frac{1}{\sqrt{N}}\bigg{(}\frac{({\mathfrak{C}}^{\prime}-1)}{\Delta^{2}}\log
N+\log\frac{1}{\delta}\bigg{)}^{1/2}\bigg{)}.$

The proof of Proposition 2, along with expanded expressions of the rate
displaying its dependence on $\delta$, is given in Section A of the
Supplement. This result shows that convex combinations of rule sets enjoy an
important advantage over a single rule set; the generalization penalty (the
second term) grows only with the largest complexity of its constituents, while
allowing for extra _continuous_ parameters (the $v_{k}$) to minimize the
training error (the first term). Hence, for a fixed complexity budget, convex
combinations of rule sets can significantly improve training performance over
single rule sets while displaying the generalization performance of its
largest constituent rule set.

Importantly, the introduction of the continuous parameters now alleviates the
binary loss issue faced in implementing DRO training of rule sets, since the
_margin_ of misclassification takes (depending on the $v_{k}$) more
informative values. We modify the loss function to highlight this margin, and
for a given datapoint $x$ and its label $y$ set it as

$l(F(x),y)=\bigg{|}\sum_{k}v_{k}h_{k}(x)-\frac{1}{2}\bigg{|}\cdot\mathbb{I}(\hat{y}\neq
y).$ (2)

The algorithm presented in the next section advocates for a DRO approach of
training convex combinations of rule sets with the loss function (2).

##  3 ENSEMBLING ALGORITHM

We now present our algorithm for identifying the optimal sparse convex DNF
combination $F(x)$. Its task is divided into two optimization formulations
that are solved in sequence. The first identifies a collection of good DNF
models using a distributionally robust combinatorial formulation. The DNF
models considered are constrained to be of complexity less than or equal to
${\mathfrak{C}}^{\prime}$. The second formulation subsequently selects an
optimal sparse convex combination of DNF models from the identified
collection. This formulation enforces that convex combinations of only
$\lfloor{\mathfrak{C}}/{\mathfrak{C}}^{\prime}\rfloor$ DNF models are
considered, where ${\mathfrak{C}}$ is the user-desired sparsity parameter.

###  3.1 DRO-Based Collection of Rule Sets

The aim of this DRO formulation is to select a collection
$F(x)=\sum_{n}(1/n)\cdot h_{n}(x)$ of DNF models $h_{n}(\cdot)$ that minimizes
the robust loss objective:

$R(F)\,=\,\max_{P\in{\mathcal{P}}}\,\left\\{\,\mathbb{E}_{P}l(F(x),y)\right\\},$
(3)

where the probability ball ${\mathcal{P}}$ is centered at $P^{0}$ and admits
all probability mass functions $P=(P_{i})$ that are a distance
$D_{\phi}(P,P^{0})$ of $\rho>0$ from $P^{0}$. The loss objective is as given
in (2). Note that the collection gives equal weight to its constituent DNF
models.

The formal distributionally robust approach to constructing an optimal
collection consists of solving for $\min_{F}R(F)$ over all collections of any
arbitrary size. The feasible space of this formulation allows for a countably
infinite collection of DNF models, and is hard to search over. A similar
problem has been addressed in the ensembling literature using efficient
approximations. For instance, Friedman et al., (2000) show that the AdaBoost
algorithm for growing ensembles of models can be viewed as a greedy approach
to solving a generalized loss minimization problem over the set of all
ensembles for a specific global loss objective. Moreover, Blanchet et al.,
2019b  show that this generalized loss objective for AdaBoost can itself be
written as a robust objective function of the form (3).

This motivates the following sequential approach to growing the collection
using a greedy addition of one new DNF member in each iteration. We start with
the initial weights $P^{0}$ as the empirical (equal weights) pmf and set the
initial collection of DNF models to be empty, i.e., $F_{0}(x)=0$. Then, in
each iteration indexed by $n=1,\ldots$, we have the following steps.

  * 1.

Solve for

$h_{n}\,{:=}\,\arg\min_{h\in{\mathcal{H}}({\mathfrak{C}}^{\prime})}\,\,\mathbb{E}_{P^{n-1}}l(h(x),y).$
(4)

  * 2.

Add $h_{n}$ to collection: space$\quad
F_{n}(x)\leftarrow\frac{n-1}{n}F_{n-1}(x)+\frac{1}{n}h_{n}(x)$.

  * 3.

Solve for the robust loss maximizing pmf

$P^{n}\,\,{:=}\,\,\arg\max_{P\in{\mathcal{P}}}\;\mathbb{E}_{P}\,l(F_{n}(x),y).$
(5)

In each iteration, a new DNF model is selected that minimizes the
misclassification loss over ${\mathcal{S}}$ as weighted by the worst-case pmf
$P^{n-1}$. The DNF loss function of $l(h)=|h(x)-y|$ is used in (4). This DNF
model is then added to the collection. The last step then refreshes the worst
pmf by maximizing for the robust loss of the current collection $F_{n}(x)$,
using the form (2).

Column Generation for (4): This formulation searches over the (exponentially
large) space ${\mathcal{H}}({\mathfrak{C}}^{\prime})$ of rule sets. We follow
the approach of Dash et al., (2018) to efficiently explore this space using a
column generation reformulation, and extend their approach to the general case
where the data weights are any arbitrary pmf $P^{n-1}$. Let $K$ be the set of
all possible rules. For each $k\in K$, let $Z_{k}$ be the set of indices of
data points with label 0. For each data index $i$, let $K_{i}$ be the set of
indices $k$ such that rule $k$ is satisfied by data point $x_{i}$. Let $c_{k}$
be the complexity of rule $k$. We then solve the formulation

$\displaystyle\min\quad$ $\displaystyle\sum_{i}\xi_{i}+\sum_{k\in K}\sum_{i\in
Z_{k}}P_{i}w_{k}$ (6) $\displaystyle s.t.\quad$
$\displaystyle\xi_{i}+P_{i}\sum_{k\in K_{i}}w_{k}\geq P_{i},\;\;\xi_{i}\geq
0,$ (7) $\displaystyle\sum_{k\in K}c_{k}w_{k}\leq C,\;w_{k}\in\\{0,1\\},\;k\in
K.$ (8)

The formulation chooses a set of rules (given by nonzero values of $w_{k}$)
that does not minimize the $0$-$1$ loss but instead a weighted Hamming loss:
for each data point $x_{i}$ with label $1$, $\xi_{i}$ provides the exact loss
weighted by $P_{i}$, while the second objective term with $P_{i}$ values set
to $1$ provides the number of rules in the chosen rule set that give incorrect
values to $0$ label data points. Since the number of rules is exponential in
number, we solve this formulation via column generation as in (Dash et al.,,
2018) except that terms in the dual constraints for the LP relaxation of (6) –
(8) are scaled by the values $P_{i}$. Section B in the Supplement provides
further details on the column generation method.

Solving for worst case pmf in (5): The objective function of (5) is expressed
as $\sum_{i\in{\mathcal{S}}}l(F_{n}(x_{i}),y_{i})\cdot P_{i}$, where the
individual loss values $z_{i}=l(F_{n}(x_{i}),y_{i})$ are constructed as per
(2). The $z_{i}$ take values in the discrete set $\\{0,1/n,\ldots,1\\}$, and
so as the iterations $n$ grow, the issue noted in Section 2.3 is progressively
alleviated. The formulation (5) can be equivalently written as the concave
problem:

$\displaystyle{R}(F_{n})$
$\displaystyle\;=\;\max_{{P}=({P}_{i})\in{\mathcal{P}}}\sum_{i\in{\mathcal{S}}}{P}_{i}z_{i}$
(9) s.t. $\displaystyle\sum_{i\in{\mathcal{S}}}\phi(N{P}_{i})\leq
N\rho,\;\;\sum_{i\in{\mathcal{S}}}{P}_{i}=1,\,{P}_{i}\geq 0.$

Here, the first constraint represents the restriction on the $\phi$-divergence
$D_{\phi}(P,P^{0})$ being within $\rho$. The objective of (9) is linear in the
variable $P$ and the constraints are convex, thus rendering a concave program.
Ghosh et al., (2021) utilize Lagrangian methods to provide an algorithm that
can solve this formulation efficiently. The details of this formulation are
provided in Section C of the supplement. We close this description by adapting
a result from Ghosh et al., (2021) for a worst-case bound on the computational
effort required to obtain an $\epsilon$-optimal solution to (9).

######  Proposition 3 ((Ghosh et al.,, 2021, Proposition 2)).

For any $\phi$-divergence, the algorithm in Section C finds a feasible
solution $(\tilde{P}^{*})$ to (9) with an objective function value $\tilde{R}$
such that $|R(F)-\tilde{R}|<\epsilon$ with a worst-case computational effort
bounded by $O(N\log N+(\log\frac{1}{\epsilon})^{2})$, where $\epsilon$ is a
small precision parameter.

###  3.2 Selecting Sparse Convex Ensembles

The formulation in the previous section sequentially generates a collection
$\\{h_{k}(x),\,\,k=1,\ldots,n\\}$ of DNF models each of which are the optimal
solutions to minimizing the expected loss under a sequence of worst-case pmf
values $P^{n}$. Each $h_{k}$ has been further chosen such that their maximum
complexity is at most ${\mathfrak{C}}^{\prime}$. We now address the question
of selecting an optimal sparse convex combination (or weighted ensemble) from
within this collection.

As defined in Section 2.4, the misclassification loss (2) of the convex
combinations measures the margin of misclassification of $\hat{y}$. The
optimal sparse convex combination $\sum_{k=1}^{n}v_{k}h_{k}$ is then
determined by solving

$\min_{v_{k}}\;\mathbb{E}_{P}\,l\bigg{(}\sum_{k=1}^{n}v_{k}h_{k}(x),y\bigg{)}\qquad\mbox{s.t.}\sum_{\\{k\,|\,v_{k}>0\\}}c(h_{k})\leq{{\mathfrak{C}}}.$
(10)

This formulation is solved after every iteration of the collection-growing
formulation is finished. The data weight vector used in (10) can be $P^{n}$,
the worst-case pmf of the current iteration. The performance of this
formulation does not seem to be sensitive to this choice of $P^{n}$ in our
experiments of Section 4; and a selection of $P=P^{0}$, the initial empirical
distribution, produces similar results. (Additional details and discussion of
this issue are provided in Sections C and D of the supplement.) The complexity
parameter ${\mathfrak{C}}$ represents the user desired complexity of the
overall classifier output by our method. Section 4 further discusses how
${\mathfrak{C}}^{\prime}$ should be set, and also addresses the case where
both parameters ${\mathfrak{C}}$ and ${\mathfrak{C}}^{\prime}$ can be set
adaptively.

We turn our attention to solving (10) efficiently. Let ${\mathcal{S}}^{0}$ and
${\mathcal{S}}^{1}$ be the subsets of the training data indices that represent
data labeled with $y=0$ and $y=1$, respectively; so,
${\mathcal{S}}^{0}\cup{\mathcal{S}}^{1}={\mathcal{S}}$. For each data index
$i$, let ${\mathcal{K}}_{i}$ represent the set of all DNF models that provide
a value of $1$ at $x_{i}$:
${\mathcal{K}}_{i}\,=\,\\{k\,|\,h_{k}(x_{i})=1,\,\,k=1,\ldots,n\\}$. Define
auxiliary variables $\xi_{i}$ for each data index $i$ to represent the amount
by which the data is misclassified. We then rewrite the problem (10) as the
IP:

$\displaystyle\min_{\xi^{0},\xi^{1}}\qquad\quad\quad\sum_{i\in{\mathcal{S}}^{0}}p_{i}\xi_{i}\,\,+\,\,\sum_{i\in{\mathcal{S}}^{1}}p_{i}\xi_{i}$
(11) $\displaystyle\mbox{subject
to:}\quad\sum_{k\in{\mathcal{K}}_{i}}v_{k}\,\,\leq\,\,\frac{1}{2}-\Delta+\xi_{i},\;\;\;\forall
i\in{\mathcal{S}}^{0},$ (12)
$\displaystyle\quad\xi_{i}+\sum_{k\in{\mathcal{K}}_{i}}v_{k}\,\,\geq\,\,\frac{1}{2}+\Delta,\;\;\;\quad\;\forall
i\in{\mathcal{S}}^{1},$ (13) $\displaystyle\quad
v_{k}\in[0,1],\;\quad\sum_{k}v_{k}=1,\;\;\;\;\xi^{0}_{i},\xi^{1}_{i}\geq 0,$
$\displaystyle\quad v_{k}\leq w_{k},\;\quad
w_{k}\in\\{0,1\\},\;\;\sum_{k}w_{k}c(h_{k})\;\leq\;{\mathfrak{C}}.$ (14)

The constraints (12) and (13) define the magnitude of the misclassification
margin for each data index $i$ using the auxiliary variables $\xi_{i}$,
depending on by how much the threshold of $1/2$ is breached for each label
type. Note that an extra parameter $\Delta$ is introduced to strictly separate
the two data subsets. Our experiments in Section 4 show that a small
separation parameter, e.g.,
$\Delta={\mathfrak{C}}/\left(2\cdot\min_{k}c(h_{k})\right)$, suffices. The
binary variables $w_{k}$ are introduced in (14) to enforce the maximum allowed
complexity of the chosen convex combination weight values $v_{k}\in[0,1]$,
$\sum_{k}v_{k}=1$.

The formulation (11) is a relatively small IP with only as many binary
variables as the number of DNF models $h_{k}$ in the current collection. In
Section 4, we present an adaptive procedure for implementing our algorithm
that solves this IP formulation a few times after each iteration of the
collection-growing formulation of Section 3.1 in order to adaptively identify
the best sparse convex combination to return.

##  4 NUMERICAL EXPERIMENTS

We follow the numerical experiments conducted by Dash et al., (2018) in
evaluating their column-generation based sparsifier and conduct tests on seven
classification datasets from the UCI repository (Lichman,, 2013). In addition,
we use the dataset from the FICO Explainable Machine Learning Challenge
(FICO,, 2018), which contains $23$ numerical features of the credit history of
$10,459$ individuals ($9,871$ after removing records with all entries missing)
for predicting repayment risk (good/bad). Section D of the supplement provides
additional details on the various datasets and how they were pre-processed.

Table 1: (Middle-Left Four Columns:) Classification Performance over Unseen Test Data of Models Created by Four Competing Methods for Seven UCI Repository Datasets and the FICO Dataset. (Middle-Right Three Columns:) Complexity of Models Constructed by the Three Explainable Methods. (Last Column:) Run-Length (number of iterations $n$) of DR Method. All Results presented with a Mean and a $95\%$ Confidence Interval. Dataset | Test Performance (%) | Model Complexity |  $n$, Number
---|---|---|---
Name | DR | CG | CART | RF | DR | CG | CART | Iterates
heart |  81.5(1.7) | 78.9(4.7) |  81.6(4.7) |  82.5(1.4) | 26.0(2.7) |  11.3(3.5) | 32.0(15.9) | 35.7(3.6)
ILPD |  69.1(1.1) |  69.6(2.4) | 67.4(3.1) |  69.8(1.0) |  13.3(2.3) |  10.9(5.3) | 56.5(21.4) | 27.0(1.8)
FICO |  72.1(0.3) | 71.7(1.0) | 70.9(0.6) |  73.1(0.2) | 25.0(2.5) |  13.3(8.0) | 155.0(53.9) | 33.9(3.3)
ionosphere |  93.3(1.4) | 90.0(3.5) | 87.2(3.5) |  93.6(1.4) | 20.3(2.9) |  12.3(5.9) | 46.1(8.2) | 38.9(4.0)
liver |  58.6(1.6) |  59.7(4.7) | 55.9(2.7) |  60.0(1.6) | 29.3(1.4) |  5.2(2.4) | 60.2(30.6) | 38.5(3.2)
pima |  74.6(1.4) | 74.1(3.7) | 72.1(2.5) |  76.1(1.6) | 25.0(2.8) |  4.5(2.5) | 34.7(11.4) | 34.5(3.6)
transfusion |  78.4(0.8) |  77.9(2.7) |  78.7(2.2) |  77.3(0.6) | 14.3(1.3) |  5.6(2.4) | 14.3(4.5) | 30.2(3.1)
WDBC |  95.5(0.5) | 94.0(2.4) | 93.3(1.8) |  97.2(0.4) | 28.3(1.9) |  13.9(4.7) | 15.6(4.3) | 36.8(3.9)

Our procedure to identify the best sparse convex combinations of DNF models
implements our above algorithm by adaptively determining the best values for
the complexity parameters ${\mathfrak{C}}^{\prime}$ (from (4)) and
${\mathfrak{C}}$ (from (10)). The generalization result in Proposition 2,
bolstered by our expanded set of results in Section 4.2, shows that setting a
low value for ${\mathfrak{C}}^{\prime}$ grants more generalization flexibility
to the sparse ensembling formulation (10), and so we set
${\mathfrak{C}}^{\prime}=5$ throughout the results of this section. The two
formulations are then solved sequentially in iterations $n=1,\ldots$ as
follows.

  1. 1.

Solve the DRO rule set formulation (4) using the column-generation approach to
identify an additional member to be added to the collection
$\\{h_{k},\,k=1,\ldots,n\\}$ of good DNF models; each model has complexity
$c(h_{k})\leq{\mathfrak{C}}^{\prime}$. Then solve the maximization formulation
(5) to obtain the pmf $P^{n}$ that maximizes robust loss.

  2. 2.

Solve the sparse convex combination formulation (10) in the IP form (11) –
(14) for increasing values of ${\mathfrak{C}}=m\cdot{\mathfrak{C}}^{\prime}$
with $m=2,\ldots$. Stop when the value of the identified optimal solution does
not decrease beyond a given threshold.

The inner iterations in Step 2 are cutoff when ${\mathfrak{C}}=30$. The outer
iterations $n$ are stopped when the optimal solution returned by the second
step does not improve significantly over $20$ successive iterations. We set
the improvements threshold to show $0.5\%$ or more change. The dataset weights
used in (10) are set to the empirical (equal weights) pmf $P^{0}$.

Implementation Details. We implement the procedure described above in Java,
and use version $22.1$ of the commercial solver IBM ILOG CPLEX to solve the
linear programs and mixed-integer linear programs. Our column generation
implementation is an extension of the approach described in Dash et al.,
(2018) and we refer to that paper for more details. The main parameters of the
implementation are a time limit on the solution of the pricing problem in each
iteration, which we set to $30$ seconds, a second time limit on the overall
column generation process, which we set to $300$ seconds, and a limit on the
number of column generation iterations, which we set to $5$. We set the DRO
parameter $\rho$ to $0.05$ as measured using $\chi^{2}$-divergence for the
distance function $D_{\phi}(P,P^{0})$. The value of the parameter $\rho$ is
chosen to match the recommendation of $O({\mathfrak{C}}/N)$ from Namkoong and
Duchi, (2017) for continuous convex learning settings. The implementation of
the sparse convex combination IP formulation (11) – (14) limits the time spent
in the solver to $600$ seconds.

The results presented here are based on splitting each input dataset uniformly
at random into a training dataset, containing $70\%$ of the data, and a test
dataset with the remaining $30\%$ of the data. The implementation then uses
the training dataset to execute the sequence of our two formulations and
obtain the set of final rules, which is followed by evaluating the performance
of those rules on the test dataset and reporting the results. A total of $20$
such permutations are considered to produce the confidence intervals
presented.

Competing Methods. We compare the test performance of our DRO-based rule set
induction approach (DR) with the sparse rule set algorithm (CG) proposed by
Dash et al., (2018). In addition, we include results for two decision tree
models presented in (Dash et al.,, 2018): the CART method (Breiman et al.,,
1984) and the random forest RF method (Breiman,, 2001). The latter method
produces dense ensembles of decision trees and is not appropriate from the
explainability point of view, but often produces models with good
generalization for comparison.

Table 1 provides both the test accuracy and the complexity of the models
identified by the various methods. For both these sets of results, the method
with the best average outcomes among interpretable models is highlighted in
bold and the best average outcome overall is underlined. Any method with
average outcomes within the confidence interval of the overall best outcome is
highlighted in italics.

The results in Table 1 (middle-left) show that DR comes close to and in half
of the eight cases even exceeds the generalization performance of the
competing explainable models CG and CART. Moreover, the solution quality of DR
is statistically comparable to the non-interpretable RF method, which produces
the best results in all but one dataset (transfusion). In fact, across all the
datasets, the mean performance of the DR method lies within the confidence
intervals of the overall best method, thus either representing the best method
or statistically identical to the best method.

We also compare the complexity of the models from our DR approach with those
obtained by three of the methods; the omitted RF method produces a dense
ensemble with orders of magnitude higher complexity and is not appropriate for
creating interpretable models. The results in Table 1 (middle-right) show that
the complexity of the estimated models from DR is within a reasonable range of
the CG method, even statistically identical in one case (ILPD), while
consistently better than CART and often significantly so. These results
illustrate the efficacy predicted in Proposition 2, in that sparse convex
mixtures of ${\mathfrak{C}}^{\prime}=5$ rule sets have greater flexibility in
reducing training error and simultaneously keep a low generalization gap. In
contrast, the individual rule sets produced by the CG method use higher
complexity than $5$ to obtain a reasonable balance between reduced training
error and lower test efficacy, but overall do not perform as well on test data
as the convex combinations generated by DR. While the CG method may improve
its training power with a higher ${\mathfrak{C}}$ bound, it may suffer a
higher generalization gap as shown in Proposition 1. Hence, it is not readily
clear that the gap in test performance can be closed for any value of
${\mathfrak{C}}$ by the CG method.

The last column also shows how many iterations of the DR method were required
to attain the corresponding results. Note that each iteration contains one
call to solve the expensive formulation (4). These results compare very
favourably with the order of hundreds of calls to similar model estimation
formulations required by the cross-validation based fine tuning of parameters
of the competing methods. For example, CG solves $150$ of a similar
combinatorial optimization formulation with the aid of their column generation
approach.

Additional results are provided below in Section 4.2 and in the supplement
where the parameters chosen in this section are varied. A larger sparse
mixture complexity cutoff of $60$ is considered and, as can be expected, this
slightly improves the generalization performance at the expense of more
complex (and harder to interpret) models. Results for a larger DRO parameter
$\rho=0.5$ are presented, largely illustrating that this value makes the DRO
formulation too conservative. Parameters of the stopping criteria are also
varied to show that with a higher computational budget, slightly higher
quality solutions can be obtained by DR for the same complexity cutoff of
$30$.

###  4.1 Cyclical Behavior of DRO for DNF Models

We noted in Section 2.3 that basing the performance on the misclassification
losses $l(h(x),y)=|h(x)-y|$ can lead to undesirable cyclical behavior. This is
because $l$ takes only two values $\\{0,1\\}$, and thus the solution
$P^{\ast}(h)=\arg\max_{P\in{\mathcal{P}}}\mathbb{E}_{P}l(h(x),y)$ is
correspondingly degenerate, partitioning the dataset ${\mathcal{S}}$ into two
subsets of the correctly classified and the misclassified. The data points in
each of the two subsets are assigned the same probability. The subset over
which $l(h(x),y)=1$ has a large value (since all points in this subset are
misclassified) and the points where the loss is zero have a small, possibly
zero, probability (since all points in this subset are correctly classified).
This is quite uninformative and easily leads to cycling when solving the
minimax formulation in alternating iterations over $h$ and $P$, in that
subsequent iterations in which $P$ is updated can cause the same two subsets
of points to exchange loss values of $0$ and $1$, and thus the same subsets
are over- and under-emphasized in a cyclical manner. This phenomenon is
demonstrated in Figure 1, which presents results from solving the DRO
formulation of minimizing the robust loss objective $R(h)$ defined in Section
2.3. A simplified version of the heuristic presented in Section 3.1 was
implemented to solve this formulation: Step 2 of growing the collection of DNF
models was omitted and in Step 3, we solve for the worst-case solution
$P^{\ast}(h_{n})$ that maximizes the weighted expected loss of the current DNF
iterate $h_{n}$ only.

Figure 1: Histogram over Training Dataset of Four Consecutive Iterations
($\\#47-\\#50$) of the DRO Formulation in Section 2.3 with $\\{0,1\\}$ Loss of
DNF Models. Horizontal axis is over Loss values in $[0,1]$ and Vertical Axis
provides Histogram Counts of datapoints. Experiments were over the ‘Heart’
instance; DRO Constraints of $\rho=0.05$ Measured by $\chi^{2}$-distances were
Imposed, and all $h_{n}\in{\mathcal{H}}$ satisfy ${\mathfrak{C}}\leq 5$.

Figure 1 plots in orange histograms of the instantaneous losses
$l(h_{n}(x_{i}),y_{i})$ observed for the $n$-th iterate $h_{n}$ over the
datapoints (indexed by $i$) in the training dataset for $n=47,48,49$ and $50$.
The experiments were conducted for the ‘heart’ problem instance, and DNF
models with complexity ${\mathfrak{C}}\leq 5$ were considered. Since losses
$l$ are binary, this gives us a two-bin histogram in each iteration, with bins
at loss values $0$ and $1$. Additionally provided in blue are histograms of
the observed running average loss values at the data points:
$(1/n)\sum_{j=1}^{n}l(h_{j}(x_{i}),y_{i})$. The histograms of the average
values are instructive in that they show identical histograms in every other
iteration, which indicates that the chosen DNF iterates $h_{n}$ flip the
instantaneous loss values of a handful of datapoints and leaves the vast
majority unchanged. This leads to a roughly equal split of the dataset into
three subsets, the first where no chosen $h_{n}$ iterate can classify the
datapoints correctly (and so has an average loss of $1$), the second where all
chosen $h_{n}$ iterates classify the points correctly (with an average loss of
$0$), and the third set where alternate iterates flip the correct/incorrect
classification (leading to an average loss of $\approx 0.5$). These results
illustrate why the DRO philosophy of training continuously valued parameter
models is of limited direct use when the models can produce only binary
$\\{0,1\\}$ loss information.

###  4.2 Comparisons of Parameter Settings

This section expands on the performance of the DR algorithm. We present a
number of tables documenting the results of numerical experiments where
specific parameters of the algorithm are varied. In each table presented
below, we present two sets of columns each for the classification performance
over unseen test data, the complexity of models constructed and the run length
(number of iterations $n$) of the DR method. The first column is the results
of the DR method with hyperparameters as presented in Table 1, and the results
of the second column are obtained by changing one hyperparameter at a time.
All Results are presented with a Mean and a $95\%$ Confidence Interval, and
the result with the best mean performance value among the two columns are
highlighted in bold.

Table 2: Sparse Ensembles Constructed from Low (${\mathfrak{C}}^{\prime}=5$) Complexity DNF Models are Compared to Ensembles Constructed from Moderate (${\mathfrak{C}}^{\prime}=10$) Complexity DNF Models.  Dataset | Test Performance (%) of DR | Model Complexity of DR |  $n$, number of Iterates of DR
---|---|---|---
Name | ${\mathfrak{C}}^{\prime}=5$ | ${\mathfrak{C}}^{\prime}=10$ | ${\mathfrak{C}}^{\prime}=5$ | ${\mathfrak{C}}^{\prime}=10$ | ${\mathfrak{C}}^{\prime}=5$ | ${\mathfrak{C}}^{\prime}=10$
heart |  81.5(1.7) | 79.9(1.7) |  26.0(2.7) | 30.0(0.0) |  35.7(3.6) | 37.1(3.4)
ILPD | 69.1(1.1) |  69.8(0.8) |  13.3(2.3) | 30.0(0.0) |  27.0(1.8) | 33.9(3.5)
FICO |  72.1(0.3) | 71.5(0.3) |  25.0(2.5) | 30.0(0.0) | 33.9(3.3) |  28.4(2.0)
ionosphere |  93.3(1.4) | 92.8(0.9) |  20.3(2.9) | 30.0(0.0) |  38.9(4.0) | 39.6(2.3)
liver |  58.6(1.6) | 55.6(2.0) |  29.3(1.4) | 30.0(0.0) | 38.5(3.2) |  34.4(2.8)
pima |  74.6(1.4) | 72.9(0.9) |  25.0(2.8) | 30.0(0.0) |  34.5(3.6) | 35.9(3.8)
transfusion |  78.4(0.8) | 77.5(0.8) |  14.3(1.3) | 29.0(1.3) | 30.2(3.1) |  27.6(0.8)
WDBC | 95.5(0.5) |  95.6(0.5) |  28.3(1.9) | 30.0(0.0) |  36.8(3.9) | 36.9(2.7)

The results presented in Table 2 contrast the ensembling of DNF models with a
low value for ${\mathfrak{C}}^{\prime}=5$ to ensembling models with double the
complexity ${\mathfrak{C}}^{\prime}=10$. All other parameters are held to
their value used in Table 1: $P=P^{0}$ in the ensembling formulation (11),
${\mathfrak{C}}$ is cutoff at $30$, the DRO formulation (5) sets $\rho=0.05$,
and the outer iterations stop when no progress is made in $20$ consecutive
iterations. From the results, it is clear that the low complexity DNF models
grant more generalization flexibility to the sparse ensembling formulation
(10), beating the best ensembles assembled from ${\mathfrak{C}}^{\prime}=10$
in six of the eight cases (and significantly so, for example the ‘liver’
instance) and is within the $95\%$ confidence interval of the best in all
cases. Moreover, the total model complexity ${\mathfrak{C}}$ produced for
${\mathfrak{C}}^{\prime}=5$ ensembles are significantly smaller in all
instances. The generalization guarantee provided by Proposition 2 anticipates
this, since it predicts that performance may degrade as
${\mathfrak{C}}^{\prime}$ increases.

Table 3: Aggregate Complexity of Sparse Ensemble Determined by (11) with a Cutoff at ${\mathfrak{C}}\leq 30$ Compared to a Cutoff at ${\mathfrak{C}}\leq 60$ Dataset | Test Performance (%) of DR | Model Complexity of DR |  $n$, number of Iterates of DR
---|---|---|---
Name | cutoff $=30$ | cutoff $=60$ | cutoff $=30$ | cutoff $=60$ | cutoff $=30$ | cutoff $=60$
heart | 81.5(1.7) |  82.2(1.7) |  26.0(2.7) | 36.0(5.0) |  35.7(3.6) | 38.7(3.7)
ILPD |  69.1(1.1) |  69.1(1.1) |  13.3(2.3) |  13.3(2.3) |  27.0(1.8) |  27.0(1.8)
FICO | 72.1(0.3) |  72.2(0.3) |  25.0(2.5) | 33.8(5.0) |  33.9(3.3) | 38.1(3.7)
ionosphere |  93.3(1.4) |  93.3(1.4) |  20.3(2.9) |  20.3(2.9) |  38.9(4.0) |  38.9(4.0)
liver | 58.6(1.6) |  59.3(1.3) |  29.3(1.4) | 47.0(5.) |  38.5(3.2) | 42.1(3.4)
pima |  74.6(1.4) | 74.5(1.4) |  25.0(2.8) | 30.3(3.7) |  34.5(3.6) | 35.5(3.6)
transfusion |  78.4(0.8) |  78.4(0.8) |  14.3(1.3) |  14.3(1.3) |  30.2(3.1) |  30.2(3.1)
WDBC | 95.5(0.5) |  95.7(0.5) |  28.3(1.9) | 33.8(4.0) |  36.8(3.9) | 38.2(3.6)

Table 3 varies the cutoff allowed in the inner sparse ensembling iterations
where ${\mathfrak{C}}$ is increased as
${\mathfrak{C}}=m\cdot{\mathfrak{C}}^{\prime},\,\,m=1,2,\ldots$. The results
presented in Table 1 above cuts off at ${\mathfrak{C}}\leq 30$ or earlier,
which is compared to a larger cutoff of $60$. Other parameters are held the
same as in Table 1: ${\mathfrak{C}}^{\prime}=5$, the formulation (11) uses
$P=P^{0}$, the DRO formulation (5) sets $\rho=0.05$, and outer iterations stop
when $20$ consecutive iterations do not improve the training performance. As
can be expected, the higher cutoff slightly improves the generalization
performance, but this is at the expense of significantly more complex (and
harder to interpret) models. Overall the cutoff of $30$ appears to be
preferable from the view of all three factors of good generalization
performance at a low computational cost while producing models with good
interpretability.

Table 4 presents another illustration of the better generalization of the
ensembles constructed by the DR method than the DNF models produced by the CG
method. Test performance of the CG and DR methods are provided where the total
complexity budget ${\mathfrak{C}}$ is kept fixed to the value determined by
the DR method when run with the settings that produced the results in Table 1.
In every case, we observe that the test performance produced by the DR method
exceeds (by a considerable margin in some cases) that produced by the CG
method. Proposition 1 states that the DNF models produced by the CG method may
suffer generalization error that is larger than the training error by a
$O(\sqrt{{\mathfrak{C}}})$ term. So, while the DNF models have lower training
error with the higher ${\mathfrak{C}}$ budget, they will generalize less well.
This is also attested to by the best model complexity chosen for CG with
cross-validation as presented in Table 1: the average complexity is smaller
than that of the DR output in nearly every case (though statistically
identical for the ILPD dataset) and each is indeed closer to the complexity
${\mathfrak{C}}^{\prime}=5$ of the models that are aggregated by the DR
method. On the other hand, by virtue of the generalization bound in
Proposition 2, an ensemble of models produced by DR can reduce training error
while still enjoying test performance within an
$O(\sqrt{{\mathfrak{C}}^{\prime}})$ term.

Table 4: Results for CG when Complexity Observed in DR Results are Imposed as Constraints on CG Dataset | Test Performance (%)
---|---
Name | DR | CG
heart |  81.5(1.7) | 77.4(1.4)
ILPD |  69.1(1.1) | 68.7(1.4)
FICO |  72.1(0.3) | 70.9(0.3)
ionosphere |  93.3(1.4) | 88.4(1.2)
liver |  58.6(1.6) | 55.1(1.7)
pima |  74.6(1.4) | 71.7(1.3)
transfusion |  78.4(0.8) | 77.5(0.7)
WDBC |  95.5(0.5) | 93.2(0.8)

## Acknowledgements

We thank Alex Gray for suggesting the problem and the potential benefits of
combining DRO and column generation to address the problem.

## References

  * (1) Blanchet, J., Kang, Y., and Murthy, K. (2019a).  Robust wasserstein profile inference and applications to machine learning.  Journal of Applied Probability, 56(3):830–857.
  * (2) Blanchet, J. H., Kang, Y., Zhang, F., and Hu, Z. (2019b).  A distributionally robust boosting algorithm.  2019 Winter Simulation Conference (WSC), pages 3728–3739.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine Learning, 45:5–32.
  * Breiman et al., (1984) Breiman, L., Friedman, J. H., Olshen, R. A., and Stone, C. J. (1984).  Classification and regression trees.  Biometrics, 40:874.
  * Dash et al., (2018) Dash, S., Gunluk, O., and Wei, D. (2018).  Boolean decision rules via column generation.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc.
  * Dhurandhar et al., (2020) Dhurandhar, A., Shanmugam, K., and Luss, R. (2020).  Enhancing simple models by exploiting what they already know.  In Proceedings of the 37th International Conference on Machine Learning, ICML’20. JMLR.org.
  * FICO, (2018) FICO (2018).  FICO explainable machine learning challenge.  https://community.fico.com/community/xml.  Last accessed 2018-05-16.
  * Freund and Schapire, (1997) Freund, Y. and Schapire, R. E. (1997).  A decision-theoretic generalization of on-line learning and an application to boosting.  Journal of Computer and System Sciences, 55(1):119–139.
  * Friedman et al., (2000) Friedman, J., Hastie, T., and Tibshirani, R. (2000).  Additive logistic regression: a statistical view of boosting (With discussion and a rejoinder by the authors).  The Annals of Statistics, 28(2):337 – 407.
  * Friedman and Popescu, (2008) Friedman, J. H. and Popescu, B. E. (2008).  Predictive learning via rule ensembles.  Annals of Applied Statistics, 2(3):916–954.
  * Ghosh et al., (2021) Ghosh, S., Squillante, M. S., and Wollega, E. (2021).  Efficient generalization with distributionally robust learning.  In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W., editors, Advances in Neural Information Processing Systems, volume 34, pages 28310–28322. Curran Associates, Inc.
  * Golea et al., (1997) Golea, M., Bartlett, P., Lee, W. S., and Mason, L. (1997).  Generalization in decision trees and dnf: Does size matter?  In Jordan, M., Kearns, M., and Solla, S., editors, Advances in Neural Information Processing Systems, volume 10. MIT Press.
  * Hu et al., (2019) Hu, X., Rudin, C., and Seltzer, M. (2019).  Optimal sparse decision trees.  In Advances in Neural Information Processing Systems (NeurIPS).
  * Lakkaraju et al., (2016) Lakkaraju, H., Bach, S. H., and Leskovec, J. (2016).  Interpretable decision sets: A joint framework for description and prediction.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, page 1675–1684.
  * Lam and Zhou, (2017) Lam, H. and Zhou, E. (2017).  The empirical likelihood approach to quantifying uncertainty in sample average approximation.  Operations Research Letters, 45(4):301–307.
  * Lichman, (2013) Lichman, M. (2013).  UCI machine learning repository.
  * Luenberger, (1969) Luenberger, D. G. (1969).  Optimization by Vector Space Methods.  John Wiley & Sons.
  * Muggleton et al., (2018) Muggleton, S. H., Schmid, U., Zeller, C., Tamaddoni-Nezhad, A., and Besold, T. (2018).  Ultra-strong machine learning: comprehensibility of programs learned with ilp.  Machine Learning, 107:1119–1140.
  * Namkoong and Duchi, (2017) Namkoong, H. and Duchi, J. C. (2017).  Variance-based regularization with convex objectives.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems 30, pages 2971–2980. Curran Associates, Inc.
  * Ribeiro et al., (2016) Ribeiro, M. T., Singh, S., and Guestrin, C. (2016).  ”why should i trust you?”: Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, page 1135–1144.
  * Rudin, (2019) Rudin, C. (2019).  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  Nat Mach Intell., 1(5):206––215.
  * Rudin et al., (2022) Rudin, C., Chen, C., Chen, Z., Huang, H., Semenova, L., and Zhong, C. (2022).  Interpretable machine learning: Fundamental principles and 10 grand challenges.  Statist. Surv., 16:1–85.
  * Schapire et al., (1998) Schapire, R. E., Freund, Y., Bartlett, P., and Lee, W. S. (1998).  Boosting the margin: A new explanation for the effectiveness of voting methods.  The Annals of Statistics, 26(5):1651–1686.
  * Schmid et al., (2017) Schmid, U., Zeller, C., Besold, T., Tamaddoni-Nezhad, A., and Muggleton, S. H. (2017).  How does predicate invention affect human comprehensibility?  In ILP 2016.
  * Sollich and Krogh, (1995) Sollich, P. and Krogh, A. (1995).  Learning with ensembles: How overfitting can be useful.  In Touretzky, D., Mozer, M., and Hasselmo, M., editors, Advances in Neural Information Processing Systems, volume 8. MIT Press.
  * Wang et al., (2017) Wang, T., Rudin, C., Doshi-Velez, F., Liu, Y., Klampfl, E., and MacNeille, P. (2017).  A bayesian framework for learning rule sets for interpretable classification.  Journal of Machine Learning Research, pages 1–37.

##  Appendix A PROOFS OF THEORETICAL ANALYSIS

The DNF models $h\in{\mathcal{H}}$ produce binary output $h\,:\;x\rightarrow
y\in\\{0,1\\}$, where ${\mathcal{H}}$ denotes the space of all DNF rule sets
of interest. The data $(x,y)$ follow the true (unknown) distribution $P^{t}$,
and we are given only an $N$-sized sample dataset
${\mathcal{S}}=\\{(x_{n},y_{n}),\,\,n\in[N]\\}$ generated from $P^{t}$, using
the notation $[N]=\\{1,\ldots,N\\}$. The empirical equal-weight probability
mass function (pmf) on ${\mathcal{S}}$ is denoted by $P^{0}$. In the sequel,
we use the shorthand $\mathbb{E}_{t}$ and $\mathbb{P}_{t}$ for
$\mathbb{E}_{P^{t}}$ and $\mathbb{P}_{P^{t}}$, respectively, and similarly
$\mathbb{E}_{0}$ and $\mathbb{P}_{0}$ for $\mathbb{E}_{P^{0}}$ and
$\mathbb{P}_{P^{0}}$.

Size of ${\mathcal{H}}$. The finite set
${\mathcal{H}}={\mathcal{H}}({\mathfrak{C}})$ contains all the DNF terms of
complexity at most ${\mathfrak{C}}$. There are $\binom{2d}{c}$ choices for
conjunctions of size $c$, and the complexity cost of each such term is $c+1$.
Hence, each $h\in{\mathcal{H}}$ is constituted from conjunctions of up to size
${\mathfrak{C}}-1$.

For the number of conjunctions of size $c$, we have the lower bound

$\binom{2d}{c}=\frac{2d}{c}\cdot\frac{2d-1}{c-1}\ldots\frac{2d-c+1}{1}\geq\bigg{(}\frac{2d}{c}\bigg{)}^{c}.$

Meanwhile, we obtain the corresponding upper bound

$\binom{2d}{c}=\frac{2d}{c}\cdot\frac{2d-1}{c-1}\ldots\frac{2d-c+1}{1}\leq\frac{(2d)^{c}}{c!}\leq\bigg{(}\frac{2de}{c}\bigg{)}^{c},$

where the last inequality uses the fact that
$\exp(c)=\sum_{j=0}^{\infty}c^{j}/j!$.

Since the number of conjunctions grows exponentially with $c$, there exists a
constant $V>1$ such that $|{\mathcal{H}}|\leq
V\cdot\binom{2d}{{\mathfrak{C}}-1}$. To see this, consider the subset of
${\mathcal{H}}$ with one term of size $c<{\mathfrak{C}}-1$. Then, this subset
is of size at most
$\binom{2d}{c}\cdot\binom{2d}{{\mathfrak{C}}-1-(c+1)}\lll\binom{2d}{{\mathfrak{C}}-1}$.
Thus, we also have that

$({\mathfrak{C}}-1)\log\frac{2d}{{\mathfrak{C}}-1}\;\leq\;\log|{\mathcal{H}}|\;\leq\;\log
V+({\mathfrak{C}}-1)\log\frac{2de}{{\mathfrak{C}}-1}.$ (15)

###  A.1 Proposition 1: Generalization of DNF Models

We now establish the following proposition from the main body of the paper.

Proposition 1. Let $P^{t}$ be the underlying distribution and $P^{0}$ the
empirical pmf over a training set ${\mathcal{S}}$ of size $N$, and the loss is
measured as $l(h(x),y)=|h(x)-y|$. Then, for any $\delta\in(0,1)$, there exists
constants such that for all rule sets $h$ with max-complexity
$c(h)\leq{\mathfrak{C}}$, we have with probability at least $1-\delta$:

$\mathbb{E}_{P^{t}}[l(h(x),y)]\;\leq\;{\mathbb{E}_{P^{0}}[l(h(x),y)]}\;+\;\bigg{(}\frac{2}{N}\bigg{(}({\mathfrak{C}}-1)\log\frac{2d}{{\mathfrak{C}}-1}+\log\frac{1}{\delta}\bigg{)}\bigg{)}^{\frac{1}{2}}.$

Proof. For the true (unknown) data distribution $(x,y)\sim P^{t}$ over the
space of features and labels, let $p_{h}=\mathbb{P}_{t}(h(x)\neq y)$. Since
both $h$ and $y$ are binary, the loss function
$l(h(x),y)=|h(x)-y|=\mathbb{I}(h(x)\neq y)$. Define the random variable
$Z_{h}=\mathbb{I}(h(x)\neq y)$ as a Bernoulli random variable with the
probability of success (misclassification) $p_{h}$. The training dataset
${\mathcal{S}}$ can also be viewed as a set of $N$ independent draws of
Bernoulli random variables $Z_{n}$, and
$\hat{p}_{h,N}=\frac{1}{N}\sum_{n\in[N]}\mathbb{I}(Z_{n}=1)$ is an $N$ sample
approximation of the probability $p_{h}$.

For each $h\in{\mathcal{H}}$, the right-tail Chernoff bound for Bernoulli
random variables yields

$\mathbb{P}(p_{h}-\hat{p}_{h,N}\geq\lambda)\leq\exp(-N\lambda^{2}/2),$

where the probability is taken with respect to the sampling of ${\mathcal{S}}$
from $P_{t}$. To obtain a uniform bound that applies to _all_
$h\in{\mathcal{H}}$, we employ the union bound principle that
$\mathbb{P}(\cup_{i}A_{i})\leq\sum_{i}\mathbb{P}(A_{i})$ for events $A_{i}$.
Thus, we have

$\mathbb{P}\left(\exists
h\in{\mathcal{H}}\;|\;p_{h}-\hat{p}_{h,N}\geq\lambda\right)\leq|{\mathcal{H}}|\exp(-N\lambda^{2}/2).$

Using the lower bound in (15), we obtain the desired generalization result
when $\lambda$ is set as

$\lambda=\bigg{(}\frac{2}{N}\bigg{(}({\mathfrak{C}}-1)\log\frac{2d}{{\mathfrak{C}}-1}+\log\frac{1}{\delta}\bigg{)}\bigg{)}^{\frac{1}{2}}.$
∎

###  A.2 Proposition 2: Generalization of Convex Combinations

We consider convex combination models of the form
$F(x)=\sum_{h\in{\mathcal{H}}}v_{h}h(x)$, where the real-valued $v_{h}$ are
non-negative and $\sum_{h}v_{h}=1$. The set representing misclassification by
$F(x)$ is compactly written as

$\bigg{\\{}\bigg{(}y-\frac{1}{2}\bigg{)}\bigg{(}F(x)-\frac{1}{2}\bigg{)}\leq
0\bigg{\\}}\;\Longleftrightarrow\;\bigg{\\{}F(x)\leq\frac{1}{2}\,\,,\,\,y=1\bigg{\\}}\,\,\bigcup\,\,\bigg{\\{}F(x)>\frac{1}{2}\,\,,\,\,y=0\bigg{\\}}.$

Define $L(f,x,y){:=}(y-1/2)(f(x)-1/2)$. We wish to obtain a bound on
$\mathbb{P}_{t}(L(F,x,y)\leq 0)$ in terms of the corresponding probability
under the empirical pmf $P^{0}$. The next result provides a bound on the
generalization performance of such convex models. Note that, since the $v_{h}$
are real-valued, any uniform bound on the worst performance of _any_ $F$
cannot be stated as cleanly as in Proposition 1. We introduce an additional
parameter $\Delta>0$ in the following proposition from the main body of the
paper, which we now establish.

Proposition 2. For any convex rule set combination $F=\sum_{k}v_{k}h_{k}$,
define ${\mathfrak{C}}^{\prime}\;{:=}\;\max_{k}\\{c(h_{k})\\}$ to be the
maximum complexity of its constituent conjunctive terms, and let the predicted
label be $\hat{y}=\mathbb{I}(F(x)\geq 1/2)$. Then, for any given
$\delta\in(0,1)$, with probability at least $1-\delta$, we have

$\displaystyle\mathbb{P}_{t}(L(F,x,y)\leq 0)\quad$
$\displaystyle\leq\quad\mathbb{P}_{0}(L(F,x,y)\leq\Delta)\quad+\quad
O\bigg{(}\frac{1}{\sqrt{N}}\bigg{(}\frac{({\mathfrak{C}}^{\prime}-1)}{\Delta^{2}}\log
N+\log\frac{1}{\delta}\bigg{)}^{1/2}\bigg{)}.$

The parameter $\Delta$ plays an important role in the analysis of
generalization as presented in Proposition 2. The upper bound described above
becomes increasingly lax as $\Delta\searrow 0$. It is not merely a convenient
device to ease the construction of the rigorous guarantee, but it also has an
important practical implication in the sparse ensembling formulation (11). Our
experiments show that the $\Delta$ in this formulation, which plays the same
role as in Proposition 2, must be a small positive value for the optimal
sparse convex combinations identified by the formulation to be of higher
generalization power than the individual DNF models.

Proof of Proposition 2.

Let us start by noting that the set of convex weights
$\\{v_{h},\,h\in{\mathcal{H}}={\mathcal{H}}({\mathfrak{C}}^{\prime})\\}$
themselves represent a pmf. The analysis of the generalization error of $F$
utilizes this fact and approximates $F$ using an $M$-sized collection of DNFs
sampled from $(v_{h})$, as introduced by Schapire et al., (1998). Define the
approximation set

${\mathcal{G}}_{M}:=\bigg{\\{}g\;:x\rightarrow\bigg{\\{}0,\frac{1}{M},\ldots,1\bigg{\\}}\;\bigg{|}\;g(x)=\frac{1}{M}\sum_{m\in[M]}h_{m}(x),\;\mbox{where}\,\,\,h_{m}\in{\mathcal{H}}({\mathfrak{C}}^{\prime})\,\,\,\forall\,m\in[M]\bigg{\\}}.$

(We will specify the right size $M$ to use later in the proof.) Note that
$|{\mathcal{G}}_{M}|=|{\mathcal{H}}|^{M}$. The analysis of
$g\in{\mathcal{G}}_{M}$ is easier than $F$ directly because $g$ only takes
values in the lattice $\\{0,1/M,2/M,\ldots,1\\}$. Let ${\mathcal{V}}$ be the
distribution induced on ${\mathcal{G}}_{M}$ by sampling the components of a
$g\in{\mathcal{G}}_{M}$ from the pmf formed from the coefficients $(v_{h})$ of
$F$. Our main result is established in two steps. In the first step, we start
by uniformly bounding the probability of the misclassification event
$\\{L(F,x,y)\leq 0\\}$ of $F$ under the true distribution $P^{t}$ with the
probability that $L(g,x,y)\leq\Delta/2$ under the empirical pmf $P^{0}$ for an
_arbitrarily_ chosen $g\in{\mathcal{G}}_{M}$ that is sampled from the
distribution ${\mathcal{V}}$. We then bound the latter with the probability
that $L(F,x,y)\leq\Delta$ under the same empirical pmf $P^{0}$.

Let $A=\\{L(F,x,y)\leq 0\\}$ and $B=\\{L(g,x,y)\leq\Delta/2\\}$, where $g$ is
any (fixed) element from ${\mathcal{G}}_{M}$ and $\Delta>0$ is an arbitrary
constant. Then, by the law of total probability and the definition of
conditional probability, we have that
$\mathbb{P}_{t}A\,\,=\,\,\mathbb{P}_{t}[A\cap B]\,\,+\,\,\mathbb{P}_{t}[A\cap
B^{c}]=\mathbb{P}_{t}[A|B]\,\mathbb{P}_{t}B\,\,+\,\,\mathbb{P}_{t}[B^{c}|A]\,\mathbb{P}_{t}A\,\,\leq\,\,\mathbb{P}_{t}B\,\,+\,\,\mathbb{P}_{t}[B^{c}|A]$.
Thus, we conclude

$\mathbb{P}_{t}(L(F,x,y)\leq
0)\;\leq\;\mathbb{P}_{t}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}\;+\;\mathbb{P}_{t}\bigg{(}L(g,x,y)>\frac{\Delta}{2}\,\,\bigg{|}\,\,L(F,x,y)\leq
0\bigg{)}.$

Since this applies to any $g\in{\mathcal{G}}_{M}$, we can take expectations on
the right hand side with respect to the distribution ${\mathcal{V}}$ on
${\mathcal{G}}_{M}$ to obtain that

$\displaystyle\mathbb{P}_{t}(L(F,x,y)\leq 0)$
$\displaystyle\leq\mathbb{E}_{{\mathcal{V}}}\bigg{[}\mathbb{P}_{t}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}\bigg{]}+\mathbb{E}_{t}\bigg{[}\mathbb{P}_{{\mathcal{V}}}\bigg{(}L(g,x,y)>\frac{\Delta}{2}\,\,\bigg{|}\,\,L(F,x,y)\leq
0\bigg{)}\bigg{]}.$ (16)

We shall bound in turn the inner probability terms within the expectations of
each summand on the right hand side of (16). For the first summand, we treat
the $g$ as a fixed sum of DNF terms, and thus apply the same union bound
analysis from the proof of Proposition 1 to the Bernoulli random variable
$\mathbb{I}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}$. This yields that the
probability (over the sampling of ${\mathcal{S}}$ from $P^{t}$) that there
exists at least one $g$ and $\Delta>0$ for which

$\mathbb{P}_{t}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}>\mathbb{P}_{0}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}\,\,+\,\,\lambda_{M}$
(17)

is at most

$(M+1)|{\mathcal{H}}|^{M}\exp(-2N\lambda^{2}_{M}).$

The term $\exp(-2N\lambda^{2}_{M})$ comes from the Chernoff bound for _any_
fixed $g$ and $\Delta$. There are $|{\mathcal{G}}_{M}|=|{\mathcal{H}}|^{M}$
such $g$ over which to take the union. Moreover, each $g\in{\mathcal{G}}_{M}$
takes one of only $M+1$ distinct values $\\{0,1/M\ldots,1\\}$, and there are
at most $(M+1)|{\mathcal{H}}|^{M}$ combinations of $g$ and $\Delta$ to
consider. Hence, for any fixed $M$, we proceed as in the proof of Proposition
1 to set

$\displaystyle\lambda_{M}=\bigg{(}\frac{1}{2N}\log\bigg{(}\frac{(M+1)|{\mathcal{H}}|^{M}}{\delta_{M}}\bigg{)}\bigg{)}^{\frac{1}{2}},$
(18)

and thus obtain that, with at least probability $1-\delta_{M}$, the bound in
(17) is satisfied for _every_ $g$ and $\Delta$. Taking the expectation needed
in (16) with respect to the distribution ${\mathcal{V}}$, we have with
probability $1-\delta_{M}$ that

$\mathbb{P}_{t,{\mathcal{V}}}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}\leq\mathbb{P}_{0,{\mathcal{V}}}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}+\lambda_{M}.$

To complete the analysis of this first summand of (16), we need to provide an
upper bound for the first term on the right hand side of the above expression.
Utilizing the same conditioning argument that led to (16), we can write that

$\mathbb{P}_{0,{\mathcal{V}}}\bigg{(}L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}\;\leq\;\mathbb{P}_{0}\bigg{(}L(F,x,y)\leq{\Delta}\bigg{)}\,\,+\,\,\mathbb{E}_{0}\bigg{[}\mathbb{P}_{{\mathcal{V}}}\bigg{(}L(F,x,y)>{\Delta}\,\,\bigg{|}\,\,L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}\bigg{]}.$
(19)

Here, the $F$ in the first term on the right hand side is unaffected by the
distribution ${\mathcal{V}}$ of $g\in{\mathcal{G}}_{M}$. The probability
inside the expectation in the second term is similar in spirit to that of the
second summand in (16) and deals with the sampling of $g\in{\mathcal{G}}_{M}$
from the induced distribution ${\mathcal{V}}$. Note that
$F(x)=\mathbb{E}_{{\mathcal{V}}}g(x)$ and the pmf $(v_{h})$ from which the
$g(x)\sim{\mathcal{V}}$ is sampled has finite support of size
$|{\mathcal{H}}|$. This allows us to exploit the Chernoff bound for finite
support distributions. In particular, for the second term in (19), we obtain
that

$\displaystyle\mathbb{P}_{{\mathcal{V}}}\bigg{(}L(F,x,y)>{\Delta}\,\,\bigg{|}\,\,L(g,x,y)\leq\frac{\Delta}{2}\bigg{)}\;\leq\;\exp\bigg{(}-\frac{M\Delta^{2}}{8}\bigg{)}.$
(20)

For the second summand in (16), we analogously have

$\displaystyle\mathbb{P}_{{\mathcal{V}}}\bigg{(}L(g,x,y)>\frac{\Delta}{2}\,\,\bigg{|}\,\,L(F,x,y)\leq
0\bigg{)}\;\leq\;\exp\bigg{(}-\frac{M\Delta^{2}}{8}\bigg{)}.$ (21)

Upon combining (16), (17), (19), (20) and (21), we have established, for an
arbitrary $M$ and with probability at least $1-\delta_{M}$, that

$\displaystyle\mathbb{P}_{t}(L(F,x,y)\leq
0)\;\;\leq\;\;\mathbb{P}_{0}\bigg{(}L(F,x,y)\leq{\Delta}\bigg{)}\,\,+\,\,2\exp\bigg{(}-\frac{M\Delta^{2}}{8}\bigg{)}+\lambda_{M},$
(22)

where $\lambda_{M}$ is as chosen in (18). To obtain a bound that applies
uniformly for _any_ $M$, we once again appeal to the union bound analysis by
setting $\delta_{M}=\delta\bigg{(}\frac{1}{M}-\frac{1}{M+1}\bigg{)}$. Thus,
with probability at least $1-\sum_{M\geq 1}\delta_{M}=1-\delta$, we have that
the bound in (22) is uniformly satisfied for any $M\geq 1$.

There are two terms that depend on $M$ in the right side of (22). The middle
term exponentially decreases with $M$, while $\lambda_{M}$ (with form given in
(18) after substituting $\delta_{M}$) has a leading term that increases as a
small polynomial of $M$. In order to strike a balance between these terms, we
choose a value for $M$ by setting
$\exp(-M\Delta^{2}/8)=\sqrt{\log|{\mathcal{H}}|/N}$. Thus, we obtain

$M=\frac{4}{\Delta}^{2}\log\bigg{(}\frac{N}{\log|{\mathcal{H}}({\mathfrak{C}}^{\prime})|}\bigg{)}.$

Further, using the lower bound on
$\log|{\mathcal{H}}({\mathfrak{C}}^{\prime})|$ from (15) yields the desired
result. ∎

##  Appendix B COLUMN GENERATION FOR (4)

The size of the set $K$ in the formulation to search over the space of rule
sets ${\mathcal{H}}({\mathfrak{C}}^{\prime})$ clearly grows as an exponential
function of the number of binary features. Solving IPs with exponentially many
variables is impractical; even the associated LP relaxation is hard to solve.
Notice that the number of constraints (that are not variable constraints) in
this formulation is bounded by the number of data points with label 1. One
approach to solving large LPs that have few constraints but very many
variables is to use column generation. This approach is motivated by the
following: (1) there is a solution to the LP in which the number of nonzero
variable values is bounded above by the number of constraints; (2) when
solving the LP relaxation with the simplex method, it is fairly common that
the total number of variables that are ever set to nonzero values in
intermediate steps is a small multiple of the number of constraints. In other
words, only a small subset of all possible $w_{k}$ variables (rules) is
generated explicitly and the optimality of the LP is guaranteed by iteratively
solving a pricing problem that implicitly checks the dual constraints.

To apply this framework to the IP (which we refer to as the “master IP” or
MIP, with associated LP relaxation denoted by “master LP” or MLP), we first
replace $K$ with a very small subset of its elements and explicitly solve the
LP relaxation of the resulting smaller problem, which we call the Restricted
master LP or RMLP. Any optimal solution of RMLP can be extended to a solution
of MLP with the same objective value by setting all missing $w_{k}$ variables
to zero; this solution provides an upper bound on the optimal solution value
of MLP. One can potentially improve such a solution by augmenting RMLP with
additional variables corresponding to some of the missing rules. The second
step is to identify such rules without explicitly considering all of them.
Repeating these steps until there are no improving rules (i.e., variables
missing from the Restricted MLP that can reduce the cost) solves the MLP to
optimality.

To find the missing improving rules, one needs to check for variables in MLP
missing from RMLP that have negative reduced cost (with respect to the dual
solution from RMLP applied to variables in MLP not present in RMLP). The
reduced cost of a missing variable (or variable set to zero in MLP) gives the
maximum possible change in objective value per unit increase in that
variable’s value. Therefore, if all missing variables from RMLP have non-
negative reduced cost, then the optimal solution for RMLP yields an optimal
solution of MLP. Furthermore, missing variables that have large negative
reduced costs are more likely to improve the objective value of RMLP. We next
formulate an optimization problem that uses the optimal dual solution of RMLP
and checks for negative reduced cost columns in MLP. Let $\mu_{i}\geq 0$
denote the dual variable associated with constraint (7) corresponding to the
$i$th data point with label 1. Let $\lambda\geq 0$ be the dual variable
associated with the complexity bound constraint. Let $\delta_{i}\in\\{0,1\\}$
denote whether the $i$th sample satisfies rule $k$. Recall that $c_{k}$ is the
complexity of rule $k$; then its reduced cost is equal to

$\sum_{i:y_{i}=0}P_{i}\delta_{i}-\sum_{i:y_{i}=1}P_{i}\mu_{i}\delta_{i}+\lambda
c_{k}.$ (23)

The first term in (23) is the cost of rule $k$ in the objective function (6),
expressed in terms of $\delta_{i}$. The second term is the sum of the dual
variables associated with constraints (7). The last term is the dual variable
associated with the complexity constraint multiplied by the complexity of rule
$k$.

We now formulate an IP to express rules as conjunctions of the original
features. Let the decision variable $z_{j}\in\\{0,1\\}$ denote if feature $j$
is selected in the rule. Let $S_{i}$ correspond to the zero-valued features in
sample $i$. Let $J$ be the index set of all features. Then the Pricing Problem
below identifies the rule $k$ missing from RMLP that has the lowest reduced
cost.

$\displaystyle z_{CG}~{}=~{}\textbf{min}$
$\displaystyle\lambda\Biggl{(}1+\sum_{j\in
J}z_{j}\Biggr{)}-\sum_{i:y_{i}=1}P_{i}\mu_{i}\delta_{i}$
$\displaystyle+\sum_{i:y_{i}=0}P_{i}\delta_{i}$ (24) s.t.
$\displaystyle\delta_{i}+z_{j}$ $\displaystyle\leq 1,$ $\displaystyle\quad
j\in S_{i},\;\forall i:y_{i}=1$ (25) $\displaystyle\delta_{i}$
$\displaystyle\geq 1-\sum_{j\in S_{i}}z_{j},\quad\delta_{i}\geq 0,\;$
$\displaystyle\forall i:y_{i}=0$ (26) $\displaystyle\sum_{j\in J}z_{j}$
$\displaystyle\leq c,$ (27) $\displaystyle z_{j}$ $\displaystyle\in\\{0,1\\},$
$\displaystyle{j\in J}.$ (28)

The first term in (24) expresses the complexity $c_{k}$ in terms of the number
of selected features. Constraints (25), (26) ensure that the rule acts as a
conjunction, i.e., it is satisfied ($\delta_{i}=1$) only if no zero-valued
features are selected ($z_{j}=0$ for $j\in S_{i}$). Similar to $\xi_{i}$ in
MIP, the variables $\delta_{i}$ do not have to be explicitly defined as binary
due to the objective function. Constraint (27) bounds the number of features
allowed in any rule.

The optimal solution to the Pricing Problem above gives the rule with the
minimum reduced cost that is missing from RMLP. The reduced cost of this rule
equals $z_{CG}$ and if $z_{CG}<0$, then the corresponding variable is added to
RMLP. More generally, any feasible solution to the Pricing Problem that has a
negative objective function value gives a rule with a negative reduced cost
and therefore can be added to RMLP to improve its value.

##  Appendix C SOLVING DRO FORMULATION (5)

We now present an algorithm due to Ghosh et al., (2021) that is adapted in
this paper to solve the DRO formulation (5) for the worst-case pmf. Recall
that the formulation (5) can be equivalently written as the concave problem
given in (9), where the objective function of (5) is expressed as
$\sum_{i\in[N]}l(F_{n}(x_{i}),y_{i})\cdot P_{i}$ and the individual loss
values $z_{i}=l(F_{n}(x_{i}),y_{i})$ are constructed as per (2) taking values
in the discrete set $\\{0,1/n,\ldots,1\\}$. The Lagrangian objective of (9)
can be expressed as

${\mathcal{L}}(\alpha,\lambda,P)=\sum_{i\in{\mathcal{S}}}P_{i}z_{i}+\lambda\left(1-\sum_{i\in{\mathcal{S}}}P_{i}\right)+\frac{\alpha}{N}\left(N\rho-\sum_{i\in{\mathcal{S}}}\Phi(NP_{i})\right).$
(29)

From optimization theory (Luenberger,, 1969), this leads to the optimal
objective value of (9) given by $\tilde{R}^{\ast}(F_{n})=\min_{\alpha\geq
0,\lambda}\max_{P_{i}\geq 0}{\mathcal{L}}(\alpha,\lambda,P)$. We then solve
the Lagrangian formulation in (29) for a given iteration $n$ using Algorithm
1, adapted from (Ghosh et al.,, 2021), to obtain a feasible primal-dual
solution $(P^{*},\alpha^{*},\lambda^{*})$ to (9) with an objective value
$\tilde{R}^{*}$ such that $|R^{*}(F_{n})-\tilde{R}^{*}|<\epsilon$ according to
Proposition 3.

Algorithm 1 Robust Loss Maximization

Given: loss values ${\mathcal{Z}}=\\{z_{1},\ldots,z_{N}\\}$; dataset support
indices $[N]=\\{1,\ldots,N\\}$; $D_{\phi}$ constraint $\rho$.

  1. 1.

Case: $\alpha^{*}=0$ along with constraint $D_{\phi}({P}^{*},P^{0})\leq\rho$.

    1. (a)

Let ${\mathcal{N}}=\\{i\in[N]:z_{i}=\max_{j\in[N]}z_{j}\\}$. Set
$\alpha^{*}=0$ in (29), and then an optimal solution is ${P}^{*}$ where
${P}^{*}_{i}=\frac{1}{|{\mathcal{N}}|},\;\forall i\in{\mathcal{N}}$, and
${P}^{*}_{i}=0,\;\forall i\notin{\mathcal{N}}$.

    2. (b)

If $D_{\phi}({P}^{*},P^{0})\leq\rho$, then stop and return ${P}^{*}$.

  2. 2.

Case: constraint $D_{\phi}({P}^{*},P^{0})=\rho$ with $\alpha^{*}\geq 0$.

    1. (a)

Keeping $\lambda,\alpha$ fixed, solve for the optimal ${P}^{*}$ (as a function
of $\lambda,\alpha$) that maximizes ${\mathcal{L}}(\alpha,\lambda,{P})$,
applying the constraint ${P}_{i}\geq 0$.

    2. (b)

Keeping $\alpha$ fixed, solve for the optimal $\lambda^{*}$ using the first
order optimality condition on ${\mathcal{L}}(\alpha,\lambda,{P}^{*})$. Note
that this is equivalent to satisfying the equation
$\sum_{i\in[N]}{P}^{*}_{i}=1$. The proof of Proposition 3 (Ghosh et al.,,
2021) shows that this step is at worst a bisection search in one dimension,
but in some cases (e.g., KL-divergence) a solution $\lambda^{*}$ is available
in closed form. The proof of Proposition 3 also provides finite bounds
$[\underline{\lambda},\bar{\lambda}]$ on the range over which we need to
search for $\lambda^{*}$.

    3. (c)

Apply the first order optimality condition to the one-dimensional function
${\mathcal{L}}(\alpha,\lambda^{*}(\alpha),{P}^{*})$ to obtain the optimal
$\alpha^{*}\geq 0$. This is equivalent to requiring that $\alpha^{*}$
satisfies the equation $\sum_{i\in[N]}\phi({P}^{*}_{i})=\rho$. The proof of
Proposition 3 (Ghosh et al.,, 2021) shows that this is at worst a one-
dimensional bisection search which embeds the previous step in each function
call of the search.

    4. (d)

Let $\phi^{\prime}(s)=\mathrm{d}\phi(s)/\,\mathrm{d}s$ denote the derivative
of $\phi(s)$ with respect to $s$, and $(\phi^{\prime})^{-1}$ its inverse.
Define the index set ${\mathcal{N}}:=\\{i\in[N]\,\,\mid\,\lambda^{*}\leq
z_{i}-\alpha^{*}\phi^{\prime}(0)\\}$, with ${\mathcal{N}}=\emptyset$ if
$\phi^{\prime}(s)\rightarrow-\infty$ as $s\rightarrow 0+$. Set

$P^{*}_{i}=\left\\{\begin{array}[]{ll}\frac{1}{N}(\phi^{\prime})^{-1}\left(\frac{z_{i}-\lambda^{*}}{\alpha^{*}}\right),&\,\,\,i\in{\mathcal{N}}\\\
0&\,\,\,i\notin{\mathcal{N}}\end{array}\right..$ (30)

Return ${P}^{*}$.

##  Appendix D DETAILS ON NUMERICAL EXPERIMENTS

The datasets were pre-processed as described by Dash et al., (2018). We used
the standard “dummy”/“one-hot” coding to transform categorical features into
multiple binary features. Each category is represented by two binary features,
one indicating the presence of that category and one indicating the absence of
that category. The non-binary numerical features are transformed into a
sequence of binary features where each feature represents the satisfaction or
not of a certain threshold and its negation. As an example, for a certain
numerical feature $X$, we could create the features $X\leq 1$, $X\leq 2$ and
$X>1$, $X>2$. The thresholds used were deciles computed from the samples.

Two of the UCI repository datasets were slightly modified as follows:

  * •

For the ‘liver’ dataset, we used the number of drinks as the output variable
instead of the selector variable. The number of drinks $Y$ was represented by
the two features $Y\leq 2$ and $Y>2$.

  * •

For the ‘heart’ dataset, we used only the Cleveland data and removed 4 samples
with ‘ca’ = ?, yielding 299 samples. The label was binarized as either $X=0$
or $X>0$.

The FICO dataset was also slightly modified due to missing values and the
presence of special values. The 588 records that only contain missing values
(represented by -9) were removed. Entries meaning ‘no inquiries or
delinquencies observed’ (represented by -7) were replaced by the maximum
number of elapsed months in the data plus 1. The remaining missing values
(represented by -9) and entries meaning ‘not applicable’ (represented by -8)
were combined into a single null category. During binarization, a special
indicator was created for these null values and all other comparisons with the
null value return False. Values greater than 7 (other) in
‘MaxDelq2PublicRecLast12M’ were imputed as 7 (current and never delinquent)
based on the corresponding values in ‘MaxDelqEver’.

This section expands on the performance of the DR (Distributionally robust
Rule sets) algorithm. We present a number of tables documenting the results of
numerical experiments where specific parameters of the algorithm are varied.
In each table presented below, we present two sets of columns each for the
classification performance over unseen test data, the complexity of models
constructed and the run length (number of iterations $n$) of the DR method.
The first column is the results of the DR method with hyperparameters as
presented in Table 1, and the results of the second column are obtained by
changing one hyperparameter at a time. All Results are presented with a Mean
and a $95\%$ Confidence Interval, as in the main body of the paper. The result
with the best mean performance value among the two columns are highlighted in
bold.

Table 5: Data Weights $P$ in the Sparse Ensembling Formulation (11) are varied between the Empirical PMF $P^{0}$ and the DRO-identified Worst Weights $P^{n}$ from (5) Dataset | Test Performance (%) of DR | Model Complexity of DR |  $n$, number of Iterates of DR
---|---|---|---
Name | $P^{0}$ | $P^{n}$ | $P^{0}$ | $P^{n}$ | $P^{0}$ | $P^{n}$
heart | 81.5(1.7) |  81.9(1.5) | 26.0(2.7) |  25.3(2.6) |  35.7(3.6) | 37.6(4.0)
ILPD | 69.1(1.1) |  69.6(1.0) |  13.3(2.3) | 15.3(2.4) |  27.0(1.8) | 28.4(2.1)
FICO |  72.1(0.3) | 71.9(0.3) |  25.0(2.5) | 27.8(2.0) |  33.9(3.3) | 36.2(3.3)
ionosphere |  93.3(1.4) | 93.1(1.2) | 20.3(2.9) |  17.3(2.0) | 38.9(4.0) |  38.0(3.4)
liver |  58.6(1.6) |  58.6(1.3) | 29.3(1.4) |  28.3(1.9) |  38.5(3.2) | 39.9(3.2)
pima |  74.6(1.4) | 73.8(1.2) | 25.0(2.8) |  21.3(3.3) |  34.5(3.6) | 37.0(3.7)
transfusion | 78.4(0.8) |  78.6(0.9) |  14.3(1.3) |  14.3(2.5) | 30.2(3.1) |  29.9(2.8)
WDBC |  95.5(0.5) | 95.4(0.6) | 28.3(1.9) |  25.3(2.9) | 36.8(3.9) |  35.4(3.4)

Table 5 analyzes the sensitivity of the sparse ensembling formulation (11) to
the data weights $P$ used. We specifically compare the setting of these
weights to the empirical pmf ($P=P^{0}$) with the setting of these weights to
the current worst-case pmf ($P=P^{n}$) identified by the $n$-th iteration of
the DRO formulation (5). All other parameters are held to their value used in
Table 1: ${\mathfrak{C}}^{\prime}=5$, ${\mathfrak{C}}$ is cutoff at $30$, the
DRO formulation (5) sets $\rho=0.05$, and the stopping criterion of the outer
iteration waits to ensure no improvements are noticed in $20$ consecutive
iterations. The results establish that neither setting dominates the other,
and so in our main results we use $P=P^{0}$. The $P^{n}$ may not bring
additional benefits in the ensembling formulation (11) because the impact of
the worst-case weights are already encoded in the choice of the DNF iterates
produced by the formulation (4) that determines the members of the collection
of iterates among which the sparse ensemble is chosen.

Table 6: DRO Radius Parameter $\rho$ is varied from $0.05$ to $0.5$ Dataset | Test Performance (%) of DR | Model Complexity of DR |  $n$, number of Iterates of DR
---|---|---|---
Name | $\rho=0.05$ | $\rho=0.5$ | $\rho=0.05$ | $\rho=0.5$ | $\rho=0.05$ | $\rho=0.5$
heart | 81.5(1.7) |  81.6(1.1) |  26.0(2.7) | 29.5(1.0) |  35.7(3.6) | 40.3(3.5)
ILPD | 69.1(1.1) |  69.3(1.1) | 13.3(2.3) |  12.8(1.8) | 27.0(1.8) |  26.3(0.3)
FICO |  72.1(0.3) | 71.4(0.3) | 25.0(2.5) |  20.8(3.0) |  33.9(3.3) | 35.3(3.2)
ionosphere |  93.3(1.4) | 91.5(1.1) |  20.3(2.9) | 30.0(0.0) |  38.9(4.0) | 42.2(3.0)
pima |  74.6(1.4) | 73.8(1.5) | 25.0(2.8) |  24.3(3.5) |  34.5(3.6) | 38.7(3.4)
transfusion | 78.4(0.8) |  78.6(0.7) | 14.3(1.3) |  13.0(1.3) | 30.2(3.1) |  26.9(0.8)
WDBC |  95.5(0.5) | 95.0(0.4) | 28.3(1.9) |  27.8(2.0) | 36.8(3.9) |  36.3(3.7)

In Table 6, results are presented for a larger DRO parameter $\rho=0.5$,
keeping all other parameters the same as in Table 1:
${\mathfrak{C}}^{\prime}=5$, ${\mathfrak{C}}$ is cutoff at $30$, the
ensembling formulation (11) uses $P=P^{0}$, and iterations $n$ stop after $20$
consecutive iterations of little improvement. Recall that this parameter was
set to match the recommendation of ${\mathfrak{C}}^{\prime}/N$ from the DRO
analysis of Namkoong and Duchi, (2017) for continuous convex models. The
results presented in Table 6 illustrate that $\rho=0.05$ is on balance better
than the larger $\rho=0.5$ which may make the DRO formulation too
conservative. (Note that, for the ‘liver’ instance, the setting of $\rho=0.5$
was too large rendering overly conservative results, and thus such results are
omitted.)

Table 7: Minimum Number of Consecutive Iterations of No-Improvement to Observe before Stopping Outer Iterations Indexed by $n$ Dataset | Test Performance (%) of DR | Model Complexity of DR |  $n$, number of Iterates of DR
---|---|---|---
Name | $10$ | $20$ | $30$ | $10$ | $20$ | $30$ | $10$ | $20$ | $30$
heart | 81.5(1.7) | 81.5(1.7) |  81.6(1.7) |  25.0(3.0) | 26.0(2.7) | 25.5(2.8) |  21.2(2.0) | 35.7(3.6) | 44.2(2.3)
ILPD | 69.1(1.1) | 69.1(1.1) |  69.2(1.1) |  13.3(2.3) |  13.3(2.3) |  13.3(2.3) |  16.2(0.3) | 27.0(1.8) | 36.8(1.3)
FICO | 72.1(0.2) | 72.1(0.3) |  72.2(0.2) |  25.0(2.9) |  25.0(2.5) | 26.5(2.1) |  22.5(2.5) | 33.9(3.3) | 43.4(2.5)
ionosphere | 93.1(1.4) |  93.3(1.4) | 93.2(1.4) |  19.3(2.5) | 20.3(2.9) | 21.3(2.9) |  22.4(2.5) | 38.9(4.0) | 46.0(2.2)
liver |  58.7(1.4) | 58.6(1.6) | 58.2(1.8) |  28.3(1.9) | 29.3(1.4) | 28.8(1.7) |  26.4(2.4) | 38.5(3.2) | 46.8(1.7)
pima |  74.6(1.2) |  74.6(1.4) | 74.3(1.3) |  23.8(2.8) | 25.0(2.8) | 24.3(2.9) |  22.2(2.4) | 34.5(3.6) | 43.4(2.5)
transfusion |  78.5(0.7) | 78.4(0.8) | 78.4(0.8) |  13.3(1.3) | 14.3(1.3) | 14.3(1.3) |  17.4(2.1) | 30.2(3.1) | 39.6(2.5)
WDBC | 95.4(0.6) | 95.5(0.5) |  95.6(0.5) |  27.3(2.4) | 28.3(1.9) |  27.3(2.4) |  21.3(1.6) | 36.8(3.9) | 44.7(2.3)

Table 7 varies the stopping criterion parameter of how many consecutive
iterations of less than $0.5\%$ aggregate improvement (in training
performance) to observe before terminating the outer iterations. We present
results for this parameter being set at $10$ and $30$ and reproduce the
results for $20$ as presented in the main body of the paper. The other
parameters are held at their previous values: ${\mathfrak{C}}^{\prime}=5$,
$\rho=0.05$, and ${\mathfrak{C}}$ is cutoff at $30$. The more aggressive
parameter value of $10$ can produce lower complexity models more quickly,
while the more conservative setting of $30$ can provide improvements in test
performance at the expense of longer computation time and less interpretable
models. The value of $20$, chosen for the results presented in the main body,
seems to attain a good balance over these three competing factors.