Supersparse linear integer models for optimized medical scoring systems

Abstract

Scoring systems are linear classification models that only require users to add, subtract and multiply a few small numbers in order to make a prediction. These models are in widespread use by the medical community, but are difficult to learn from data because they need to be accurate and sparse, have coprime integer coefficients, and satisfy multiple operational constraints. We present a new method for creating data-driven scoring systems called a Supersparse Linear Integer Model (SLIM). SLIM scoring systems are built by using an integer programming problem that directly encodes measures of accuracy (the 0–1 loss) and sparsity (the $\ell _0$-seminorm) while restricting coefficients to coprime integers. SLIM can seamlessly incorporate a wide range of operational constraints related to accuracy and sparsity, and can produce acceptable models without parameter tuning because of the direct control provided over these quantities. We provide bounds on the testing and training accuracy of SLIM scoring systems, and present a new data reduction technique that can improve scalability by eliminating a portion of the training data beforehand. Our paper includes results from a collaboration with the Massachusetts General Hospital Sleep Laboratory, where SLIM is being used to create a highly tailored scoring system for sleep apnea screening.

1 Introduction

Scoring systems are linear classification models that only require users to add, subtract and multiply a few small numbers in order to make a prediction. These models are used to assess the risk of numerous serious medical conditions since they allow physicians to make quick predictions, without extensive training, and without the use of a computer (see e.g., Knaus et al. 1991; Bone et al. 1992; Moreno et al. 2005). Many medical scoring systems that are currently in use were hand-crafted by physicians, whereby a panel of experts simply agreed on a model (see e.g., the CHADS$_2$ score of Gage et al. 2001). Some medical scoring systems are data-driven in the sense that they were created by rounding logistic regression coefficients (see e.g., the SAPS II score of Le Gall et al. 1993). Despite the widespread use of medical scoring systems, there has been little to no work that has focused on machine learning methods to learn these models from data.

Scoring systems are difficult to create using traditional machine learning methods because they need to be accurate, sparse, and use small coprime integer coefficients. This task is especially challenging in medical applications because models also need to satisfy explicit constraints on operational quantities such as the false positive rate or the number of features. Such requirements represent serious challenges for machine learning. Current methods for sparse linear classification such as the lasso (Tibshirani 1996) and elastic net (Zou and Hastie 2005) control the accuracy and sparsity of models via approximations to speed up computation, and require rounding to yield models with coprime integer coefficients. Approximations such as convex surrogate loss functions, ${\ell _{1}}$-regularization, and rounding not only degrade predictive performance but make it difficult to address operational constraints imposed by physicians. To train a model that satisfies a hard constraint on the false positive rate, for instance, we must explicitly calculate its value, which is impossible when we control accuracy by means of a surrogate loss function. In practice, traditional methods can only address multiple operational constraints through a tuning process that involves high-dimensional grid search. As we show, this approach often fails to produce a model that satisfies operational constraints, let alone a scoring system that is optimized for predictive accuracy.

In this paper, we present a new method to create data-driven scoring systems called a Supersparse Linear Integer Model (SLIM). SLIM is an integer program that optimizes direct measures of accuracy (the 0–1 loss) and sparsity (the $\ell _0$-seminorm) while restricting coefficients to a small set of coprime integers. In comparison to current methods for sparse linear classification, SLIM can produce scoring systems that are fully optimized for accuracy and sparsity, and that satisfy a wide range of complicated operational constraints without any parameter tuning.

The main contributions of our paper are as follows.

• We present a principled machine learning approach to learn scoring systems from data. This approach can produce scoring systems that satisfy multiple operational constraints without any parameter tuning. Further, it has a unique advantage for imbalanced classification problems, where constraints on class-based accuracy can be explicitly enforced.

• We derive new bounds on the accuracy of discrete linear classification models. In particular, we present discretization bounds that guarantee that we will not lose training accuracy when the size of the coefficient set is sufficiently large. In addition, we present generalization bounds that relate the size of the coefficient set to a uniform guarantee on testing accuracy.

• We develop a novel data reduction technique that can improve the scalability of supervised classification algorithms by removing a portion of the training data beforehand. Further, we show how data reduction can be applied directly to SLIM.

• We present results from a collaboration with the Massachusetts General Hospital (MGH) Sleep Laboratory where SLIM is being used to create a highly tailored scoring system for sleep apnea screening. Screening for sleep apnea is important: the condition is difficult to diagnose, has significant costs, and affects over 12 million people in the United States alone (Kapur 2010).

• We provide a detailed experimental comparison between SLIM and eight popular classification methods on publicly available datasets. Our results suggest that SLIM can produce scoring systems that are accurate and sparse in a matter of minutes.

The remainder of our paper is structured as follows. In the remainder of Sect. 1, we discuss related work. In Sect. 2, we introduce SLIM and discuss its special properties. In Sect. 3, we explain how SLIM can easily accommodate operational constraints that are important for medical scoring systems to be used in practice. In Sect. 4, we present theoretical bounds on the accuracy of SLIM scoring systems and other discrete linear classification models. In Sect. 5, we present a data reduction technique to decrease the computation associated with training SLIM scoring systems and other supervised classification models. In Sect. 6, we discuss a collaboration with the MGH Sleep Laboratory where we used SLIM to create a highly tailored scoring system for sleep apnea screening. In Sect. 7, we report experimental results to show that SLIM can create scoring systems that are accurate and sparse in minutes. In Sect. 8, we present two specialized extensions of SLIM.

1.1 Related work

Our work is related to several streams of research, namely: medical scoring systems; sparse linear classification; discrete linear classification; and mixed-integer programming (MIP) approaches for classification. In what follows, we discuss related work in each area separately.

1.1.1 Medical scoring systems

Some popular medical scoring systems include:

• SAPS I, II and III, to assess ICU mortality risk (Le Gall et al. 1984, 1993; Moreno et al. 2005);

• APACHE I, II and III, to assess ICU mortality risk (Knaus et al. 1981, 1985, 1991);

• CHADS$_2$, to assess the risk of stroke in patients with atrial fibrillation (Gage et al. 2001);

• Wells Criteria for pulmonary embolisms (Wells et al. 2000); Wells Criteria for and deep vein thrombosis (Wells et al. 1997);

• TIMI, to assess the risk of death and ischemic events (Antman et al. 2000);

• SIRS, to detect system inflammatory response syndrome (Bone et al. 1992);

All of the scoring systems listed above are sparse linear models with small coprime coefficients. The CHADS$_2$ scoring system, for instance, uses 5 coefficients with values of 1 and 2.

Many medical scoring systems were constructed without optimizing for predictive accuracy. In some cases, physicians built scoring systems by combining existing linear classification methods with heuristics. The SAPS II score, for instance, was constructed by rounding logistic regression coefficients: as Le Gall et al. (1993) write, “the general rule was to multiply the $\beta$ for each range by 10 and round off to the nearest integer.” This approach is at odds with the fact that rounding is known to produce suboptimal solutions in the field of integer programming. In other cases, scoring systems were hand-crafted by a panel of physicians, and not learned from data at all. This appears to have been the case for CHADS$_2$ as suggested by Gage et al. (2001): “To create CHADS$_2$, we assigned 2 points to a history of prior cerebral ischemia and 1 point for the presence of other risk factors because a history of prior cerebral ischemia increases the relative risk (RR) of subsequent stroke commensurate to 2 other risk factors combined. We calculated CHADS$_2$, by adding 1 point each for each of the following—recent CHF, hypertension, age 75 years or older, and DM—and 2 points for a history of stroke or TIA.” Methods that can easily produce highly tailored prediction models, such as SLIM, should eliminate the need for physicians to create models by hand.

In addition to the sleep apnea application that we present in Sect. 6, SLIM has also been used to create medical scoring system for diagnosing cognitive impairments using features derived from a clock-drawing test (Souillard-Mandar et al. 2015), and to create scoring systems for recidivism prediction (Zeng et al. 2015).

1.1.2 Sparse linear classification models

In comparison to SLIM, the majority of current methods for sparse linear classification are designed to fit models with real coefficients and would therefore require rounding to yield scoring systems. In practice, rounding the coefficients of a linear model may significantly alter its accuracy and sparsity, and may produce a scoring system that violates operational constraints on these quantities. Further, many current methods also control accuracy and sparsity by means of convex surrogate functions to preserve scalability (see e.g., Tibshirani 1996; Efron et al. 2004). As we show in Sects. 6 and 7, surrogate functions provide a poor trade-off between accuracy and sparsity. Convex surrogate loss functions, for instance, produce models that do not attain the best learning-theoretic guarantee on predictive accuracy and are not robust to outliers (Li and Lin 2007; Brooks and Lee 2010; Nguyen and Sanner 2013). Similarly, $\ell _1$-regularization is only guaranteed to recover the correct sparse solution (i.e., the one that minimizes the $\ell _0$-norm) under restrictive conditions that are rarely satisfied in practice (Zhao and Bin 2007; Liu and Zhang 2009). In fact, $\ell _1$-regularization may recover a solution that attains a significant loss in predictive accuracy relative to the correct sparse solution (see e.g., Lin et al. 2008, for a discussion). Sparse linear classifiers can also be produced using feature selection algorithms (Guyon and Elisseeff 2003; Mao 2004), though these algorithms cannot guarantee an optimal balance between accuracy and sparsity as they typically rely on greedy optimization (with some exceptions, see e.g., Bradley et al. 1999).

1.1.3 Discrete linear classification models

SLIM is part of a recent stream of research on creating linear classifiers with discrete coefficients. Specifically, Chevaleyre et al. (2013) have considered training linear classifiers with binary coefficients by rounding the coefficients of linear classifiers that minimize the hinge loss. In addition, Carrizosa et al. (2013) have considered training linear classifiers with small integer coefficients by using a MIP formulation that optimizes the hinge loss. The discretization bounds and generalization bounds in Sect. 4 are a novel contribution to this body of work and applicable to all linear models with discrete coefficients.

SLIM can reproduce the models proposed by Chevaleyre et al. (2013) and Carrizosa et al. (2013) (see e.g., our formulation to create M-of-N rule tables in Sect. 8.2.1). The converse, however, is not necessarily true because the methods of Chevaleyre et al. (2013) and Carrizosa et al. (2013) have the following weaknesses: (i) they optimize the hinge loss as opposed to the 0–1 loss; and (ii) they do not include a mechanism to control sparsity. These differences may result in better scalability compared to SLIM. However, they also eliminate the ability of these methods to produce scoring systems that are sparse, that satisfy operational constraints on accuracy and/or sparsity, and that can be trained without parameter tuning.

1.1.4 MIP approaches for classification

SLIM uses integer programming (IP) to achieve three distinct goals: (i) minimize the 0–1 loss; (ii) penalize the ${\ell _{0}}$-norm for feature selection; and (iii) restrict coefficients to a small set of integers. MIP approaches have been used to tackle each of these goals, albeit separately.

MIP formulations to minimize the 0–1 loss, for instance, were first proposed in Liittschwager and Wang (1978) and Bajgier and Hill (1982), and later refined by Mangasarian (1994), Asparoukhov and Stam (1997) and Glen (1999). Similarly, MIP formulations that penalize the $\ell _0$-norm for feature selection were proposed in Goldberg and Eckstein (2010), Goldberg and Eckstein (2012) and Nguyen and Franke (2012). To our knowledge, the only MIP formulation to restrict coefficients to a small set of integers is proposed in Carrizosa et al. (2013).

SLIM has unique practical benefits in comparison to these MIP formulations since it handles all three of these goals simultaneously. As we discuss in Sect. 2, the simultaneous approach allows SLIM to train models parameter tuning, and make use of the variables that encode the 0–1 loss and ${\ell _{0}}$-norm to accommodate important operational constraints. Further, restricting coefficients to a discrete set that is finite leads to an IP formulation whose LP relaxation is significantly tighter than other MIP formulations designed to minimize the 0–1 loss and/or penalize the $\ell _0$-norm.

The problem of finding a linear classifier that minimizes the 0–1 loss function is sometimes referred to as misclassification minimization for the linear discriminant problem in the MIP community (see e.g., Rubin 2009; Lee and Wu 2009, for an overview). Seeing how early attempts at misclassification minimization were only feasible for tiny datasets with at most 200 examples (see e.g., Joachimsthaler and Stam 1990), a large body of work has focused on improving scalability by using heuristics (Rubin 1990; Yanev and Balev 1999), decomposition procedures (Rubin 1997), cutting planes (Brooks 2011) and specialized branch-and-bound algorithms (Nguyen and Sanner 2013). The data reduction technique we present in Sect. 5 is a novel contribution to this body of work, and addresses a need for general methods to remove redundant data put forth by Bradley et al. (1999). We compare and contrast data reduction to existing approaches in greater detail in Sect. 5.3.

2 Methodology

We start with a dataset of N i.i.d. training examples $\mathcal {D}_N = \{(\varvec{x}_i,y_i)\}_{i=1}^N$ where each $\varvec{x}_i \in \mathcal {X}\subseteq \mathbb {R}^{P+1}$ denotes a vector of features $[1, x_{i,1},\ldots ,x_{i,P}]^T$ and each $y_i \in \mathcal {Y}= \{-1,1\}$ denotes a class label. We consider linear classification models of the form $\hat{y}=\text {sign}(\varvec{\lambda }^T\varvec{x})$, where $\varvec{\lambda }\subseteq \mathbb {R}^{P+1}$ represents a vector of coefficients $[\lambda _0, \lambda _1,\ldots ,\lambda _P]^T$ and $\lambda _0$ represents an intercept term. We learn the values of the coefficients from the data $\mathcal {D}_N$ by solving an optimization problem of the form:

$$\begin{aligned} \begin{aligned} \min _{\varvec{\lambda }}&\qquad {\text {Loss}}\left( \varvec{\lambda };\mathcal {D}_N\right) + C \cdot \varPhi (\varvec{\lambda }) \\ {\text {s.t.}}&\qquad \varvec{\lambda }\in \mathcal {L}. \end{aligned} \end{aligned}$$

(1)

Here: the loss function, ${\text {Loss}}\left( \varvec{\lambda };\mathcal {D}_N\right) : \mathbb {R}^{P+1} \times (\mathcal {X}\times \mathcal {Y})^N \rightarrow \mathbb {R}$, penalizes misclassifications; the interpretability penalty, $\varPhi (\varvec{\lambda }): \mathbb {R}^{P+1} \rightarrow \mathbb {R}$, induces soft qualities that are desirable but may be sacrificed for greater accuracy; the interpretability set, $\mathcal {L}$, encodes hard qualities that are absolutely required; and the trade-off parameter, C, controls the balance between accuracy and soft qualities. We assume: (i) the interpretability set contains the null vector so that $\mathbf {0} \in \mathcal {L}$; (ii) the interpretability penalty is additively separable so that $\varPhi (\varvec{\lambda }) = \sum _{j=0}^P\varPhi _j(\lambda _j)$; (iii) the intercept is never penalized so that $\varPhi _0(\lambda _0) = 0$.

A Supersparse Linear Integer Model (SLIM) is a special case of the optimization problem in (1):

$$\begin{aligned} \begin{aligned} \min _{\varvec{\lambda }}&\qquad \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0 + \epsilon \left\| \varvec{\lambda }\right\| _1 \\ {\text {s.t.}}&\qquad \varvec{\lambda }\in \mathcal {L}. \end{aligned} \end{aligned}$$

(2)

SLIM directly optimizes accuracy and sparsity by minimizing the 0–1 loss $\frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right]$ and ${\ell _{0}}$-norm $\left\| \varvec{\lambda }\right\| _0:=\sum _{j=1}^P{\mathbbm {1}\left[ \lambda _j\ne 0\right] }$ respectively. The constraints usually restrict coefficients to a finite set of discrete values such as $\mathcal {L}= \{-10,\ldots ,10\}^{P+1}$, and may include additional operational constraints such as $\left\| \varvec{\lambda }\right\| _0 \le 10$. SLIM includes a tiny ${\ell _{1}}$-penalty $\epsilon \left\| \varvec{\lambda }\right\| _1$ in the objective for the sole purpose of restricting coefficients to coprime values. ¹ To be clear, the ${\ell _{1}}$-penalty parameter $\epsilon$ is always set to a value that is small enough to avoid ${\ell _{1}}$-regularization (that is, $\epsilon$ is small enough to guarantee that SLIM never sacrifices accuracy or sparsity to attain a smaller ${\ell _{1}}$-penalty).

SLIM is designed to produce scoring systems that attain a pareto-optimal trade-off between accuracy and sparsity: when we minimize 0–1 loss and the $\ell _0$-penalty, we only sacrifice classification accuracy to attain higher sparsity, and vice versa. Minimizing the 0–1 loss produces scoring systems that are completely robust to outliers and attain the best learning-theoretic guarantee on predictive accuracy (see e.g., Brooks 2011; Nguyen and Sanner 2013). Similarly, controlling for sparsity via ${\ell _{0}}$-regularization prevents the additional loss in accuracy due to $\ell _1$-regularization (see Lin et al. 2008, for a discussion). We can make further use of the variables that encode the 0–1 loss and ${\ell _{0}}$-penalty to formulate operational constraints related to accuracy and sparsity (see Sect. 3).

One unique benefit in minimizing an approximation-free objective function over a finite set of discrete coefficients is that the free parameters in the objective of (2) have special properties.

Remark 1

If $\epsilon <\frac{\min {(1/N,C_0})}{\max _{\varvec{\lambda }\in \mathcal {L}}\left\| \varvec{\lambda }\right\| _1}$ and $\mathcal {L}$ is a finite subset of $\mathbb {Z}^{P+1}$ then the optimization of (2) will produce a scoring system with coprime integer coefficients without affecting accuracy or sparsity. That is,

$$\begin{aligned}&\mathop {{{\mathrm{argmin}}}}\limits _{\varvec{\lambda }\in \mathcal {L}} \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0 + \epsilon \left\| \varvec{\lambda }\right\| _1 \subseteq \\&\quad \mathop {{{\mathrm{argmin}}}}\limits _{\varvec{\lambda }\in \mathcal {L}}\frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0 \end{aligned}$$

and,

$$\begin{aligned} \text {gcd}(\{\lambda _j^*\}_{j=0}^P)=1 \text { for all } \varvec{\lambda }^* \in \mathop {{{\mathrm{argmin}}}}\limits _{\varvec{\lambda }\in \mathcal {L}} \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0 + \epsilon \left\| \varvec{\lambda }\right\| _1. \end{aligned}$$

Remark 2

The trade-off parameter $C_0$ represents the maximum accuracy that SLIM will sacrifice to remove a feature from the optimal scoring system.

Remark 3

If $C_0 < \frac{1}{NP}$ and $\epsilon < \frac{\min {(1/N,C_0})}{\max _{\varvec{\lambda }\in \mathcal {L}}\left\| \varvec{\lambda }\right\| _1} = \frac{C_0}{\max _{\varvec{\lambda }\in \mathcal {L}}\left\| \varvec{\lambda }\right\| _1}$ then the optimization of (2) will produce a scoring system with coefficients $\varvec{\lambda }\in \mathcal {L}$ that attains the highest possible training accuracy. That is,

$$\begin{aligned} \mathop {{{\mathrm{argmin}}}}\limits _{\varvec{\lambda }\in \mathcal {L}} \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0 + \epsilon \left\| \varvec{\lambda }\right\| _1 \subseteq \mathop {{{\mathrm{argmin}}}}\limits _{\varvec{\lambda }\in \mathcal {L}} \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] . \end{aligned}$$

Remark 4

If $C_0 > 1 - \frac{1}{N}$ and $\epsilon <\frac{\min {(1/N,C_0)}}{\max _{\varvec{\lambda }\in \mathcal {L}}\left\| \varvec{\lambda }\right\| _1} = \frac{1/N}{\max _{\varvec{\lambda }\in \mathcal {L}}\left\| \varvec{\lambda }\right\| _1}$ then the optimization of (2) will produce a scoring system with coefficients $\varvec{\lambda }\in \mathcal {L}$ that attains the highest possible sparsity. That is,

$$\begin{aligned} \mathop {{{\mathrm{argmin}}}}\limits _{\varvec{\lambda }\in \mathcal {L}} \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0 + \epsilon \left\| \varvec{\lambda }\right\| _1 \subseteq \mathop {{{\mathrm{argmin}}}}\limits _{\varvec{\lambda }\in \mathcal {L}} C_0 \left\| \varvec{\lambda }\right\| _0. \end{aligned}$$

The aforementioned properties are only possible using the formulation in (2). In particular, the properties in Remarks 2–4 require that we control accuracy using the 0–1 loss and control sparsity using an $\ell _0$-penalty. In addition, the property in Remark 1 requires that we restrict coefficients to a finite set of discrete values.

2.1 SLIM IP Formulation

We train SLIM scoring systems using the following IP formulation:

$$\begin{aligned} \min _{\varvec{\lambda },{\varvec{\psi }},{\varvec{\varPhi }},{\varvec{\alpha }},{\varvec{\beta }}}&\qquad \frac{1}{N}\sum _{i=1}^{N} \psi _i + \sum _{j=1}^{P} \varPhi _j \nonumber \\ {\text {s.t.}}&\qquad M_i \psi _i \ge \gamma -\sum _{j=0}^P y_i \lambda _j x_{i,j}&{i}={1}\text {,...,}{N}&\qquad \textit{ 0-1 loss} \end{aligned}$$

(3a)

$$\begin{aligned}&\qquad \varPhi _j = C_0\alpha _j + \epsilon \beta _j&{j}={1}\text {,...,}{P}&\qquad \textit{ int. penalty} \end{aligned}$$

(3b)

$$\begin{aligned}&\qquad -\varLambda _j\alpha _j \le \lambda _j \le \varLambda _j\alpha _j&{j}={1}\text {,...,}{P}&\qquad {\ell _{0}}\textit{ -norm} \end{aligned}$$

(3c)

$$\begin{aligned}&\qquad -\beta _j \le \lambda _j \le \beta _j&{j}={1}\text {,...,}{P}&\qquad {\ell _{1}}\textit{ -norm} \nonumber \\&\qquad \lambda _j \in \mathcal {L}_j&{j}={0}\text {,...,}{P}&\qquad \textit{ coefficient set} \nonumber \\&\qquad \psi _i \in \{0,1\}&{i}={1}\text {,...,}{N}&\qquad \textit{ loss variables} \nonumber \\&\qquad \varPhi _j \in \mathbb {R}_+&{j}={1}\text {,...,}{P}&\qquad \textit{ penalty variables} \nonumber \\&\qquad \alpha _j \in \{0,1\}&{j}={1}\text {,...,}{P}&\qquad {\ell _{0}}\textit{ variables} \nonumber \\&\qquad \beta _j \in \mathbb {R}_+&{j}={1}\text {,...,}{P}&\qquad {\ell _{1}}\textit{ variables} \end{aligned}$$

(3d)

Here, the constraints in (3a) set the loss variables $\psi _i = \mathbbm {1}\left[ y_i \varvec{\lambda }^T\varvec{x}_i \le 0\right]$ to 1 if a linear classifier with coefficients $\varvec{\lambda }$ misclassifies example i. This is a Big-M constraint for the 0–1 loss that depends on scalar parameters $\gamma$ and $M_i$ (see e.g., Rubin 2009). The value of $M_i$ represents the maximum score when example i is misclassified, and can be set as $M_i = \max _{\varvec{\lambda }\in \mathcal {L}} (\gamma - y_i\varvec{\lambda }^T\varvec{x}_i)$ which is easy to compute since $\varvec{\lambda }$ is restricted to a finite discrete set. The value of $\gamma$ represents the “margin” and should be set as a lower bound on $y_i\varvec{\lambda }^T\varvec{x}_i$. When the features are binary, $\gamma$ can be set to any value between 0 and 1. In other cases, the lower bound is difficult to calculate exactly so we set $\gamma =0.1$, which makes an implicit assumption on the values of the features. The constraints in (3b) set the total penalty for each coefficient to $\varPhi _j = C_0 \alpha _j + \epsilon \beta _j$, where $\alpha _j := \mathbbm {1}\left[ \lambda _j\ne 0\right]$ is defined by Big-M constraints in (3c), and $\beta _j := |\lambda _j|$ is defined by the constraints in (3d). We denote the largest absolute value of each coefficient as $\varLambda _j := \max _{\lambda _j\in \mathcal {L}_j} |\lambda _j|$.

Restricting coefficients to a finite discrete set results in significant practical benefits for the SLIM IP formulation, especially in comparison to other IP formulations that minimize the 0–1-loss and/or penalize the $\ell _0$-norm. Many IP formulations compute the 0–1 loss and $\ell _0$-norm by means of Big-M constraints that use on Big-M constants (see e.g., Wolsey 1998). Restricting the coefficients to a finite discrete set allows us to bound Big-M constants that are typically set heuristically to a large value. Specifically, the Big-M constant for computing the 0–1 loss in constraints (3a) is bounded as $M_i \le \max _{\varvec{\lambda }\in \mathcal {L}} (\gamma - y_i\varvec{\lambda }^T\varvec{x}_i)$ (compare with Brooks 2011, where the same parameter has to be approximated by a “sufficiently large constant”). Similarly, the Big-M constant used to compute the $\ell _0$-norm in constraints (3c) is bounded at $\varLambda _j \le \max _{\lambda _j\in \mathcal {L}_j} |\lambda _j|$ (compare with Guan et al. 2009, where this same parameter has to be approximated by a “sufficiently large constant”). These differences lead to a tighter LP relaxation, which narrows the integrality gap, and subsequently improves the ability of commercial IP solvers to obtain a proof of optimality.

3 Operational constraints

In this section, we discuss how SLIM can accommodate a wide range of operational constraints related to the accuracy and sparsity of predictive models. The following techniques provide users with a practical approach to create tailored prediction models. They are made possible by the facts that: (i) the variables that encode the 0–1 loss and ${\ell _{0}}$-penalty in the SLIM IP formulation can also be used to handle accuracy and sparsity; and (ii) the free parameters in the SLIM objective can be set without tuning (see Remarks 2–4).

3.1 Loss constraints for imbalanced data

The majority of classification problems in the medical domain are imbalanced. Handling imbalanced data is incredibly difficult for most classification methods since maximizing classification accuracy often produces a trivial model (i.e., if the probability of heart attack is 1 %, a model that never predicts a heart attack is still 99 % accurate).

SLIM has unique benefits when training scoring systems on imbalanced problems. Specifically, it is the only method that can train a supervised classification model at any point on the ROC curve without any parameter tuning. When physicians specify hard constraints on sensitivity (or specificity), we can encode these as loss constraints into the IP formulation, and solve the IP to obtain the least specific (or most sensitive) model without parameter tuning. To train the most sensitive scoring system with a maximum error of $\gamma \in [0,1]$ on negatively-labeled examples we solve an IP with the form:

$$\begin{aligned} \min _{\varvec{\lambda }}&\qquad {\frac{W}{N}}^+ \sum _{i\in {\mathcal {I}}^{+}} \mathbbm {1}\left[ y_i \varvec{\lambda }^T\varvec{x}_i \le 0\right] + \frac{W}{N}^{-} \sum _{i\in {\mathcal {I}}^{-}} \mathbbm {1}\left[ y_i \varvec{\lambda }^T\varvec{x}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0 + \epsilon \left\| \varvec{\lambda }\right\| _1 \nonumber \\ {\text {s.t.}}&\qquad \frac{1}{N^-} \sum _{i\in {\mathcal {I}}^{-}} \mathbbm {1}\left[ y_i \varvec{\lambda }^T\varvec{x}_i \ge 0\right] \le \gamma \nonumber \\&\qquad \varvec{\lambda }\in \mathcal {L}. \end{aligned}$$

(4)

Here, we have used a weighted 0–1 loss function where $W^{+}$ and $W^{-}$ are weights that control the accuracy on the $N^{+}$ positively-labeled examples from the set ${\mathcal {I}}^{+}= \{i:y_i = +1\}$, and $N^{-}$ negatively-labeled examples from the set ${\mathcal {I}}^{-}= \{i:y_i = -1\}$, respectively. Assuming that $W^{+}+ W^{-}= 1$, we can set $W^{+}> \frac{N^{-}}{1+N^{-}}$ and $W^{-}= 1 - W^{+}$ so that the optimization aims to find a scoring system that classifies all of the positively-labeled examples correctly, at the expense of misclassifying all of the negatively-labeled examples. Constraint (4) prevents this from happening by limiting the error on negatively-labeled examples to $\gamma$. Thus, the optimal scoring system attains the highest sensitivity among models with a maximum error of $\gamma$ on negatively-labeled examples.

3.2 Feature-based constraints for input variables

SLIM provides fine-grained control over the composition of input variables in a scoring system by formulating feature-based constraints. Specifically, we can make use of the indicator variables that encode the $\ell _0$-norm $\alpha _j := \mathbbm {1}\left[ \lambda _j \ne 0\right]$ to formulate arbitrarily complicated logical constraints between features such as “either-or” conditions or “if-then” conditions (see e.g., Wolsey 1998). This presents a practical alternative to create classification models that obey structured sparsity constraints (see e.g., Jenatton et al. 2011) or hierarchical constraints (see e.g., Bien et al. 2013).

The indicator variables $\alpha _j$ can be used to limit the number of input variables to at most $\varTheta$ by adding the constraint,

$$\begin{aligned} \sum _{j = 1}^P{\alpha _j} \le \varTheta . \end{aligned}$$

More complicated constraints include, for example, “if-then” constraints to ensure that a model will include hypertension and $heart\_attack$ if it also includes stroke:

$$\begin{aligned} \alpha _{heart\_attack} + \alpha _{hypertension} \le 2 \alpha _{stroke}, \end{aligned}$$

or hierarchical constraints to ensure that an input variable in the leaves can only be used when all features above it in the hierarchy are also used:

$$\begin{aligned} \alpha _{leaf} \le \alpha _{node} {\text { for all nodes above the leaf}}. \end{aligned}$$

3.3 Feature-based preferences

Physicians often have soft preferences between different input variables. SLIM allows practitioners to encode these preferences by specifying a distinct trade-off parameter for each coefficient $C_{0,j}$.

Explicitly, when our model should use feature j instead of feature k, we set $C_{0,k} = C_{0,j} + \delta$, where $\delta > 0$ represents the maximum additional training accuracy that we are willing to sacrifice in order to use feature j instead of feature k. Thus, setting $C_{0,k} = C_{0,j} + 0.02$ would ensure that we would only be willing to use feature k instead of feature j if it yields an additional 2 % gain in training accuracy over feature k.

This approach can also be used to handle problems with missing data. Consider training a model where feature j contains $M<N$ missing points. Instead of dropping these points, we can impute the values of the M missing examples, and adjust the trade-off parameter $C_{0,j}$ so that our model only uses feature j if it yields an additional gain in accuracy of more than M examples:

$$\begin{aligned} C_{0,j} = C_0 + \frac{M}{N}. \end{aligned}$$

The adjustment factor is chosen so that: if $M=0$ then $C_{0,j} = C_0$ and if $M=N$ then $C_{0,j} = 1$ and the coefficient is dropped entirely (see Remark 4). This ensures that features with lots of imputed values are more heavily penalized than features with fewer imputed values.

4 Bounds on training and testing accuracy

In this section, we present bounds on the training and testing accuracy of SLIM scoring systems.

4.1 Discretization bounds on training accuracy

Our first result shows that we can always craft a finite discrete set of coefficients $\mathcal {L}$ so that the training accuracy of a linear classifier with discrete coefficients $\varvec{\lambda }\in \mathcal {L}$ (e.g. SLIM) is no worse than the training accuracy of a baseline linear classifier with real-valued coefficients $\varvec{\rho }\in \mathbb {R}^P$ (e.g. SVM).

Theorem 1

(Minimum margin resolution bound) Let $\varvec{\rho }= [\rho _1,\ldots ,\rho _P]^T \in \mathbb {R}^P$ denote the coefficients of a baseline linear classifier trained using data $\mathcal {D}_N = (\varvec{x}_i,y_i)_{i=1}^N$. Let $X_{\max } = \max _i \Vert \varvec{x}_i\Vert _2$ and $\gamma _{\min } = \min _i \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2}$ denote the largest magnitude and minimum margin achieved by any training example, respectively.

Consider training a linear classifier with coefficients $\varvec{\lambda }= [\lambda _1,\ldots ,\lambda _P]^T$ from the set $\mathcal {L}= \{-\varLambda ,\ldots ,\varLambda \}^P$. If we choose a resolution parameter $\varLambda$ such that:

$$\begin{aligned} \varLambda&> \frac{X_{\max }\sqrt{P}}{2 \gamma _{\min }}, \end{aligned}$$

(5)

then there exists $\varvec{\lambda }\in \mathcal {L}$ such that the 0–1 loss of $\varvec{\lambda }$ is less than or equal to the 0–1 loss of $\varvec{\rho }$:

$$\begin{aligned} \sum _{i=1}^N \mathbbm {1}\left[ y_i\varvec{\lambda }^T\varvec{x}_i\le 0\right] \le \sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\rho }^T \varvec{x}_i \le 0\right] . \end{aligned}$$

Proof

See Appendix. $\square$

The proof of Theorem 1 uses a rounding procedure to choose a resolution parameter $\varLambda$ so that the coefficient set $\mathcal {L}$ contains a classifier with discrete coefficients $\varvec{\lambda }$ that attains the same the 0–1 loss as the baseline classifier with real coefficients $\varvec{\rho }$. If the baseline classifier $\varvec{\rho }$ is obtained by minimizing a convex surrogate loss, then the optimal SLIM classifier trained with the coefficient set from Theorem 1 may attain a lower 0–1 loss than $\mathbbm {1}\left[ y_i\varvec{\rho }^T\varvec{x}_i\le 0\right]$ because SLIM directly minimizes the 0–1 loss.

The next corollary yields additional bounds on the training accuracy by considering progressively larger values of the margin. These bounds can be used to relate the resolution parameter $\varLambda$ to a worst-case guarantee on training accuracy.

Corollary 1

(kth margin resolution bound) Let $\varvec{\rho }= [\rho _1,\ldots ,\rho _P]^T \in \mathbb {R}^P$ denote the coefficients of a linear classifier trained with data $\mathcal {D}_N = (\varvec{x}_i,y_i)_{i=1}^N$. Let $\gamma _{(k)}$ denote the value of the kth smallest margin, $\mathcal {I}_{(k)}$ denote the set of training examples with $\frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2} \le \gamma _{(k)}$, and $X_{(k)} = \max _{i \not \in \mathcal {I}_{(k)}} \Vert \varvec{x}_i\Vert _2$ denote the largest magnitude of any training example $\varvec{x}_i \in \mathcal {D}_N$ for $i \not \in \mathcal {I}_{(k)}$.

Consider training a linear classifier with coefficients $\varvec{\lambda }= [\lambda _1,\ldots ,\lambda _P]^T$ from the set $\mathcal {L}= \{-\varLambda ,\ldots ,\varLambda \}^P$. If we choose a resolution parameter $\varLambda$ such that:

$$\begin{aligned} \varLambda > \frac{X_{(k)}\sqrt{P}}{2 \gamma _{(k)}}, \end{aligned}$$

then there exists $\varvec{\lambda }\in \mathcal {L}$ such that the 0–1 loss of $\varvec{\lambda }$ and the 0–1 loss of $\varvec{\rho }$ differ by at most $k-1$:

$$\begin{aligned} \sum _{i=1}^N \mathbbm {1}\left[ y_i\varvec{\lambda }^T\varvec{x}_i\le 0\right] - \sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\rho }^T \varvec{x}_i \le 0\right] \le k - 1. \end{aligned}$$

Proof

The proof follows by applying Theorem 1 to the reduced dataset $\mathcal {D}_N \backslash \mathcal {I}_{(k)}$. $\square$

We have now shown that good discretized solutions exist and can be constructed easily. This motivates that optimal discretized solutions, which by definition are better than rounded solutions, will also be good relative to the best non-discretized solution.

4.2 Generalization bounds on testing accuracy

According to the principle of structural risk minimization (Vapnik 1998), fitting a classifier from a simpler class of models may lead to an improved guarantee on predictive accuracy. Consider training a classifier $f:\mathcal {X}\rightarrow \mathcal {Y}$ with data $\mathcal {D}_N = (\varvec{x}_i,y_i)_{i=1}^N$, where $\varvec{x}_i \in \mathcal {X}\subseteq \mathbb {R}^P$ and $y_i \in \mathcal {Y}= \{-1,1\}$. In what follows, we provide uniform generalization guarantees on the predictive accuracy of all functions, $f \in \mathcal {F}$. These guarantees bound the true risk $R^\text {true}(f) = \mathbb {E}_{\mathcal {X},\mathcal {Y}} \mathbbm {1}\left[ f(\varvec{x}) \ne y\right]$ by the empirical risk $R^\text {emp}(f) = \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ f(\varvec{x}_i) \ne y_i\right]$ and other quantities important to the learning process.

Theorem 2

(Occam’s Razor) Let $\mathcal {F}$ denote the set of linear classifiers with coefficients $\varvec{\lambda }\in \mathcal {L}$:

$$\begin{aligned} \mathcal {F}= \left\{ f :\mathcal {X}\rightarrow \mathcal {Y}\;\big |\; f(\varvec{x})={\text {sign}}\left( \varvec{\lambda }^T\varvec{x}\right) \quad \hbox { and }\quad \varvec{\lambda }\in \mathcal {L}\right\} . \end{aligned}$$

For every $\delta > 0,$ with probability at least $1-\delta$, every classifier $f\in \mathcal {F}$ obeys:

$$\begin{aligned} R^\mathrm{true}(f) \le R^\mathrm{emp}(f) + \sqrt{\frac{\log (|\mathcal {L}|) - \log (\delta )}{2N}}. \end{aligned}$$

A proof of Theorem 2 can be found in Section 3.4 of Bousquet et al. (2004). The result that more restrictive hypothesis spaces can lead to better generalization provides motivation for using discrete models without necessarily expecting a loss in predictive accuracy. The bound indicates that we include more coefficients in the set $\mathcal {L}$ as the amount of data N increases.

In Theorem 3, we improve the generalization bound from Theorem 2 by exploiting the fact that we can bound the number of non-zero coefficients in a SLIM scoring system in terms of the value of $C_0$.

Theorem 3

(Generalization of Sparse Discrete Linear Classifiers) Let $\mathcal {F}$ denote the set of linear classifiers with coefficients $\varvec{\lambda }$ from a finite set $\mathcal {L}$ such that:

$$\begin{aligned} \mathcal {F}&= \left\{ f :\mathcal {X}\rightarrow \mathcal {Y}\;\big |\; f(\varvec{x})={\text {sign}}\left( \varvec{\lambda }^T\varvec{x}\right) \right\} \\ \varvec{\lambda }&\in \mathop {{{\mathrm{argmin}}}}\limits _{\lambda \in \mathcal {L}} \frac{1}{N} \sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0\left\| \varvec{\lambda }\right\| _0 \end{aligned}$$

For every $\delta > 0,$ with probability at least $1-\delta$, every classifier $f\in \mathcal {F}$ obeys:

$$\begin{aligned} R^{\mathrm{true}}(f)&\le R^{\mathrm{emp}}(f) + \sqrt{\frac{\log (|\mathcal {H}_{P,C_0}|) - \log (\delta )}{2N}}. \end{aligned}$$

where

$$\begin{aligned} \mathcal {H}_{P,C_0}&= \bigg \{\varvec{\lambda }\in \mathcal {L}\;\Big |\; \left\| \varvec{\lambda }\right\| _0 \le \left\lfloor \frac{1}{C_0} \right\rfloor \bigg \}. \end{aligned}$$

Proof

See Appendix. $\square$

This theorem relates the trade-off parameter $C_0$ in the SLIM objective to the generalization of SLIM scoring systems. It indicates that increasing the value of the $C_0$ parameter will produce a model with better generalization properties.

In Theorem 4, we produce a better generalization bound by exploiting the fact that SLIM scoring systems use coprime integer coefficients (see Remark 1). In particular, we express the generalization bound from Theorem 2 using the P-dimensional Farey points of level $\varLambda$ (see Marklof 2012, for a definition).

Theorem 4

(Generalization of discrete linear classifiers with coprime coefficients) Let $\mathcal {F}$ denote the set of linear classifiers with coprime integer coefficients, $\varvec{\lambda }$, bounded by $\varLambda$:

$$\begin{aligned} \mathcal {F}&= \Big \{f :\mathcal {X}\rightarrow \mathcal {Y}\;\big |\; f(\varvec{x})={\text {sign}}\left( \varvec{\lambda }^T\varvec{x}\right) \hbox { and } \varvec{\lambda }\in \mathcal {L}\Big \},\\ \mathcal {L}&= \Big \{\varvec{\lambda }\in \mathbb {\hat{Z}}^P \;\big |\; |\lambda _j| \le \varLambda \quad \hbox { for }\quad j=1,\ldots ,P\Big \},\\ \hat{\mathbb {Z}}^P&= \Big \{\varvec{z}\in \mathbb {Z}^P \;\big |\; \hbox {gcd}(\varvec{z}) = 1\Big \}. \end{aligned}$$

For every $\delta > 0,$ with probability at least $1-\delta$, every classifier $f\in \mathcal {F}$ obeys:

$$\begin{aligned} R^{\mathrm{true}}(f)&\le R^{\mathrm{emp}}(f) + \sqrt{\frac{\log (|\mathcal {C}_{P,\varLambda }|) - \log (\delta )}{2N}}, \end{aligned}$$

where $\mathcal {C}_{P,\varLambda }$ denotes the set of Farey points of level $\varLambda$:

$$\begin{aligned} \mathcal {C}_{P,\varLambda }&= \left\{ \frac{\varvec{\lambda }}{q} \in [0,1)^P: (\varvec{\lambda },q) \in \mathbb {\hat{Z}}^{P+1}\quad \text { and }\quad 1 \le q \le \varLambda \right\} . \end{aligned}$$

The proof involves a counting argument over coprime integer vectors, using the definition of Farey points from number theory.

In Fig. 1, we plot the relative density of coprime integer vectors in $\mathbb {Z}^P$ bounded by $\varLambda$ (i.e., $|\mathcal {C}_{P,\varLambda }|/(2\varLambda +1)^P$) and the relative improvement in the generalization bound due to the use of coprime coefficients. This shows that using coprime coefficients can significantly reduce the number of classifiers based on the dimensionality of the data and the value of $\varLambda$. The corresponding improvement in the generalization bound may be significant when the data are high dimensional and $\varLambda$ is small.

Fig. 1

Relative density of coprime integer vectors in $\mathbb {Z}^P$ (left), and the relative improvement in the generalization bound due to the use of coprime coefficients for $\delta =0.01$ (right)

5 Data reduction

Data reduction is a technique that can decrease the computation associated with training a supervised classification model by discarding redundant training data. This technique can be applied to any supervised classification method where the training procedure is carried out by solving an optimization problem. However, it is best suited for methods such as SLIM, where the underlying optimization problem may be difficult to solve for large instances. In this section, we first describe how data reduction works in a general setting, and then show how it can be applied to SLIM.

5.1 Data reduction for optimization-based supervised classification

Consider training a classifier $f: \mathcal {X}\rightarrow \mathcal {Y}$ by solving a computationally challenging optimization problem,

$$\begin{aligned} \min _{f} ~ Z(f;\mathcal {D}_N) ~\text {s.t.}~ f\in \mathcal {F}. \end{aligned}$$

(6)

We refer to the optimization problem in (6) as the original problem. Here, $\mathcal {F}$ represents the set of feasible classifiers and $Z:\mathcal {F}\times (\mathcal {X}\times \mathcal {Y})^N \rightarrow \mathbb {R}$ represents its objective function.

Data reduction aims to decrease the computation associated with solving the original problem by removing redundant examples from $\mathcal {D}_N = (\varvec{x}_i,y_i)_{i=1}^N$ (i.e., data points that can be safely discarded without changing the optimal solution to (6)). The technique requires users to specify a surrogate problem that is considerably easier to solve. Given the initial training data $\mathcal {D}_N = (\varvec{x}_i,y_i)_{i=1}^N$, and the surrogate problem, data reduction solves $N+1$ variants of the surrogate problem to identify redundant examples. These examples are then removed from the initial training data to leave behind a subset of reduced training data $\mathcal {D}_M \subseteq \mathcal {D}_N$ that is guaranteed to yield the same optimal classifier as $\mathcal {D}_N$. Thus, the computational gain from data reduction comes from training a model with $\mathcal {D}_M$ (i.e., solving an instance of the original problem with $N-M$ fewer examples).

We provide an overview of data reduction in Algorithm 1. To explain how the algorithm works, let us denote the surrogate problem as:

$$\begin{aligned} \min _{f} ~ \tilde{Z}(f;\mathcal {D}_N) ~\text {s.t.}~ f\in \tilde{\mathcal {F}}. \end{aligned}$$

(7)

Here $\tilde{Z}: \tilde{\mathcal {F}} \times (\mathcal {X}\times \mathcal {Y})^N \rightarrow \mathbb {R}$ denotes the objective function of the surrogate problem, and $\tilde{\mathcal {F}}$ denotes its set of feasible classifiers. Data reduction can be used with any surrogate problem so long as the $\varepsilon$-level set of the surrogate problem contains all optimizers to the original problem. That is, we can use any feasible set $\tilde{\mathcal {F}}$ and any objective function $\tilde{Z}(.)$ as long as we can specify a value of $\varepsilon$ such that

$$\begin{aligned} \tilde{Z}(f^*) \le \tilde{Z}(\tilde{f}^*) + \varepsilon \quad \forall f^* \in \mathcal {F}^* \text { and }\tilde{f}^* \in \tilde{\mathcal {F}}^*. \end{aligned}$$

(8)

Here, $f^*$ denotes an optimal classifier to the original problem from the set $\mathcal {F}^* = {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f)$, and $\tilde{f}^*$ denotes an optimal classifier to the surrogate problem from the set $\tilde{\mathcal {F}}^* = {{\mathrm{argmin}}}_{f\in \tilde{\mathcal {F}}} \tilde{Z}(f)$. The width of the the surrogate level set $\varepsilon$ is related to the amount of data that will be removed. If $\varepsilon$ is too large, the method will not remove very many examples and will be less helpful for reducing computation (see Fig. 3).

In the first stage of data reduction, we solve the surrogate problem to: (i) compute the upper bound on the objective value of classifiers in the surrogate level set $\tilde{Z}(\tilde{f}^*)+\varepsilon$; and (ii) to identify a baseline label $\tilde{y}_i := {\text {sign}}\left( \tilde{f}^*(\varvec{x}_i)\right)$ for each example $i = 1,\ldots ,N$. In the second stage of data reduction, we solve a variant of the surrogate problem for each example $i = 1,\ldots ,N$. The ith variant of the surrogate problem includes an additional constraint that forces example i to be classified as $-\tilde{y}_i$:

$$\begin{aligned} \min _{f} ~ \tilde{Z}(f;\mathcal {D}_N) ~\text {s.t.}~ f\in \tilde{\mathcal {F}} \text { and } \tilde{y}_i f(\varvec{x}_i) < 0 \end{aligned}$$

(9)

We denote the optimal classifier to the ith variant as $\tilde{f}^*_{\text {-}i}$. If $\tilde{f}^*_{\text {-}i}$ lies outside of the surrogate level set (i.e., $\tilde{Z}(\tilde{f}^*_{\text {-}i})>\tilde{Z}(\tilde{f}^*)+\varepsilon$) then no classifier in the surrogate level set will label example i as $-\tilde{y}_i$. In other words, all classifiers in the surrogate level set must label this example as $\tilde{y}_i$. Since the surrogate level set contains the optimal classifiers to the original problem by the assumption in (8), we can therefore remove example i from the reduced dataset $\mathcal {D}_M$ because we know that an optimal classifier to the original problem will label this point as $\tilde{y}_i$. We illustrate this situation in Fig. 2.

In Theorem 5, we prove that we obtain the same set of optimal classifiers if we train a model with the initial data $\mathcal {D}_N$ or the reduced data $\mathcal {D}_M$. In Theorem 6, we provide sufficient conditions for a surrogate loss function to satisfy the level set condition in (8).

Fig. 2

We initialize data reduction with $\varepsilon$ large enough so that $\tilde{Z}(f^*)<\tilde{Z}(\tilde{f}^*)+\varepsilon$ for all $f^* \in \mathcal {F}^*$ and all $\tilde{f}^* \in \tilde{\mathcal {F}}^*$. Here, $f^*$ is the optimal classifer to the original problem from the set of optimal classifiers $\mathcal {F}^*$, and $\tilde{f}^*$ is the optimal classifier to the surrogate problem from the set of optimal classifiers $\tilde{\mathcal {F}}^*$. Data reduction fits a classifier $\tilde{f}^*_{\text {-}i}$ for each example in the initial training data $\mathcal {D}_N$ by solving a variant of the surrogate problem with an additional constraint that forces $\tilde{f}^*_{\text {-}i}$ to classify i in a different way than $\tilde{f}^*$. If $\tilde{Z}(\tilde{f}^*_{\text {-}i})>\tilde{Z}(\tilde{f}^*)+\varepsilon$, then we know the predicted class of example i under $f^*$ and can remove it from the reduced training data $\mathcal {D}_M$

Theorem 5

(Equivalence of the Reduced Data) Consider an optimization problem to train a classifier $f \in \mathcal {F}$ with data $\mathcal {D}_N$,

$$\begin{aligned} \min _{f} ~ Z(f;\mathcal {D}_N) ~s.t.~ f\in \mathcal {F}, \end{aligned}$$

as well as a surrogate optimization problem to train a classifier $f \in \tilde{\mathcal {F}}$ with data $\mathcal {D}_N$,

$$\begin{aligned} \min _{f} ~ \tilde{Z}(f;\mathcal {D}_N) ~s.t.~ f\in \tilde{\mathcal {F}}. \end{aligned}$$

Let $f^* \in \mathcal {F}^* := {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_N)$ and $\tilde{f} \in \tilde{\mathcal {F}}^* := {{\mathrm{argmin}}}_{f\in \tilde{\mathcal {F}}} \tilde{Z}(f;\mathcal {D}_N)$. If we choose a value of $\varepsilon$ so that

$$\begin{aligned} \tilde{Z}(f^*;\mathcal {D}_N)&\le \tilde{Z}(\tilde{f}^*;\mathcal {D}_N) + \varepsilon \quad \forall f^* \in \mathcal {F}^* \text { and } \tilde{f}^* \in \tilde{\mathcal {F}}^*, \end{aligned}$$

(10)

then Algorithm 1 will output a reduced dataset $\mathcal {D}_M \subseteq \mathcal {D}_N$ such that

$$\begin{aligned} {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_N) = {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_M). \end{aligned}$$

(11)

Proof

See Appendix.

Theorem 6

(Sufficient conditions to satisfy the level set condition) Consider an optimization problem where the objective minimizes the 0–1 loss function $Z_{01}: \mathbb {R}^P\rightarrow \mathbb {R}$,

$$\begin{aligned} \min _{\varvec{\lambda }\in \mathbb {R}^P} Z_{01}\left( \varvec{\lambda }\right) , \end{aligned}$$

as well as a surrogate optimization problem where the objective minimizes a surrogate loss function $\psi :\mathbb {R}^P\rightarrow \mathbb {R}$,

$$\begin{aligned} \min _{\varvec{\lambda }\in \mathbb {R}^P} Z_{\psi }\left( \varvec{\lambda }\right) . \end{aligned}$$

If the surrogate loss function $\psi$ satisfies the following properties for all $\varvec{\lambda }\in \mathbb {R}^P$, $\varvec{\lambda }^{*}_{01}\in {{\mathrm{argmin}}}_{\varvec{\lambda }\in \mathbb {R}^P} Z_{01}\left( \varvec{\lambda }\right)$, and $\varvec{\lambda }^{*}_{\psi }\in {{\mathrm{argmin}}}_{\varvec{\lambda }\in \mathbb {R}^P} Z_{\psi }\left( \varvec{\lambda }\right)$:

• I. Upper bound on the 0–1 loss: $Z_{01}\left( \varvec{\lambda }\right) \le Z_{\psi }\left( \varvec{\lambda }\right)$

• II. Lipschitz near $\varvec{\lambda }^{*}_{01}$: $\Vert \varvec{\lambda }- \varvec{\lambda }^{*}_{\psi }\Vert < A \implies Z_{\psi }\left( \varvec{\lambda }\right) - Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) < L \Vert \varvec{\lambda }- \varvec{\lambda }^{*}_{\psi }\Vert$

• III. Curvature near $\varvec{\lambda }^{*}_{\psi }$: $\Vert \varvec{\lambda }- \varvec{\lambda }^{*}_{\psi }\Vert > C_{\varvec{\lambda }} \implies Z_{\psi }\left( \varvec{\lambda }\right) - Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) > C_{\psi }$

• IV. Closeness of loss near $\varvec{\lambda }^{*}_{01}$: $| Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) - Z_{01}\left( \varvec{\lambda }^{*}_{01}\right) | < \varepsilon$

then it will also satisfy a level-set condition required for data reduction,

$$\begin{aligned} Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) \le Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) +\varepsilon \quad \forall \varvec{\lambda }^{*}_{01}\text { and } \varvec{\lambda }^{*}_{\psi }, \end{aligned}$$

whenever $\varepsilon = LC_{\varvec{\lambda }}$ obeys $C_{\psi } > 2\varepsilon$.

Proof

See Appendix. $\square$

5.2 Off-the-shelf data reduction for SLIM

Data reduction can easily be applied to SLIM by using an off-the-shelf approach where we use the LP relaxation of the SLIM IP as the surrogate problem.

When we use the LP relaxation to the SLIM IP as the surrogate problem, we can determine a suitable width for the surrogate level set $\varepsilon$ by using a feasible solution to the SLIM IP. To see this, let us denote the SLIM IP as $\min _f Z(f) \;{\text {s.t.}}\; f \in \mathcal {F}$, and denote its LP relaxation as $\min _f Z(f) \;{\text {s.t.}}\; f \in \tilde{\mathcal {F}}$. In addition, let us denote the optimal solution to the SLIM IP as $f^*$ and the optimal solution to the LP relaxation as $\tilde{f}^*$. Since $\mathcal {F}\subseteq \tilde{\mathcal {F}}$, we have that $Z(\tilde{f}^*) \le Z(f^*)$. For any feasible solution to the SLIM IP $\hat{f} \in \mathcal {F}$, we also have that $Z(f^*) \le Z(\hat{f})$. Combining both inequalities, we see that,

$$\begin{aligned} Z(\tilde{f}^*) \le Z(f^*) \le Z(\hat{f}). \end{aligned}$$

Thus, we can satisfy the level set condition (8) using a feasible solution to the SLIM IP $\hat{f} \in \mathcal {F}$ by setting the width of the surrogate level set as,

$$\begin{aligned} \varepsilon (\hat{f}) := Z(\hat{f}) - Z(\tilde{f}^*). \end{aligned}$$

In Fig. 3, we show much training data can be discarded using off-the-shelf data reduction when we train a SLIM scoring system on the bankruptcy dataset (see Table 4). Specifically, we plot the percentage of data removed by Algorithm 1 for values of $\varepsilon \in [\varepsilon _{\min },\varepsilon _{\max }]$ where $\varepsilon _{\min }$ and $\varepsilon _{\max }$ represent the smallest and largest widths of the surrogate level set that could be used in practice. In particular, $\varepsilon _{\min }$ is computed using the optimal solution to the IP as:

$$\begin{aligned} \varepsilon _{\min } := Z(f^*) - Z(\tilde{f}^*), \end{aligned}$$

and $\varepsilon _{\max }$ is computed using a feasible solution to the IP that can be guessed without any computation (i.e., a linear classifier with coefficients $\varvec{\lambda }= \mathbf 0$):

$$\begin{aligned} \varepsilon _{\max } := Z(\mathbf 0 ) - Z(\tilde{f}^*). \end{aligned}$$

In this case, we can discard over 40 % of the training data by using the trivial solution $\varvec{\lambda }= \mathbf 0$, and discard over 80 % of the training data by using a higher quality feasible solution.

Fig. 3

Proportion of training data filtered as a function of the width of the level set, $\varepsilon$ for the bankruptcy dataset. Here, the original problem is an instance of the SLIM IP with $C_0 = 0.01$ and $\mathcal {L}= \{-10,\ldots ,10\}^{P+1}$

5.3 Discussion and related work

Data reduction is a technique that can be applied to a wide range of supervised classification methods, including methods that minimize the 0–1 loss function.

Data reduction is fundamentally different from many techniques for improving the scalability of 0–1 loss minimization, such as oscillation heuristics (Rubin 1990), decomposition (Rubin 1997), cutting planes (Brooks 2011) and specialized branch-and-bound algorithms (Nguyen and Sanner 2013). In fact, data reduction is most similar to the scenario reduction methods in the stochastic programming literature (see e.g., Dupačová et al. 2000, 2003). In comparison to scenario reduction techniques, data reduction does not require the objective function to satisfy stability or convexity assumptions, and is designed to recover the true optimal solution as opposed to an $\varepsilon$-optimal solution.

Data reduction has the advantage that it easily be applied to SLIM by using SLIM’s LP relaxation as the surrogate problem. This off-the-shelf approach may be used as a preliminary procedure before the training process, or as an iterative procedure that is called by the IP solver during the training process as feasible solutions are found. Off-the-shelf data reduction makes use of the integrality gap to identify examples that have to be classified a certain way. Accordingly, the effectiveness of this approach may be improved using techniques that narrow the integrality gap—specifically by using higher quality feasible solutions and/or strengthening SLIM’s LP relaxation (e.g., by using the cutting planes of Brooks 2011).

6 Application to sleep apnea screening

In this section, we discuss a collaboration with the MGH Sleep Laboratory where we used SLIM to create a scoring system for sleep apnea screening (see also Ustun et al. 2015). Our goal is to highlight the flexibility and performance of our approach on a real-world problem that requires a tailored prediction model.

6.1 Data and operational constraints

The dataset for this application contains $N = 1922$ records of patients and $P = 33$ binary features related to their health and sleep habits. Here, $y_i=+1$ if patient i has obstructive sleep apnea (OSA), and there is significant class imbalance as Pr$(y_i=+1) = 76.9\,\%$.

To ensure that the scoring system we produced would be used and accepted by physicians, our collaborators specified three simple operational constraints:

• 1. Limited FPR The model had to achieve the highest possible true positive rate (TPR) while maintaining a maximum false positive rate (FPR) of 20 %. This would ensure that the model could diagnose as many cases of sleep apnea as possible but limit the number of faulty diagnoses.

• 2. Limited model size The model had to be transparent and use at most 5 features. This would ensure that the model was could be explained and understood by other physicians in a short period of time.

• 3. Sign constraints The model had to obey established relationships between well-known risk factors and the incidence of sleep apnea (e.g. it could not suggest that a patient with hypertension had a lower risk of sleep apnea since hypertension is a positive risk factor for sleep apnea).

6.2 Training setup and model selection

We trained a SLIM scoring system with integer coefficients between $-$10 and 10. We addressed all operational constraints without parameter tuning or model selection, as follows:

• We added a loss constraint using the loss variables to limit the maximum FPR at 20 %. We then set $W^{+}= N^{-}/(1+N^{-})$ to guarantee that the optimization would yield a classifier with the highest possible TPR with a maximum FPR less than 20 % (see Sect. 3.1).

• We added a feature-based constraint using the loss variables to limit the maximum number of features to 5 (see Sect. 3.2). We then set $C_0 = 0.9W^{-}/NP$ so that the optimization would yield a classifier that did not sacrifice accuracy for sparsity (see Remark 3).

• We added sign constraints to the coefficients to ensure that our model would not violate established relationships between features and the predicted outcome (i.e., we set $\lambda _j \ge 0$ if there had to be a positive relationship, and $\lambda _j \le 0$ if there had to be a negative relationship).

With this setup, we trained 10 models with subsets of the data to assess predictive accuracy via tenfold cross validation (10-CV), and 1 final model with all of data to hand over to our collaborators. We solved each IP for 1 hour, in parallel, on 12-core 2.7 GHz machine with 48GB RAM. Thus, the training process for SLIM required 1 hour of computing time.

As a comparison, we trained models with 8 baseline classification methods shown in Table 1. We dealt with the class imbalance by using a cost-sensitive approach, where we used a weighted loss function and varied its sensitivity parameter $W^{+}$ across a large range. When possible, we addressed the remaining operational constraints by searching over a fine grid of free parameters. Model selection was difficult for baseline methods because they could not accommodate operational constraints in the same way as SLIM. For each baseline method, we chose the best model that satisfied all operational constraints by: (i) dropping any instance of the free parameters where operational constraints were violated; (ii) choosing the instance that maximized the 10-CV mean test TPR. We ruled that an instance of the free parameters violated an operational constraint if any of the following conditions were met: (1) the 10-CV mean test FPR of the model produced with the instance was greater than the 10-CV mean test FPR of the SLIM model (20.9 %); (2) the model size ² of the final model produced with the instance was greater than 5; (3) the final model produced did not obey sign constraints. This model selection procedure may have biased the results in favor of the baseline methods because we mixed testing and training data by looking at the final model to ensure that operational constraints were satisfied.

Table 1 Training setup for all methods

6.3 Results and observations

In what follows, we report our observations related to operational constraints, predictive performance and interpretability. We show the performance of the best model from each method in Table 2, and summarize the operational constraints they satisfied in Table 3.

Table 2 TPR, FPR and model size for all methods

6.3.1 On the difficulties of handling operational constraints

Among the 9 classification methods that we used, only SLIM, Lasso and Elastic Net could produce a model that satisfied all of operational constraints given to us by physicians. Tree and rule-based methods such as CART, C5.0 Tree and C5.0 Rule were unable to produce a model with a maximum FPR of 20 % (see Fig. 4). Methods that used ${\ell _{2}}$-regularization such as Ridge, SVM Lin. and SVM RBF were unable to produce a model with the required level of sparsity. While we did not expect all methods to satisfy all of the operational constraints, we included them to emphasize the following important points. Namely, state-of-the-art methods for applied predictive modeling do not:

• Handle simple operational constraints that are crucial for models to be used and accepted. Implementations of popular classification methods do not have a mechanism to adjust important model qualities. That is, there is no mechanism to control sparsity in C5.0T (Kuhn et al. 2012) and no mechanism to incorporate sign constraints in SVM (Meyer et al. 2012). Finding a method with suitable controls is especially difficult when a model has to satisfy multiple operational constraints.

• Have controls that are easy-to-use and/or that work correctly. When a method has suitable controls to handle operational constraints, producing a model often requires a tuning process over a high-dimensional free parameter grid. Even after extensive tuning, however, it is possible to never find a model that satisfies all operational constraints (e.g. CART, C5.0R, C5.0T for the Max FPR constraint in Fig. 4).

• Allow tuning to be portable when the training set changes. Consider a standard K-CV model selection procedure where we choose free parameters to maximize predictive accuracy. To adapt this procedure on a problem with operational constraints, we would train models on K validation folds for each instance of the free parameters, choose an instance of the free parameters that maximizes the mean K-CV test accuracy among the instances that satisfied all operational constraints, and train a “final model” for this instance. Unfortunately, there is no guarantee that the final model will obey all operational constraints.

Table 3 Percentage of instances that fulfilled operational constraints

Fig. 4

10-CV mean test FPR for models trained with CART, C5.0, C5.0T across the full range of $W^{+}$. These methods cannot produce a model that satisfies the max FPR $\le$ 20 % constraint

6.3.2 On the sensitivity of acceptable models

Among the three methods that produced acceptable models, the scoring system produced by SLIM had significantly higher sensitivity than the linear models produced by Lasso and Elastic Net—a result that we expected given that SLIM minimizes the 0–1 loss and an ${\ell _{0}}$-penalty while Lasso and Elastic Net minimize convex surrogates of these quantities. This result held true even when we relaxed various operational constraints. In Fig. 5, for instance, we plot the sensitivity and sparsity of models that satisfied the max FPR and sign constraints. Here, we see that Lasso and Elastic Net need at least 8 coefficients to produce a model with the same degree of sensitivity as SLIM. In Fig. 6, we plot the TPR and FPR of models that satisfied the sign and model size constraints. As shown, SLIM scoring systems dominate Lasso and Elastic Net models across the entire ROC curve. These sensitivity advantages are also evident in Table 2: in particular, SLIM yields a model with a similar level of sensitivity and specificity as Ridge and SVM Lin. even as it is fitting models from a far smaller hypothesis space (i.e. linear classifiers with 5 features, sign constraints and integer coefficients vs. linear classifiers with real coefficients).

Fig. 5

Sensitivity and model size of Lasso and Elastic Net models that satisfy the sign and FPR constraints. For each method, we plot the instance that attains the highest 10-CV mean test TPR at model sizes between 0 and 8. Lasso and Elastic Net need at least 8 coefficients to produce a model with the same sensitivity as SLIM

Fig. 6

ROC curve for SLIM, Lasso and Elastic Net instances that satisfy the sign and model size constraints. For each method, we plot the instance that attains the highest 10-CV mean test TPR for 10-CV mean FPR values of $5,10,\ldots ,95\,\%$. Note that we had to train 19 additional instances of SLIM to create this plot

6.3.3 On the usability and interpretability of acceptable models

We include a head-to-head comparison between the most sensitive models that satisfied all operational constraints in Fig. 7 and show the SLIM model as a scoring system in Fig. 8.

Fig. 7

Score functions of the most sensitive predictive models that satisfied all three operational constraints. The baseline models have very poor sensitivity as shown in Table 2

Fig. 8

SLIM scoring system for sleep apnea screening. This model achieves a 10-CV mean test TPR/FPR of 61.4/20.9 %, obeys all operational constraints, and was trained without parameter tuning. It also generalizes well due to the simplicity of the hypothesis space: here the training TPR/FPR of the final model is 62.0/19.6 %

Our collaborators commented that all three models were, in theory, usable by physicians because they were aligned with domain knowledge. Specifically, all models obeyed sign constraints and had large coefficients for well-known risk factors for sleep apnea such as bmi, female, age, snoring and/or hypertension. Unfortunately, the Lasso and Elastic Net models could not realistically be used as screening tools due to their poor sensitivity (29.3 % for Lasso and 44.2 % for Elastic Net). This was not the case for the SLIM model, which had a much higher sensitivity (61.4 %).

Table 4 Datasets used in the numerical experiments

Our results highlight some of the unique interpretability benefits of SLIM scoring systems—that is, their ability to provide “a qualitative understanding of the relationship between joint values of the input variables and the resulting predicted response value” (Hastie et al. 2009). SLIM scoring systems are well-suited to provide this kind of qualitative understanding due to their high level of sparsity and small integer coefficients. These qualities help users gauge the influence of each input variable with respect to the others, which is especially important because humans can only handle a few cognitive entities at once ($7\pm 2$ according to Miller 1984), and are seriously limited in estimating the association between three or more variables (Jennings et al. 1982). Accordingly, these qualities may also help users gauge the influence Sparsity and small integer coefficients also allow users to make quick predictions without a computer or a calculator, which may help them understand how the model works by actively using it to classify prototypical examples. Here, this process helped our collaborators come up with the following simple rule-based explanation for our model predicted that a patient has OSA (i.e., when SCORE $>$ 1): “if the patient is male, predict OSA if age $\ge$ 60 OR hypertension OR bmi $\ge$ 30; if the patient is female, predict OSA if bmi $\ge$ 40 AND (age $\ge$ 60 OR hypertension).”

7 Numerical experiments

In this section, we present numerical experiments to compare the accuracy and sparsity of SLIM scoring systems to other popular classification models. Our goal is to illustrate the off-the-shelf performance of SLIM and show that we can train accurate scoring systems for real-sized datasets in minutes.

7.1 Experimental setup

Datasets: We ran numerical experiments on 8 datasets from the UCI Machine Learning Repository (Bache and Lichman 2013) summarized in Table 4. We chose these datasets to explore the performance of each method as we varied the size and nature of the training data. We processed each dataset by binarizing all categorical features and some real-valued features. For the purposes of reproducibility, we include all processed datasets in Online Resource 1.

Methods: We summarize the training setup for each method in Table 5. We trained SLIM scoring systems using the CPLEX 12.6.0.0 API and models with baseline methods using publicly available packages in R 3.1.1 (R Core Team 2014). For each method, each dataset, and each unique combination of free parameters, we trained 10 models using subsets of the data to estimate predictive accuracy via tenfold cross-validation (10-CV), and 1 final model using all of the data to assess sparsity and interpretability. We ran all baseline methods without time constraints over a large grid of free parameters. We produced an $\ell _0$-regularization path for SLIM by solving $6 \times 11$ IPs for each dataset (6 values of $C_0$ $\times$ 11 training runs per $C_0$). We allocated at most 10 min to solve each IP, and solved 12 IPs in parallel on a 12-core 2.7 GHz machine with 48 GB RAM. Thus, it took at most 1 hour to train SLIM scoring systems for each dataset. Since the adult and haberman datasets were imbalanced, we trained all methods on these datasets with a weighted loss function where we set $W^{+}= N^{-}/N$ and $W^{-}= N^{+}/N$.

Table 5 Training setup for classification methods used for the numerical experiments

Table 6 Accuracy and sparsity of all methods on all datasets

7.2 Results and observations

We summarize the results of our experiments in Table 6 and Figs. 13, 14. We report the sparsity of models using a metric that we call model size. Model size represents the number of coefficients for linear models (Lasso, Ridge, Elastic Net, SLIM, SVM Lin.), the number of leaves for decision tree models (C5.0T, CART), and the number of rules for rule-based models (C5.0R). For completeness, we set the model size for black-box models (SVM RBF) to the number of features in each dataset.

We show the accuracy and sparsity of all methods on all dataset in Figs. 13 and 14. For each dataset, and each method, we plot a point at the 10-CV mean test error and final model size, and surround this point with an error bar whose height corresponds to the 10-CV standard deviation in test error. In addition, we include $\ell _0$-regularization paths for SLIM and Lasso on the right side of Figs. 13 and 14 to show how the test error varies at different levels of sparsity for sparse linear models.

7.2.1 On accuracy, sparsity, and computation

Our results show that many methods are unable to produce models that attain the same levels of accuracy and sparsity as SLIM. As shown in Figs. 13 and 14, SLIM always produces a model that is more accurate than Lasso at some level of sparsity, and sometimes more accurate at all levels of sparsity (e.g., spambase, haberman, mushroom, breastcancer). Although the optimization problems we solved to train SLIM scoring systems were $\mathcal {NP}$-hard, we did not find any evidence that computational issues hurt the performance of SLIM on any of the datasets. We obtained accurate and sparse models for all datasets in 10 minutes using CPLEX 12.6. Further, the solver provided a proof of optimality (i.e. a relative MIPGAP of 0.0 %) for all scoring systems we trained for mammo, mushroom, bankruptcy, breastcancer.

7.2.2 On the regularization effect of discrete coefficients

We expect that methods that directly optimize accuracy and sparsity will achieve the best possible accuracy at every level of sparsity (i.e. the best possible trade-off between accuracy and sparsity). SLIM directly optimizes accuracy and sparsity. However, it may not necessarily achieve the best possible accuracy at each level of sparsity because it restricts coefficients to a finite discrete set $\mathcal {L}$.

By comparing SLIM to Lasso, we can identify a baseline regularization effect due to this $\mathcal {L}$ set restriction. In particular, we know that when Lasso’s performance dominates that of SLIM, it is very arguably due to the use of a small set of discrete coefficients. Our results show that this tends to happen mainly at large model sizes (see e.g., the regularization path for breastcancer, heart, mammo). This suggests that the $\mathcal {L}$ set restriction has a more noticeable impact on accuracy at larger model sizes.

One interesting effect of the $\mathcal {L}$ set restriction is that the most accurate SLIM scoring system may not use all of the features in the dataset. In our experiments, we always trained SLIM with $C_0 = 0.9/NP$ to obtain a scoring system with the highest training accuracy among linear models with coefficients in $\varvec{\lambda }\in \mathcal {L}$ (see Remark 3). In the bankruptcy dataset, for example, we find that this model only uses 3 out of 6 features. This is due to the $\mathcal {L}$ set restriction: if the $\mathcal {L}$ restriction were relaxed, then the method would use all features to improve its training accuracy (as is the case with Ridge or SVM Lin.).

Fig. 9

SLIM scoring system for mushroom. This model has a 10-CV mean test error of $0.0\pm 0.0\,\%$

Fig. 10

Lasso score function for mushroom. This model has a 10-CV mean test error of $0.0 \pm 0.0\,\%$

Fig. 11

C5.0 decision tree for mushroom. This model has a 10-CV mean test error of $0.0\pm 0.0\,\%$

Fig. 12

C5.0 rule list for mushroom. This model has a 10-CV mean test error of $0.0\pm 0.0\,\%$

7.2.3 On interpretability

To discuss interpretability, we focus on the mushroom dataset, which provides a nice basis for comparison as many methods produce a model that attains perfect predictive accuracy. In Figs. 9, 10, 11 and 12, we show the sparsest models that achieve perfect predictive accuracy on the mushroom dataset. We omit models from some methods because they do not attain perfect accuracy (CART) or use far more features (Ridge, SVM Lin, SVM RBF) (Figs. 13, 14).

Fig. 13

Accuracy and sparsity of all classification methods on all datasets. For each dataset, we plot the performance of models when free parameters are set to values that minimize the 10-CV mean test error (left), and plot the performance of SLIM and Lasso across the full ${\ell _{0}}$-regularization path (right)

Fig. 14

Accuracy and sparsity of all classification methods on all datasets. For each dataset, we plot the performance of models when free parameters are set to values that minimize the 10-CV mean test error (left) and plot the performance of SLIM and Lasso across the full ${\ell _{0}}$-regularization path (right)

In this case, the SLIM scoring system uses 7 integer coefficients. However, it can be expressed as a 5 line scoring system due to the fact that odor=none, odor=almond, and odor=anise are mutually exclusive variables that all use the same coefficient. The model benefits from the fact that it not only allows users make predictions by hand but uses a linear form that helps users to gauge the influence of each input variable with respect to the others. We note that only some of these qualities are found in other models. The Lasso model, for instance, has a linear form but uses far more features. Similarly, the C5.0 models allow users to make predictions by hand, but use a hierarchical structure that makes it difficult to assess the influence of each input variable with respect to the others (see Freitas 2014, for a discussion).

We believe that these qualities represent “baseline” interpretability benefits of SLIM scoring systems. In practice, interpretability is a subjective and multifaceted notion (i.e., it depends on who will be using the model but also depends on multiple model qualities as discussed by Kodratoff 1994; Pazzani 2000; Freitas 2014). In light of this fact, SLIM has the unique benefit in that it allows practitioners to work closely with the target audience and directly encode all interpretability requirements by means of operational constraints.

8 Specialized models

In this section, we present three specialized models related to SLIM. These models are all special instances of the optimization problem in (1).

8.1 Personalized models

A Personalized Integer Linear Model (PILM) is a generalization of SLIM that provides soft control over the coefficients in a scoring system. To use this model, users define $R+1$ interpretability sets,

$$\begin{aligned} \mathcal {L}_r = \{ l_{r,1},\ldots ,l_{r,K_r}\}\quad \hbox { for }\quad r = 0,\ldots , R, \end{aligned}$$

as well as a “personalized” interpretability penalty,

$$\begin{aligned} \varPhi _j(\lambda _j)&= {\left\{ \begin{array}{ll} C_0 &{} if\; \lambda _j \in \mathcal {L}_0 \\ &{} \vdots \\ C_R &{} if\; \lambda _j \in \mathcal {L}_R. \end{array}\right. } \end{aligned}$$

In order to penalize coefficients from less interpretable sets more heavily, we need that: (i) $\mathcal {L}_1,\ldots ,\mathcal {L}_R$ are mutually exclusive; (ii) $\mathcal {L}_r$ is more interpretable than $\mathcal {L}_{r+1}$; (iii) the trade-off parameters are monotonically increasing in r, so that $C_0<\cdots <C_R$. The values of the parameters $C_r$ can be set as the minimum gain in training accuracy required for the optimal classifier to use a coefficient from $\mathcal {L}_r$.

As an example, consider training a PILM scoring system with the penalty:

$$\begin{aligned} \varPhi _j(\lambda _j)&= {\left\{ \begin{array}{ll} C_0 = 0.00 &{} \text {if} \quad \lambda _j \in {0}\\ C_1 = 0.01 &{} \text {if} \quad \lambda _j \in \pm \{1,\ldots ,10\}\\ C_2 = 0.05 &{} \text {if} \quad \lambda _j \in \pm \{11,\ldots ,100\}. \end{array}\right. } \end{aligned}$$

Here, the optimal classifier will use a coefficient from $\mathcal {L}_1$ if it yields at least a 1% gain in training accuracy, and a coefficient from $\mathcal {L}_2$ if it yields at least a 5% gain in training accuracy.

We can train a PILM scoring system by solving the following IP:

$$\begin{aligned} \min _{\varvec{\lambda },{\varvec{\psi }},{\varvec{\varPhi }},\mathbf {u}}&\frac{1}{N}\sum _{i=1}^{N} \psi _i + \sum _{j=1}^{P} \varPhi _j\nonumber \\ {\text {s.t.}}&\quad M_i\psi _i \ge \gamma -\sum _{j=0}^P y_i \lambda _j x_{i,j}\quad \quad \quad \quad \quad \quad \quad {i}={1}\text {,...,}{N} \quad \quad \quad \quad \quad \quad \ \textit{ 0-1 loss} \end{aligned}$$

(12a)

$$\begin{aligned}&\varPhi _j = \sum _{r=0}^R \sum _{k=1}^{K_r} C_r u_{j,k,r} \ \quad \quad \quad \quad \quad \quad \quad \quad \quad {j}={1}\text {,...,}{P} \quad \quad \quad \quad \quad \ \textit{ int. penalty} \end{aligned}$$

(12b)

$$\begin{aligned}&\lambda _j = \sum _{r=0}^R \sum _{k=1}^{K_r} l_{r,k} u_{j,k,r}\ \quad \quad \quad \quad \quad \ \quad \quad \quad \quad {j}={0}\text {,...,}{P} \quad \quad \quad \textit{ coefficient values} \end{aligned}$$

(12c)

$$\begin{aligned}&1 = \sum _{r=0}^R \sum _{k=1}^{K_r} u_{j,k,r} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {j}={0}\text {,...,}{P} \quad \quad \quad \!\textit{ 1 int. set per coef.}\\&\psi _i \in \{0,1\}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {i}={1}\text {,...,}{N} \quad \quad \quad \, \quad \textit{ loss variables}\nonumber \\&\varPhi _j \in \mathbb {R}_+ \ \,\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {j}={1}\text {,...,}{P} \ \,\quad \textit{ int. penalty variables} \nonumber \\&u_{j,r,k} \in \{0,1\}\quad {j}={0}\text {,...,}{P} \quad {r}={0}\text {,...,}{R} \quad {k}={1}\text {,...,}{K_r} \quad \,\textit{ coef. value variables}\nonumber \end{aligned}$$

(12d)

Here, the loss constraints and Big-M parameters in (12a) are identical to those from the SLIM IP formulation in Sect. 2. The $u_{j,k,r}$ are binary indicator variables that are set to 1 if $\lambda _j$ is equal to $l_{k,r}$. Constraints (12d) ensure that each coefficient uses exactly one value from one interpretability set. Constraints (12c) ensure that each coefficient $\lambda _j$ is assigned a value from the appropriate interpretability set $\mathcal {L}_r$. Constraints (12b) ensure that each coefficient $\lambda _j$ is assigned the value specified by the personalized interpretability penalty.

8.2 Rule-based models

SLIM can be extended to produce specialized “rule-based” models when the training data are composed of binary rules. In general, any real-valued feature can be converted into a binary rule by setting a threshold (e.g., we can convert age into the feature $age \ge 25 := \mathbbm {1}\left[ age \ge 25\right]$). The values of the thresholds can be set using domain expertise, rule mining, or discretization techniques (Liu et al. 2002).

In what follows, we assume that we train models with a binarized dataset that contains $T_j$ binary rules $\varvec{h}_{j,t} \in \{0,1\}^N$ for each feature $\varvec{x}_j \in \mathbb {R}^N$ in the original dataset. Thus, we consider models with the form:

$$\begin{aligned} \hat{y} = {\text {sign}}\left( \lambda _0 + \sum _{j=1}^P \sum _{t=1}^{T_j} \lambda _{j,t}h_{j,t}\right) . \end{aligned}$$

We make the following assumptions about the binarization process. If $\varvec{x}_j$ is a binary variable, then it is left unchanged so that $T_j = 1$ and $\varvec{h}_{j,T_j} := \varvec{x}_j$. If $\varvec{x}_j$ is a categorical variable $\varvec{x}_j \in \{1,\ldots ,K\}$, the binarization yields a binary rule for each category so that $T_j=K$ and $\varvec{h}_{j,t}:=\mathbbm {1}\left[ \varvec{x}_j=k\right]$ for $t = 1,\ldots ,K$. If $\varvec{x}_j$ is a real variable, then the binarization yields $T_j$ binary rules ³ of the form $\varvec{h}_{j,t} := \mathbbm {1}\left[ \varvec{x}_j \ge v_{j,t}\right]$ where $v_{j,t}$ denotes the tth threshold for feature j.

8.2.1 M-of-N rule tables

M-of-N rule tables are simple rule-based models that, given a set of N rules, predict $\hat{y} = +1$ if at least M of them are true (see e.g., Fig. 15). These models have the major benefit that they do not require the user to compute a mathematical expression. M-of-N rule tables were originally proposed as auxiliary models that could be extracted from neural nets (Towell and Shavlik 2013) but can also be trained as stand-alone discrete linear classification models as suggested by Chevaleyre et al. (2013).

Fig. 15

M-of-N rule table for the breastcancer dataset for $C_0 = 0.9/NP$. This model has 8 rules and a 10-CV mean test error of $4.8 \pm 2.5\%$. We trained this model with binary rules $h_{i,j} := \mathbbm {1}\left[ x_{i,j} \ge 3\right]$

We can produce a fully optimized M-of-N rule table by solving an optimization problem of the form:

$$\begin{aligned} \min _{\varvec{\lambda }}&\qquad \sum _{i=1}^N\mathbbm {1}\left[ y_i \hat{y}_i \le 0\right] + C_0 \left\| \varvec{\lambda }\right\| _0&\\ {\text {s.t.}}&\qquad \lambda _0 \in \{-P,\ldots ,0\}&\\&\qquad \lambda _{j,t} \in \{0,1\}&{j}={1}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}. \end{aligned}$$

The coefficients from this optimization problem yield an M-of-N rule table with $M=\lambda _0 + 1$ and $N= \sum _{j=1}^P\sum _{t=1}^{T_j} \lambda _{j,t}$. Here, we can achieve exact ${\ell _{0}}$-regularization using an ${\ell _{1}}$-penalty since $\left\| \lambda _{j,t}\right\| _0 = \left\| \lambda _{j,t}\right\| _1$ for $\lambda _{j,t} \in \{0,1\}$. Since we use the 0–1 loss, the trade-off parameter $C_0$ can be set as minimum gain in training accuracy required to include a rule in the optimal table.

We can train an M-of-N rule table by solving the following IP:

$$\begin{aligned} \min _{\varvec{\lambda },{\varvec{\psi }},{\varvec{\varPhi }}}&\frac{1}{N}\sum _{i=1}^{N} \psi _i + \sum _{j=1}^{P} \varPhi _j \nonumber \\ {\text {s.t.}}&\quad M_i \psi _i \ge \gamma -\sum _{j=0}^P\sum _{t=1}^{T_j} y_i \lambda _{j,t} h_{i,j,t}&\quad {i}={1}\text {,...,}{N}&\quad \textit{ 0-1 loss} \end{aligned}$$

(13a)

$$\begin{aligned}&\quad \varPhi _{j,t} = C_0 \lambda _{j,t}&{j}={1}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}&\quad \textit{ int. penalty} \nonumber \\&\quad \lambda _0 \in \{-P,\ldots ,0\}&\quad \quad&\quad \textit{ intercept value} \nonumber \\&\quad \lambda _{j,t} \in \{0,1\}&{j}={1}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}&\quad \textit{ coefficient values} \nonumber \\&\quad \psi _i \in \{0,1\}&{i}={1}\text {,...,}{N}&\quad \textit{ 0-1 loss indicators} \nonumber \\&\quad \varPhi _{j,t} \in \mathbb {R}_+&{j}={1}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}&\quad \textit{ int. penalty values} \end{aligned}$$

(13b)

Here, the loss constraints and Big-M parameters in (13a) are identical to those from the SLIM IP formulation in Sect. 2. Constraints (13b) define the penalty variables $\varPhi _{j,t}$ as the value of the ${\ell _{0}}$-penalty.

8.2.2 Threshold-rule models

A Threshold-Rule Integer Linear Model (TILM) is a scoring system where the input variables are thresholded versions of the original feature set (i.e. decision stumps). These models are well-suited to problems where the outcome has a non-linear relationship with real-valued features. As an example, consider the SAPS II scoring system of Le Gall et al. (1993), which assesses the mortality of patients in intensive care using thresholds on real-valued features such as $blood\_pressure>200$ and $heart\_rate<40$. TILM optimizes the binarization of real-valued features by using feature selection on a large (potentially exhaustive) pool of binary rules for each real-valued feature. Carrizosa et al. (2010), Belle et al. (2013) and Goh and Rudin (2014) take different but related approaches for constructing classifiers with binary threshold rules.

We train TILM scoring systems using an optimization problem of the form:

$$\begin{aligned} \min _{\varvec{\lambda }}&\qquad \frac{1}{N}\sum _{i=1}^N\mathbbm {1}\left[ y_i \hat{y}_i \le 0\right] + C_f \cdot \text {Features} + C_t \cdot \text {Rules per Feature} + \epsilon \left\| \varvec{\lambda }\right\| _1 \\ {\text {s.t.}}&\qquad \varvec{\lambda }\in \mathcal {L}, \\&\qquad \sum _{t=1}^{T_j} \mathbbm {1}\left[ \lambda _{j,t}\ne 0\right] \le R_{max}\quad \text { for }\quad j = 1,\ldots ,P, \\&\qquad {\text {sign}}\left( \lambda _{j,1}\right) = \cdots = {\text {sign}}\left( \lambda _{j,T_j}\right) \quad \text { for }\quad j = 1,\ldots ,P. \end{aligned}$$

TILM uses an interpretability penalty that penalizes the number of rules used in the classifier as well as the number of features associated with these rules. The small ${\ell _{1}}$-penalty in the objective restricts coefficients to coprime values as in SLIM. Here, $C_f$ tunes the number of features used in the model, $C_t$ tunes the number of rules per feature, and $\epsilon$ is set to a small value to produce coprime coefficients. TILM includes additional hard constraints to limit the number of rules per feature to $R_{max}$ (e.g., $R_{max} = 3$), and to ensure that the coefficients for binary rules from a single feature agree in sign (this ensures that each feature maintains a strictly monotonically increasing or decreasing relationship with the outcome).

We can train a TILM scoring system by solving the following IP:

$$\begin{aligned} \min _{\varvec{\lambda },{\varvec{\psi }},{\varvec{\varPhi }},{\varvec{\tau }},{\varvec{\nu }},{\varvec{\delta }}}&\frac{1}{N}\sum _{i=1}^{N} \psi _i + \sum _{j=1}^{P} \varPhi _j \nonumber \\ {\text {s.t.}}&\quad M_i \psi _i \ge \gamma -\sum _{j=0}^P\sum _{t=1}^{T_j} y_i \lambda _{j,t} h_{i,j,t}&{i}={1}\text {,...,}{N}&\quad \textit{ 0-1 loss} \end{aligned}$$

(14a)

$$\begin{aligned}&\varPhi _j = C_f\nu _j + C_t \tau _j + \epsilon \sum _{t=1}^{T_j}{\beta _{j,t}}&{j}={1}\text {,...,}{P}&\quad \textit{ int. penalty} \end{aligned}$$

(14b)

$$\begin{aligned}&T_j \nu _j = \sum _{t=1}^{T_j} \alpha _{j,t}&{j}={1}\text {,...,}{P}&\quad \textit{ feature use} \end{aligned}$$

(14c)

$$\begin{aligned}&\tau _j = \sum _{t=1}^{T_j} \alpha _{j,t}-1&{j}={1}\text {,...,}{P}&\quad \textit{ threshold/feature} \end{aligned}$$

(14d)

$$\begin{aligned}&\tau _j \le R_{max} + 1&{j}={1}\text {,...,}{P}&\quad \textit{ max thresholds} \nonumber \\&-\varLambda _j\alpha _{j,t} \le \lambda _{j,t} \le \varLambda _j\alpha _{j,t}&{j}={1}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}&\quad {\ell _{0}}\,\textit{ norm} \nonumber \\&-\beta _{j,t} \le \lambda _{j,t} \le \beta _{j,t}&{j}={1}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}&\quad {\ell _{1}}\,\textit{ norm} \end{aligned}$$

(14e)

$$\begin{aligned}&-\varLambda _j (1-\delta _{j}) \le \lambda _{j,t} \le \varLambda _j\delta _{j}&{j}={1}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}&\quad \textit{ agree in sign} \end{aligned}$$

(14f)

$$\begin{aligned}&\lambda _{j,t} \in \mathcal {L}_j&{j}={0}\text {,...,}{P} \quad {t}={1}\text {,...,}{T_j}&\quad \textit{ coefficient values} \nonumber \\&\psi _i \in \{0,1\}&{i}={1}\text {,...,}{N}&\quad \textit{ 0--1 loss indicators} \nonumber \\&\varPhi _j \in \mathbb {R}_+&{j}={1}\text {,...,}{P}&\quad \textit{ int. penalty variables} \nonumber \\&\alpha _j \in \{0,1\}&{j}={1}\text {,...,}{P}&\quad {\ell _{0}}\, \textit{ variables} \nonumber \\&\beta _j \in \mathbb {R}_+&{j}={1}\text {,...,}{P}&\quad {\ell _{1}}\, \textit{ variables} \nonumber \\&\nu _j \in \{0,1\}&{j}={1}\text {,...,}{P}&\, \textit{ feature use indicators} \nonumber \\&\tau _j \in \mathbb {Z}_+&{j}={1}\text {,...,}{P}&\, \textit{ threshold/feature variables} \nonumber \\&\delta _j \in \{0,1\}&{j}={1}\text {,...,}{P}&\quad \textit{ sign indicators} \end{aligned}$$

(14g)

Here, the loss constraints and Big-M parameters in (14a) are identical to those from the SLIM IP formulation in Sect. 2. Constraints (14b) set the interpretability penalty for each coefficient as $\varPhi _j = C_f\nu _j + C_t \tau _j + \epsilon \sum \beta _{j,t}$. The variables in the interpretability penalty include: $\nu _j$, which indicate that we use at least one threshold rule from feature j; $\tau _j$, which count the number of additional binary rules derived from feature j; and $\beta _{j,t}:= |\lambda _{j,t}|$. The values of $\nu _j$ and $\tau _j$ are set using the indicator variables $\alpha _{j,t} := \mathbbm {1}\left[ \lambda _{j,t}\ne 0\right]$ in constraints (14c) and (14d). Constraints (14e) limit the number of binary rules from feature j to $\mathbb {R}_{max}$. Constraints (14g) ensure that the coefficients of binary rules derived from feature j agree in sign; these constraints are encoded using the variables $\delta _j:=\mathbbm {1}\left[ \lambda _{j,t} \ge 0\right]$.

9 Conclusion

In this paper, we introduced a new method for creating data-driven medical scoring systems which we refer to as a Supersparse Linear Integer Model (SLIM). We showed how SLIM can produce scoring systems that are fully optimized for accuracy and sparsity, that can accommodate multiple operational constraints, and that can be trained without parameter tuning.

The major benefits of our approach over existing methods come from the fact that we avoid approximations that are designed to achieve faster computation. Approximations such as surrogate loss functions and ${\ell _{1}}$-regularization hinder the accuracy and sparsity of models as well as the ability of practitioners to control these qualities. Such approximations are no longer needed for many datasets, since using current integer programming software, we can now train scoring systems for many real-world problems. Integer programming software also caters to practitioners in other ways, by allowing them to choose from a pool of models by mining feasible solutions and to seamlessly benefit from periodic computational improvements without revising their code.

A Proofs of Theorems

1.1 Proof of Theorem 1 (Minimum margin resolution bound)

Proof

We use normalized versions of the vectors, $\varvec{\rho }/\left\| \varvec{\rho }\right\| _2$ and $\varvec{\lambda }/\varLambda$ because the 0–1 loss is scale invariant:

$$\begin{aligned} \sum _{i=1}^N \mathbbm {1}\left[ y_i\varvec{\lambda }^T\varvec{x}_i\le 0\right]&= \sum _{i=1}^N \mathbbm {1}\left[ y_i\frac{\varvec{\lambda }^T\varvec{x}_i}{\varLambda }\le 0\right] , \\ \sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\rho }^T \varvec{x}_i \le 0\right]&= \sum _{i=1}^N \mathbbm {1}\left[ y_i\frac{\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2}\le 0\right] . \end{aligned}$$

We set $\varLambda > \frac{X_{\max }\sqrt{P}}{2 \gamma _{\min }}$ as in (5). Using $\varLambda$, we then define $\varvec{\lambda }/\varLambda$ element-wise so that $\lambda _j/\varLambda$ is equal to $\rho _j/\left\| \varvec{\rho }\right\| _2$ rounded to the nearest $1/\varLambda$ for $j= 1,\ldots ,P$.

We first show that our choice of $\varLambda$ and $\varvec{\lambda }$ ensures that the difference between the margin of $\varvec{\rho }/\left\| \varvec{\rho }\right\| _2$ and the margin of $\varvec{\lambda }/\varLambda$ on all training examples is always less than the minimum margin of $\varvec{\rho }/\left\| \varvec{\rho }\right\| _2$, defined as $\gamma _{\min } = \min _i \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2}$. This statement follows from the fact that, for all i:

$$\begin{aligned} \left| \frac{\varvec{\lambda }^T \varvec{x}_i}{\varLambda }-\frac{\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2} \right| ~ \le ~&\left\| \frac{\varvec{\lambda }}{\varLambda } - \frac{\varvec{\rho }}{\left\| \varvec{\rho }\right\| _2}\right\| _2\Vert \varvec{x}_i\Vert _2 \nonumber \\ = ~&\left( \sum _{j=1}^P \left| \frac{\lambda _j}{\varLambda } - \frac{\rho _j}{\left\| \varvec{\rho }\right\| _2} \right| ^2 \right) ^{1/2}\Vert \varvec{x}_i\Vert _2 \end{aligned}$$

(15)

$$\begin{aligned} \le ~&\left( \sum _{j=1}^P \frac{1}{(2\varLambda )^2} \right) ^{1/2}\Vert \varvec{x}_i\Vert _2 \nonumber \\ = ~&\frac{\sqrt{P}}{2\varLambda }X_{\max } \end{aligned}$$

(16)

$$\begin{aligned} < ~&\frac{\sqrt{P}X_{\max }}{2\left( \frac{X_{\max }\sqrt{P}}{2\min _i \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2}} \right) } \end{aligned}$$

(17)

$$\begin{aligned} = ~&\min _i \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2}. \end{aligned}$$

(18)

Here, the inequality in (15) uses the Cauchy–Schwarz inequality; the inequality in (16) is due to the fact that the distance between $\rho _j/\left\| \varvec{\rho }\right\| _2$ and $\lambda _j/\varLambda$ is at most $1/2\varLambda$; and the inequality in (17) is due to our choice of $\varLambda$.

Next, we show that our choice of $\varLambda$ and $\varvec{\lambda }$ ensures that $\varvec{\rho }/\left\| \varvec{\rho }\right\| _2$ and $\varvec{\lambda }/\varLambda$ classify each point in the same way. We consider three cases: first, the case where $\varvec{x}_i$ lies on the margin; second, the case where $\varvec{\rho }$ has a positive margin on $\varvec{x}_i$; and third, the case where $\varvec{\rho }$ has a negative margin on $\varvec{x}_i$. For the case when $\varvec{x}_i$ lies on the margin, $\min _i|\varvec{\rho }^T\varvec{x}_i|=0$ and the theorem holds trivially. For the case where $\varvec{\rho }$ has positive margin, $\varvec{\rho }^T\varvec{x}_i>0$, the following calculation using (18) is relevant:

$$\begin{aligned} \frac{\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2} - \frac{\varvec{\lambda }^T\varvec{x}_i}{\varLambda } \le \left| \frac{\varvec{\lambda }^T\varvec{x}_i}{\varLambda }-\frac{\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2} \right| < \min _i \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2}. \end{aligned}$$

We will use the fact that for any $i^{'}$, by definition of the minimum:

$$\begin{aligned} 0 \le \frac{|\varvec{\rho }^T\varvec{x}_{i^{'}}|}{\left\| \varvec{\rho }\right\| _2}- \min _i \frac{|\varvec{\rho }^T\varvec{x}_i |}{\left\| \varvec{\rho }\right\| _2}, \end{aligned}$$

and combine this with a rearrangement of the previous expression to obtain:

$$\begin{aligned} 0\le \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2} - \min _i \frac{|\varvec{\rho }^T\varvec{x}_i |}{\left\| \varvec{\rho }\right\| _2} = \frac{\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2} - \min _i \frac{|\varvec{\rho }^T\varvec{x}_i |}{\left\| \varvec{\rho }\right\| _2} < \frac{\varvec{\lambda }^T\varvec{x}_i}{\varLambda }. \end{aligned}$$

Thus, we have shown that $\varvec{\lambda }^T\varvec{x}_i>0$ whenever $\varvec{\rho }^T\varvec{x}_i>0$.

For the case where $\varvec{\rho }$ has a negative margin on $\varvec{x}_i$, $\varvec{\rho }^T\varvec{x}_i<0$, we perform an analogous calculation:

$$\begin{aligned} \frac{\varvec{\lambda }^T\varvec{x}_i}{\Vert \varvec{\lambda }\Vert _2} - \frac{\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2} \le \left| \frac{\varvec{\lambda }^T\varvec{x}_i}{\varLambda }-\frac{\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2} \right| < \min _i \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2}. \end{aligned}$$

and then using that $\varvec{\rho }^T\varvec{x}_i<0$,

$$\begin{aligned} 0\le \frac{|\varvec{\rho }^T\varvec{x}_i|}{\left\| \varvec{\rho }\right\| _2} - \min _i \frac{|\varvec{\rho }^T\varvec{x}_i |}{\left\| \varvec{\rho }\right\| _2} = \frac{-\varvec{\rho }^T\varvec{x}_i}{\left\| \varvec{\rho }\right\| _2} - \min _i \frac{|\varvec{\rho }^T\varvec{x}_i |}{\left\| \varvec{\rho }\right\| _2} < -\frac{\varvec{\lambda }^T\varvec{x}_i}{\varLambda }. \end{aligned}$$

Thus, we have shown $\varvec{\lambda }^T\varvec{x}_i<0$ whenever $\varvec{\rho }^T\varvec{x}_i<0$.

Putting both the positive margin and negative margin cases together, we find that for all i,

$$\begin{aligned} \mathbbm {1}\left[ y_i\varvec{\rho }^T\varvec{x}_i\le 0\right] = \mathbbm {1}\left[ y_i\varvec{\lambda }^T\varvec{x}_i\le 0\right] . \end{aligned}$$

Summing over i yields the statement of the theorem. $\square$

1.2 Proof of Theorem 3 (Generalization of sparse discrete linear classifiers)

Proof

Let $Z(\varvec{\lambda };\mathcal {D}_N) = \frac{1}{N} \sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0\left\| \varvec{\lambda }\right\| _0.$ Note that $\varvec{\lambda }= 0$ is a feasible solution since we assume that $0 \in \mathcal {L}$. Since $\varvec{\lambda }= 0$ achieves an objective value of $Z(0;\mathcal {D}_N) = 1$, any optimal solution, $\varvec{\lambda }\in {{\mathrm{argmin}}}_{\lambda \in \mathcal {L}} Z(\varvec{\lambda };\mathcal {D}_N)$, must attain an objective value $Z(\varvec{\lambda };\mathcal {D}_N) \le 1$. This implies

$$\begin{aligned} Z(\varvec{\lambda };\mathcal {D}_N) \le 1,\\ C_0\left\| \varvec{\lambda }\right\| _0 \le \frac{1}{N}\sum _{i=1}^N \mathbbm {1}\left[ y_i \varvec{\lambda }^T \varvec{x}_i \le 0\right] + C_0\left\| \varvec{\lambda }\right\| _0 \le 1,\\ \left\| \varvec{\lambda }\right\| _0 \le \frac{1}{C_0},\\ \left\| \varvec{\lambda }\right\| _0 \le \left\lfloor \frac{1}{C_0} \right\rfloor . \end{aligned}$$

The last line uses that $\Vert \varvec{\lambda }\Vert _0$ is an integer.

Thus, $\mathcal {H}_{P,C_0}$ is large enough to contain all minimizers of $Z(\cdot ;\mathcal {D}_N)$ for any $\mathcal {D}_N$. The statement of the theorem follows from applying Theorem 2. $\square$

1.3 Proof of Theorem 5 (Equivalence of the reduced data)

Proof

Let us denote the set of classifiers whose objective value is less or equal to $\tilde{Z}(\tilde{f}^*;\mathcal {D}_N)$ as

$$\begin{aligned} \tilde{\mathcal {F}}^{\varepsilon } = \left\{ f \in \tilde{\mathcal {F}} \; \Big | \; \tilde{Z}(f;\mathcal {D}_N) \le \tilde{Z}(\tilde{f};\mathcal {D}_N) + \varepsilon \right\} . \end{aligned}$$

In addition, let us denote the set of points that have been removed by the data reduction algorithm

$$\begin{aligned} \mathcal {S} = \mathcal {D}_N \setminus \mathcal {D}_M. \end{aligned}$$

By definition, data reduction only removes an example if its sign is fixed. This means that ${\text {sign}}\left( f(\varvec{x}_i)\right) ={\text {sign}}\left( \tilde{f}(\varvec{x}_i)\right)$ for all $i \in \mathcal {S}$ and $f\in \tilde{\mathcal {F}}^{\varepsilon }.$ Thus, we can see that for all classifiers $f \in \tilde{\mathcal {F}}^{\varepsilon }$,

$$\begin{aligned}&Z(f;\mathcal {D}_N) = Z(f;\mathcal {D}_M)+\sum _{i\in \mathcal {S}} \mathbbm {1}\left[ y_if(\varvec{x}_i)\le 0\right] = Z(f;\mathcal {D}_M)\nonumber \\&\quad +\sum _{i\in \mathcal {S}} \mathbbm {1}\left[ y_i\tilde{f}(\varvec{x}_i)\le 0\right] =Z(f;\mathcal {D}_M)+C. \end{aligned}$$

(19)

We now proceed to prove the statement in (11). When $\mathcal {S}=\emptyset$, then $\mathcal {D}_N = \mathcal {D}_M$, and (11) follows trivially. When, $\mathcal {S}\ne \emptyset$, we note that

$$\begin{aligned} \mathcal {F}^*={{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_N)&= {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_M \cup \mathcal {S}), \nonumber \\&= {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_M) + Z(f; \mathcal {S}), \nonumber \\&= {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_M) + C, \\&= {{\mathrm{argmin}}}_{f\in \mathcal {F}} Z(f;\mathcal {D}_M). \nonumber \end{aligned}$$

(20)

Here, the statement in (20) follows directly from (19). $\square$

1.4 Proof of Theorem 6 (Sufficient conditions to satisfy the level set condition)

Proof

We assume that we have found a surrogate function, $\psi$, that satisfies conditions I–IV and choose $C_{\psi } > 2\varepsilon$.

Our proof uses the following result: if $\Vert \varvec{\lambda }^{*}_{01}- \varvec{\lambda }^{*}_{\psi }\Vert > C_{\varvec{\lambda }}$ then $\varvec{\lambda }^{*}_{01}$ cannot be a minimizer of $Z_{01}\left( \varvec{\lambda }\right)$ because this would lead to a contradiction with the definition of $\varvec{\lambda }^{*}_{01}$. To see that this result holds, we use condition III with $\varvec{\lambda }=\varvec{\lambda }^{*}_{01}$ to see that $\Vert \varvec{\lambda }^{*}_{01}- \varvec{\lambda }^{*}_{\psi }\Vert >C_{\varvec{\lambda }}$ implies $Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) - Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) >C_{\psi }$. Thus,

$$\begin{aligned} Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) + C_{\psi }&< Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) \nonumber \\ Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) + C_{\psi }&< Z_{01}\left( \varvec{\lambda }^{*}_{01}\right) + \varepsilon \end{aligned}$$

(21)

$$\begin{aligned} Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) + C_{\psi } - \varepsilon&< Z_{01}\left( \varvec{\lambda }^{*}_{01}\right) \nonumber \\ Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) + C_{\psi } - \varepsilon&< Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) \end{aligned}$$

(22)

$$\begin{aligned} Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) + \varepsilon&< Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) . \end{aligned}$$

(23)

Here the inequality in (21) follows from condition IV, the inequality in (22) follows from condition I, and the inequality in (23) follows from our choice that $C_{\psi } > 2\varepsilon$.

We proceed by looking at the LHS and RHS of (23) separately. Using condition I on the LHS of (23) we get that:

$$\begin{aligned} Z_{01}\left( \varvec{\lambda }^{*}_{\psi }\right) + \varepsilon \le Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) + \varepsilon . \end{aligned}$$

(24)

Using condition IV on the RHS of (23) we get that:

$$\begin{aligned} Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) \le Z_{01}\left( \varvec{\lambda }^{*}_{01}\right) + \varepsilon . \end{aligned}$$

(25)

Combining the inequalities in (23), (24) and (25), we get that:

$$\begin{aligned} Z_{01}\left( \varvec{\lambda }^{*}_{\psi }\right) < Z_{01}\left( \varvec{\lambda }^{*}_{01}\right) . \end{aligned}$$

(26)

The statement in (26) is a contradiction of the definition of $\varvec{\lambda }^{*}_{01}$. Thus, we know that our assumption was incorrect and thus $\Vert \varvec{\lambda }^{*}_{01}- \varvec{\lambda }^{*}_{\psi }\Vert \le C_{\varvec{\lambda }}$. We plug this into the Lipschitz condition II as follows:

$$\begin{aligned} Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right) - Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right)&\le L \Vert \varvec{\lambda }^{*}_{01}- \varvec{\lambda }^{*}_{\psi }\Vert < L C_{\varvec{\lambda }}, \\ Z_{\psi }\left( \varvec{\lambda }^{*}_{01}\right)&< L C_{\varvec{\lambda }} + Z_{\psi }\left( \varvec{\lambda }^{*}_{\psi }\right) . \end{aligned}$$

Thus, we have satisfied the level set condition with $\varepsilon = LC_{\varvec{\lambda }}$. $\square$