Loss-Optimal Classification Trees: A Generalized Framework and the Logistic Case

The Classification Tree (CT) is one of the most common models in interpretable machine learning. Although such models are usually built with greedy strategies, in recent years, thanks to remarkable advances in Mixer-Integer Programming (MIP) solvers, several exact formulations of the learning problem have been developed. In this paper, we argue that some of the most relevant ones among these training models can be encapsulated within a general framework, whose instances are shaped by the specification of loss functions and regularizers. Next, we introduce a novel realization of this framework: specifically, we consider the logistic loss, handled in the MIP setting by a linear piece-wise approximation, and couple it with $\ell_1$-regularization terms. The resulting Optimal Logistic Tree model numerically proves to be able to induce trees with enhanced interpretability features and competitive generalization capabilities, compared to the state-of-the-art MIP-based approaches.


Introduction
Classification Trees (CTs) are a popular machine learning model for classification problems [1,2].Introduced in a seminal work by Breiman et al. in 1984 [3], CTs have widely been employed for decades, especially for tasks with small-sized, tabular data.
Formally, a CT is a connected acyclic graph whose nodes are linked in a hierarchical manner.Each internal node (branch) acts by splitting data in the feature space, based on some predefined condition.Finally, nodes that have no children (i.e., there is no outbound arc connecting them to lower nodes in the hierarchy) are called leaves and are associated with a prediction label for samples that are forwarded to them by the above branches.
Most algorithms designed to construct classification trees define univariate (axis-aligned) splits [3][4][5].With univariate splits, a single feature is chosen at each node and the corresponding value of each data point is compared to a given threshold; data points are thus forwarded to one of the children nodes, based on the result of this comparison.Each prediction of the model is therefore finally obtained following a path along the tree, based on the sequence of "decisions" at each encountered node.
In order to improve the expressive power of CTs, however, more complex splitting rules have been also considered in the literature for defining branching splits.These rules can take into account linear (in this case we talk about oblique trees) or even nonlinear relations between features [6][7][8].
Yet, the intrinsic interpretability of CTs is a key factor that makes them particularly popular in domains, such as healthcare, finance or justice, where understanding the decision-making process and accounting for the choices made is important [9].For this reason, recent research has dealt with training algorithms to construct CTs having both good predictive performance and intrinsic interpretability properties.In particular, oblique trees are often considered interpretable as long as branching linear classifiers are actually sparse models [10][11][12].Then, the rules that ultimately define the decision process are basically "if-then-else" conditions with very simple clauses; hence, domain experts with possibly no technical background are still typically able to interpret and understand the logic that leads to decisions.
Within this scenario, the contribution of the present paper consists of the encapsulation of the concept of classification loss within each splitting rule in the popular optimal CTs (OCTs) framework [13].This point of view also generalizes the idea of [14], where maximum margin splits were shown to improve the generalization performance of optimal trees.In fact, both OCTs and margin-optimal trees can be seen as particular instances of the generalized loss-optimal classification tree framework discussed in this manuscript.
Following further this path, we propose to choose the log loss to construct optimal logistic classification trees (OLCTs).This way, the nice properties of logistic regression models (e.g., generalization, calibration, interpretability) can be combined with the easily readable structure of CTs.We show that the logistic loss can be effectively encapsulated within the standard mixedinteger linear programming (MILP) models for OCTs, using the piece-wise linear approximation defined in [15], and thus handled by the usual off-the-shelf MILP solvers.
In addition, we point out a simple patch for the soft feature selection strategy employed in [14] for maximum margin CT models: the sparsity requirement can be more easily handled exploiting ℓ 1 -regularization; clearly, this idea can then be extended to the case of general loss-optimal models and thus to the OCLTs framework, which finally allows to obtain sparse, yet effective classification models with the additional interpretability properties of logistic regression.
The rest of the manuscript is organized as follows: in Section 2 we report a detailed review of the literature concerning classification trees induction.Then, in Section 3, we derive the general MIP framework for optimal classification trees, based on the concepts of loss and regularizer.We consequently introduce in Section 4 a novel specification of the general framework, leading to the definition of the Optimal Logistic Classification Tree model.In Section 5, we then report the results of computational experiments aimed at assessing the actual potential of the proposed approach.Finally, we give some concluding remarks in Section 6.

Related literature
The construction of an optimal decision tree is an N P-complete problem [16] and, for this reason, traditional algorithms apply top-down greedy strategies [17], employing quality measures to define the best parameters for each branch node [3][4][5].
Easily applicable even with large datasets, unfortunately greedy approaches often generate sub-optimal structures with poor generalization capabilities [18,19].This is mainly due to the fact that decision rules at the deepest nodes are often defined taking into account very small portions of the dataset.Once the tree has been grown, these methods usually apply a post-pruning phase to reduce the complexity of the structure and alleviate the overfitting problem [3,20,21].
In order to overcome these performance drawbacks, several methods have been proposed in the literature.One of the most common approaches consists of resorting to ensembles to reduce the model variance.Particularly relevant techniques of this kind are Random Forest [22], Gradient Boosting Machines [23] and XGBoost [24] which combine, in different ways, multiple trees to increase the overall performance.
Other works are focused on oblique classification trees, i.e., on the usage of hyperplanes to recursively divide the feature space.In this case, the tree model exploits a multivariate linear function at each branch node so that multiple features are involved in the decision process [25][26][27].Although both ensembles and oblique trees are more expressive than standard univariate classification trees, both classes of approaches suffer from an evident loss of interpretability.
A conceptually different path to improve the performance concerns the improvement of the fitting procedure, rather than the development of more expressive models.In particular, exact formulations of the learning problem exploiting linear programming have been considered.Back in the early '90s, Bennet et al. proposed linear methods for constructing separation hyperplanes [28,29] and for the induction of oblique classification trees solving an LP problem at each branch node [30].
More recently, thanks to the outstanding improvements in both hardware power and software solvers capabilities [31], new formulations of the learning problem have been developed.These formulations are related to both univariate and multivariate splits, and exploit various mathematical optimization techniques [32] for an effective modelling of the learning problem.One particularly prominent approach involves the use of Mixed Integer Linear Programming (MILP) formulations.The seminal work in this context is the one by Bertsimas and Dunn [13], where an exact MILP model is proposed that can be solved up to a certifiably globally optimal CT (OCT) structure (in terms of misclassification error).Inspired by this work, a plethora of improved MILP-based approaches followed; among them, we can notably cite a formulation of the problem for the case of binary classification with categorical features [33], or a flow-based formulation exploiting linear relaxation and Bender's decomposition to efficiently handle the specific case of binary classification with binary features [34].A reformulation for the OCT model was also proposed for the case of parallel splits, allowing to significantly improve the efficiency of the approach [35].Focusing on the interpretability aspect of CTs, some works take advantage of regularization terms to limit the number of branch nodes [13] and total number of leaves in the model [36,37].
The above formulations can in principle generate certified globally optimal structures.In practice, however, there is a combinatorial explosion in the number of binary variables as the depth of the tree or the size of the dataset increases.Thus, optimality gap can be closed by branch-and-bound type procedures only for really shallow trees and on tasks with a few hundred examples at most.
For this reason, local optimization strategies have also been investigated.Starting from a given tree, trained by a greedy approach, these methods refine its structure by iteratively minimizing the misclassification loss associated with each node in the tree.Carreira-Perpinán and Tavallali [38] introduce an alternating minimization strategy, that operates level-wise on the tree, decomposing the learning problem and solving each sub-problem up to global optimality.Similarly, Dunn proposed Local Search methods [39], able to handle real world datasets and to refine greedy structures grown with CART-like algorithms, at the cost of producing suboptimal models.
In addition to exploring integer optimization, researchers have also delved into continuous optimization approaches within the context of optimal trees.Blanquero et al. presented a nonlinear programming model in their work [40], aiming to develop an optimal "randomized" classification tree with oblique splits.This approach involves making random decisions at each node based on a soft rule, which is induced by a continuous cumulative density function.
Expanding on their earlier findings, Blanquero et al. addressed the issues of global and local sparsity in the randomized optimal tree model in [41].To tackle these challenges, they incorporated regularization terms based on polyhedral norms.By employing this regularization technique, the researchers aimed to promote sparse solutions both globally and locally within the randomized tree framework.The same approach has also been investigated for the regression case in [42].
It is important to note that, in the randomized framework proposed by Blanquero et al., the assignment of a sample to a specific class is not deterministic.Instead, it is determined based on a given probability, allowing for a more flexible and probabilistic classification approach.
A different stream of research aims at improving the performance of classification trees exploiting the concept of loss in defining splits.Specifically, there are works that deal with oblique splits, greedily computing at each branch node the parameters of the separating hyperplane as a Support Vector Machine (SVM) [43].The formulation from [44] for SVM branches exploits a dual convex quadratic model that can also handle kernel functions to capture non linear patterns in data.Tibshirani and Hastie instead developed a greedy algorithm to handle high dimensional features and multiclass classification, inducing classification trees with margin properties, i.e., structures where each branch node employs a linear SVM to split the feature space [45].
With reference again to multivariate trees, approaches that exploit logistic models at nodes have also been studied.In [46] a method is presented that induces standard axis-aligned classification trees with the difference that, at each leaf node, a logistic regression model is built to provide the final prediction.Moreover, a greedy method to fit a piecewise linear logistic regression model is presented in [47], exploiting decision trees to recursively partition the feature space.The final model can be basically seen as a multivariate classification tree with a logistic regressor on each branch node.
Recently, D'Onofrio et al. [14] took advantage of an exact mixed-integer quadratic programming (MIQP) formulation of the CT learning problem including the concept of maximum margin.This approach allows to define OCTs where each branch node is a maximum-margin linear classifier.The maximum margin property, obtained minimizing the hinge loss at all branch nodes, allows for significant improvements on generalization performance.However, relying on oblique splits has the usual drawback of undermining interpretability.For this reason, the authors also proposed strategies (both hard and soft) to induce sparsity in the coefficients of each hyperplane.
3 Generalized Loss-Optimal Classification Trees Let I = {1, ..., n} and D := {(x i , y i ), x i ∈ R p , y i ∈ {−1, 1}, i ∈ I} be a finite dataset of n observations with p features and binary labels.We introduce the notation, largely based on [13,14], used for the formulation of the Loss-Optimal CTs learning problem.Formally, let T be the set of branch nodes of a generic (oblique) CT.We denote by d the number of layers of 0 Fig. 1: An oblique tree of depth 3. Branches are on white background, leaves on gray.
branching nodes of the tree, i.e., its depth, by H the set {1, . . ., d} and by T h the set of nodes at depth h, for h ∈ H.Each branch node t ∈ T is characterized by a linear function w T t x+b t , w t ∈ R p , b t ∈ R, which induces the hyperplane w T t x + b t = 0 as decision boundary.This has the effect of forwarding an examined data point x i to the left child of node t if w T t x i + b t ≤ 0 and to the right child otherwise.Note that, using this notation, a glass-box axis-aligned CT can be obtained as a special case, imposing ∥w t ∥ 0 = 1, ∀ t ∈ T , where ∥ • ∥ 0 denotes the ℓ 0 pseudo-norm.
Since we are interested in CTs where each node is a binary classifier, we can let the final prediction of the point depend solely on the decision of the last branch node.Indeed, each splitting hyperplane is a linear classifier with decision function sgn(w T t x + b).In view of this, there is no need of tracking the points up until the leaves, so that the number of nodes to be modelled can be reduced of 2 d units, or, in other words, cut to a half.An example of an oblique CT with branching classifiers of depth 3 is shown in Figure 1.
We now introduce the essential building blocks of the MILP models proposed in the literature for OCTs training.Recalling that T d ⊂ T is the set of nodes of the last branching layer, the routing of each point to the proper branch of the last layer can be modelled introducing the binary variables: Using these variables, we can force each point to be routed to exactly one of the last branching nodes, complying with the structure of the model.In particular, this can be obtained imposing the following constraints: We then introduce the constraints to take into account the path of each data point, which is fully determined by the sequence of decisions taken at the branch nodes.For this purpose, we introduce the set T d (t) ⊆ T d as the set of branches at the last layer that are successors of the node t; in other words, nodes in T d (t) belong to the last layer of the sub-tree rooted at t; moreover, we define the subsets T r d (t), T l d (t) as a partition of T d (t) such that: is the set of nodes of the last branching level that belong to the right sub-tree rooted at node t; is the set of nodes of the last branching level that belong to the left sub-tree rooted at node t.Each branch t forwards points to its left child if w T t x + b t ≤ 0 and to the right one otherwise.This logic can be enforced by means of the following forwarding constraints: where M is a suitable constant for a big-M type constraint and ϵ is a small constant preventing numerically degenerate cases.For example, given t ∈ T and a data point x i , if z i,s = 1 and s ∈ T d (t), then x i has to be routed to the last branch node s, which either belongs to T l d (t) or T r d (t).Hence, one and only one of the following pair of equations hold: If the former holds, then (2a) guarantees that The big-M strategy makes the constraint for the "wrong case" irrelevant.Also, no constraint applies for branch t in case x i is forwarded to a node s ∈ T d \ T d (t).Moreover, we shall observe that these constraints are not required if t ∈ T d since the last branching level is only responsible to make the final prediction.

The Role of Loss Functions
For a CT to be meaningful, a suitable loss function necessarily has to be used to push the last branching layers and the overall model to correctly classify the data points.This can be done by estimating the errors, or slacks ξ i , i = 1, . . ., n, committed on each data point: The definition of quantities ξ clearly depends on the particular loss function to be used.For example, most MIP models for OCTs [13,34,36,38] directly employ the misclassification loss, which is defined using additional binary variables and big-M constraints as The additional variables ŷ model the predicted class for each data point; the slack variable corresponding to each data point will be set to 0 if prediction is correct y i ŷi = 1, otherwise it will be equal to 1. Constraints (4a)-(4b) guarantee that, for a point x i arriving at the leaf t, ŷi = 1 if and only if w T t x i + b t ≥ 0. However, this is not necessarily the only option for a loss function.In fact, different choices not only might be statistically more robust, but may also avoid the introduction of the additional binary variables and logical constraints that increase problem complexity.An example of a continuous loss with these features, easily embeddable in the MIP framework, is the hinge loss, max{0, 1 − y i (w T x i + b)}, that can be modeled as: The objective function can then be defined as the sum of the total loss function L and a regularization term for weights w t , t ∈ T d , so that we basically have an empirical risk minimization problem: We note that: • By setting λ = 0 and using (4), we actually retrieve the standard OCTs problem from [13]; 2 and using (5) we get "SVM leaves".
In fact, loss terms can also be associated with the upper branching nodes; in this case, we have d vectors of slack variables, ξ 1 , . . ., ξ d : along its path to the leaves, each data point encounters only one node at each layer, and thus only the slack associated with the corresponding classifier has to be taken into account.The overall objective function for a general loss-optimal CT model is therefore given by h∈H i∈I It is worth noting at this point that the Margin Optimal CT models (MARGOT) from [14] can then be retrieved as a special case of the general framework, setting and

Exact Modelling of ℓ 0 Terms
In order to enhance the interpretability of the models, sparsity can be compelled within the weights of branching classifiers, thus inducing features selection.The ℓ 0 norm of vectors w t can easily be modeled, in a MINLP program, by introducing binary indicator variables, big-M constraints and linear expressions: Then, the value of ∥w t ∥ 0 can either be upper bounded or penalized.Following the terminology in [14] we talk about Hard Feature Selection (HFS) in the former case and Soft Feature Selection (SFS) in the latter one.

Optimal Logistic Classification Trees
In this section, we formalize a novel, particular instance of loss-optimal CT model, the Optimal Logistic CT (OLCT).In particular, we first show how to introduce the logistic loss function within the considered MIP; then, we point out the benefits of using ℓ 1 -regularization terms; we then show the resulting overall optimization model.

Logistic Loss in Mixed Integer Linear Optimization
The logistic loss function for binary linear classifiers is defined as This loss function appears in the individual terms of the summation in the negative log likelihood function of logistic regression models (see, e.g, [48]), i.e., the linear model for binary classification that is obtained by maximum likelihood estimation under the assumption that data follows a Bernoulli distribution.
In addition to often having strong performance in terms of out-of-sample prediction accuracy, other advantages of logistic regression compared to other linear classifiers, such as SVMs, include (a) the possibility of obtaining estimates of features importance by simple manipulation of the model weights [48]; (b) the opportunity of getting calibrated probability estimates associated with predictions.The above properties offer nice insights that are valuable from the perspective of model interpretability.These considerations motivate us to consider the employment of the logistic loss within the general framework for loss-optimal CTs, discussed in Section 3.
The straightforward objection, at this point, concerns the nonlinearity of the logistic loss function; indeed, contrarily to the hinge loss defining SVMs, the log loss cannot exactly be modeled by linear constraints in MILP.However, this issue can be addressed following the strategy, proposed in [15], where the best subset selection problem in logistic regression is solved by means of a MILP approach.In particular, Sato et al. proposed to approximate the logistic loss function by a piece-wise linear underestimator.The function is a convex function, thus its tangent line at a point v 0 constitutes a global underestimator; more explicitly, for all v, v 0 ∈ R, it holds To obtain an accurate approximation of the logistic loss, we can thus construct a piece-wise linear underestimator obtained as the point-wise maximum of a family of tangent lines: Sato et al. also propose a greedy strategy to select points v k where computing the tangent lines so as to minimize the approximation error: at each iteration, a tangent line is added to the piece-wise linear approximation so that the area between the exact and approximated loss function is minimized.The resulting sets of tangent points to be used for increasingly accurate approximations of the logistic loss are Obviously, a larger number of points leads to a larger number of constraints in the MILP model and thus makes it more difficult to solve.Following this methodology, we can redefine the slack variables to finally introduce the logistic loss in the loss-optimal CT model discussed in Section 3:

Lasso Regularization
Together with the loss function defining the slack values, the second element to be chosen in our generalized OCT framework is the regularizer.Using an ℓ 2regularization term, as for example in MARGOT, has the effect of making the entire problem an MIQP instance.Here, we instead propose to consider a Lasso regularizer [49], i.e., an ℓ 1 penalty term; there are two main reasons for doing so: • the ℓ 1 -norm can be easily handled by linear constraints within a linear programming model; thus, using the ℓ 1 -norm instead of the squared ℓ 2norm we can derive a fully linear model which should be easier to solve; • exploiting the well-known properties of the ℓ 1 -norm [50], we can implicitly induce sparsity within branch nodes classifiers, without the need of recurring to explicit (and expensive) ℓ 0 -norm penalization as in SFS models.The ℓ 1 -norm can be efficiently handled in a linear program, similarly as what is done in [51], setting Basically, w t is split into its positive and negative parts w + t and w − t .Indeed, constraints (11) are satisfied by ŵ+ j,t = max{0, w j,t } and ŵ− j,t = max{0, −w j,t }; this solution is such that ŵ+ j,t + ŵ− j,t = |w j,t |.We thus have 1 Any other feasible solution can be obtained by shifting ŵ+ t and ŵ− t by a vector ∆ ≥ 0, ∥∆∥ > 0; hence, for any other feasible solution we have The shift ∆, however, has no influence on the variables w t ; if the term 1 T (w + + w − ) is minimized in the objective, then ŵ+ t and ŵ− t will always be chosen and thus the actual quantity being minimized is Lasso regularization is well known not only to induce sparsity and guarantee that the optimization process is well-behaved, but also to be beneficial at tackling overfitting [48,49].Thus, if the sparsity requirement within each node is not set by a specific budget (as in HFS), but it is imprecisely imposed by means of a penalty, then the expedient of setting Ω(•) = ∥ • ∥ 1 allows to preserve, at a much lower cost, both the predictive performance of the obtained model and its interpretability.
The alternative SFS strategy from [14] exploits the objective function t∈T This kind of penalty on the ℓ 0 -norm of weights is modeled with binary variables, see eq. ( 9), highly increasing the complexity of the optimization model.Moreover, here there is an additional hyperparameter to be tuned for each layer, making the cross-validation procedure harder to manage.

The Overall Model
The overall Optimal Logistic Classification Tree (OLCT) model proposed in this paper is obtained putting together the pieces described thus far: s.t. ( 1), ( 2), ( 10), (11) (12b) The constraints referenced at (12b) contain almost all the defining elements of the logistic tree: constraint (1) defines the indicator binary variables of data points routing, constraints (2) enforce the routing logic along the tree, constraints (10) define the slacks corresponding to the logistic loss and (11) define the branching classifiers weights and model their ℓ 1 -norm.
Constraints (12c) and (12d) simply define the domain of the remaining variables.The objective function is nothing but the general loss (7) where the ℓ 1 regularizer is employed, formulated as shown in Section 4.2.
Remark 1.Note that, in our model, we use a piece-wise linear underestimator, i.e., a surrogate function, to approximate the log-loss.Thus, at the end of the training process, we can actually perform a refinement operation affecting the last layer of branching nodes.Specifically, we can exactly fit an ℓ 1 -regularized logistic regression model at each node of the last layer, using as training data the samples that actually reach that node.By this procedure, we are guaranteed that the overall exact objective function associated with the entire tree model decreases.
The same reasoning cannot be applied with higher level nodes: changes to the branching classifiers possibly change how training data are forwarded to lower nodes, thus affecting the corresponding loss terms.The overall loss associated with the tree might increase, because of the lack of global perspective.
Remark 2. One of the most appealing features of classification trees is their glass-box nature; however, the highest level of interpretability is only reached in the case of parallel splits.For this reason, we believe it to be important to explicitly point out how to retrieve a univariate optimal logistic classification trees.The model is basically equivalent to (12), except for HFS type constraints, with an upper bound on the ℓ 0 -norm of each branching classifier set to 1.This can be modeled in MILP terms by constraints (9) and setting 1 T δ t ≤ 1.
Of course, the addition of a number of new binary variables and big-M type constraints directly proportional to p × |T | is significant in terms of the computational resources needed to solve the problem and especially to certify optimality, closing the optimality gap.

Interpreting OLCTs
Inheriting, at least partially, the nice interpretability properties of logistic regression models is one of the main advantages of using the logistic loss within the OCTs framework.
Specifically, there are some aspects that can be taken into account to retrieve additional information about the model prediction mechanisms.
Evaluation of feature influence At each branching node of oblique OLCTs, the influence of each individual parameter in the splitting decision can be estimated looking at the magnitude of the corresponding coefficient and multiplying it by the standard deviation of the feature among the training data points reaching that node; formally, given the weights w t at a node t ∈ T , the influence r j, t of feature j can be estimated by: The larger the coefficient r j,t is, the higher is the connection of feature j with positive outputs; on the other hand, the largest negative values of importance are associated with features pushing the point towards the left child of the branching node.The features influence for an OLCT trained on the heart dataset is reported, as an example, in Figure 2. Fig. 2: A learned logistic oblique tree of depth 2 for the heart dataset.Branches and leaves are on white and gray background respectively.Each branch node is a sparse linear regressor; the importance coefficient for each feature involved at each decision is reported in the boxes.

Probabilistic interpretation of outputs and model calibration
The sigmoid function σ(z) = 1 1+exp(−z) is used in logistic regression to map the output of the linear function defined by w, b to (0, 1), so that a data point x can finally be classified as positive if σ(w T x + b) ≥ 0.5.
In a standalone logistic regression model, these "probability values" are related to the odds of the positive outcome over the negative one; in particular, the linear regressor is designed to model, by maximum likelihood, the log-odds (logits) of the output: .
Exponentiating we get and by simple algebraic manipulations we retrieve In other words, the logistic output is an actual probability estimate for the positive output, i.e., the logistic model is well calibrated by construction.Within the OLCTs framework, we can exploit this property both at the splitting and the classifying nodes of the tree.If t ∈ T d , then the odds associated with the outputs of w t , b t are actually interpretable as the odds of the positive class over the negative one, given that the data point belongs to the subspace deterministically defined by the splits at the higher nodes.Since each node classifier is calibrated within the corresponding space region, the overall output probability of the OLCT models should also be implicitly well-calibrated, with no need of any extra post-train calibration.
If, on the other hand, t ∈ T h , h < d, then the classifier has been trained looking not only at its corresponding slacks, but also at those of all its descendant nodes.Thus, in this case we can (somewhat improperly) interpret the probability estimates as the confidence about forwarding the data point to the right child rather than to the left one.This concepts can be visualized through the example in Figure 3. of the heart dataset.The true label y is equal to 1.In this case, the model predicts the correct class of the point and, given the sigmoid activation, we are able to get the confidence of the forwarding decision at each branch node.

Numerical Experiments
In this section, we present the results of computational experiments carried out to evaluate the performance of the proposed approach.All the experiments described in this section have been carried out on a server with an Intel ®Xeon ®Gold 6330N CPU with 28 cores and 56 threads @ 2.20GHz, but we set a limit of only 40 of the 56 available threads and the total available memory is 128GB.The code has been implemented in Python (v.3.9) and the commercial solver Gurobi [52] has been used to solve all the mixed-integer programming models considered in this work.All the code is available at https://github.com/tom1092/Optimal-Logistic-Classification-Trees.
For each instance of classification tasks, we performed an 80/20 train/test split of the data and we also standardized each feature before the training to zero mean and unit variance.Experiments are always repeated for different random seeds, resulting in different train/test splits.Hyperparameters have been tuned by cross-validation over a grid of values, where the test balanced accuracy (see below) is used as quality metric; more details will be provided for each group of experiments in the following.The parameters M and ϵ in MIP formulations have been set to 100 and 10 −5 respectively; the value of M = 100 is a reasonable value, large enough not to introduce bound constraints and tight enough not to make Gurobi branch-and-bound too inefficient; this same value was also used, for instance, in [14]; as for ϵ, the value 10 −5 is much larger than the one employed for the IntFeasTol parameter of Gurobi (10 −9 ), so that the issues highlighted in [53] did never occur, and at the same time is small enough to let the underlying logic of the constraints work properly -no feasible solution is erroneously cut off.Without loss of generality, we also decided to move the regularization parameter λ h to the slack component of the loss, to be aligned with the work in [14].The objective function now has the form: Note that, at the end of the optimization process of OLCTs, we applied the refinement strategy discussed in Remark 1.For each node in the last layer, we retrained the ℓ 1 -regularized logistic regression model using scikit-learn [54] implementation.
To compare the performance of each model, along with the running times, we used the balanced accuracy metric defined as: where TP, TN, FP, FN are the number of true positive, true negative, false positive and false negative outputs, respectively.The B Acc value is always computed on test set data.
To provide a condensed view of the results, in the following we are making use of performance profiles [55].Performance profiles provide a unified view of the relative performance of the solvers on a suite of test problems.Formally, consider a benchmark of P problem instances and a set of solvers S. For each solver σ ∈ S and problem π ∈ P, we define c π,σ = the cost for solver σ to solve problem π, where cost is the performance metric we are interested in.In particular, we will be interested in CPU time.We then consider the ratio which expresses a relative measure of the performance on problem π of solver σ against the performance of the best solver for this problem.If a solver fails to solve a problem, we shall put η π,σ = η M , with η M ≥ max{η π,σ | π ∈ P, σ ∈ S}.
Finally, the performance profile for a solver σ is given by the function which represents the estimated probability for solver σ that the performance ratio η π,σ on an arbitrary instance π is at most τ ∈ R. The function ρ σ (τ ) : [1, +∞] → [0, 1] is, in fact, the cumulative distribution of the performance ratio.Note that the value of ρ σ (1) is the fraction of problems where solver σ attained the best performance; on the other hand, lim τ →η − M ρ σ (τ ) denotes the fraction of problems solved from the given benchmark.

Preliminary experiments
The first experiments we carried out concern the assessment of the performance of our model as we vary the set V of the tangent points that are used to obtain the piece-wise linear underestimator of the logistic loss.In particular, we are interested in finding a good trade-off between the quality of the approximation and the running time of the model.
We considered a small benchmark of 5 datasets (parkinsons, wholesale, tik-tak-toe, haberman, sonar, see Table 1) from the UCI repositories [56], testing V 0 , V 1 , V 2 with refinment as configurations for the MILP problem of training an OLCT of depth 2. The experiment has been repeated for five different random seeds for each dataset for a total of 25 different problems.The Gurobi time limit has been set to 300s.In these experiments, we did not employ the ℓ 1 regularization terms, as we are mainly interested in assessing the effectiveness of log loss approximation.Note that the MIP solution process is warm-started, initializing the weights of each branch node following a strategy similar to the one adopted in [14]: starting from the root in a greedy fashion, we assign to each node weights the values obtained training a logistic regression Loss-Optimal Classification Trees classifier using the data that are forwarded to that node by the above branches.This strategy allows to obtain significant speedups in computation.
As shown by the performance profiles in Figure 4, choosing V 0 = {0, ±∞} to build the linear piece-wise approximation seems to provide a nice trade-off between running time and solution quality.Indeed, the use of V 1 does not seem to significantly improve the out-of-sample accuracy of the model; on the other hand, the more accurate approximation obtained with V 2 does result in better predictive models, but a much higher training cost has to be paid.For this reason, we decided to use the V 0 setting to assess the performance of our model in the following sections.Nonetheless, we do not rule out that, in certain settings, it might be worth exploiting the increased effectiveness provided by V 2 .Fig. 4: Comparison of the performance of OLCT models when different tangent point sets (V 0 , V 1 or V 2 ) are employed.At the end of the optimization process, the refinement of the last layer was carried out for each model.

Impact Assessment for the Global Optimization Approach
In this Section, we investigate the actual beneficial effects of conducting a global optimization phase during logistic trees training.Indeed, good performance of OLCTs might be mainly due, in principle, to the warm-start or the refinement steps.In particular, our greedy warm-start procedure constructs a logistic tree model with an iterative top-down approach, which is somewhat similar to the one proposed in [47] and therefore might actually already lead to a properly effective classifier.However, the results reported in Figure 5 suggest that solving model ( 12) does provide a substantial boost to the performance of the resulting logistic tree.
In this experiment, we examined the value of the overall (in-sample) logistic loss associated with the model after the warm-start, global optimization and refinement steps, and then we also took into account the (out-of-sample) values of balanced accuracy.Here, we solved the same 25 instances of classification problems considered in Section 5.1.Again, in order to focus on the consequences of approximating the loss in the global step, we did not employ the ℓ 1 regularization terms.We also set Gurobi time limit to 300s.Moreover, based on the results in the previous section, we chose to employ the set V 0 of tangent points in log loss approximation.We did not focus on the running times of the three phases, as the global optimization step obviously represents the main computational burden for our approach.
From Figure 5(b), we observe that solving the training problem with a global structure perspective generally allows to slightly improve the overall loss attained by the greedy model.This result has to be underlined, taking into account that the loss function is roughly approximated during the global optimization phase.Then, the refinement step on the last layer allows to really polish the model on training data, often leading to substantially lower values of loss.
The results in Figure 5(a) are even more appealing.Apparently, when it comes to test performance, learning branching rules with a global point of view is crucial to improve the effectiveness of the resulting model.In this perspective, although visible, the positive effect of refinement is much more limited.Thus, we can arguably state that the solving (12) as a global optimization has a significant effect in improving the effectiveness of logisitic tree models.

OLCTs Performance Evaluation
We now present the results of a larger computational experiment where we compare the OLCT model to the MARGOT, SFS-MARGOT [14] and OCT-H [13] approaches.To assess the performance of our method, we considered 10 standard binary classification datasets from the UCI repositories [56] that are reported in Table 1.For each dataset, we consider 5 classification problem instances obtained considering different random train/test splits, so that the overall benchmark is made up of 50 test instances.For each problem instance,   we set the depth of each tree model to 2. For all models, we set the Gurobi time limit to 300s both for validation runs and the final fit over the entire training set.Table 1: Description of the datasets used in the computational experiments.All datasets are from the UCI collection [56].
For hyperparameters tuning, this time we used a 4-fold cross-validation.For both OLCT and MARGOT, we considered the slack parameters C h ∈ {10 i , i = −2, −1, ..., 2} ∀ h ∈ H, obtaining 25 possible model configurations.On the other hand, the SFS-MARGOT model has two hyperparameters for the classifiers at each level: the slack parameter C h and the ℓ 0 penalty parameter α h .In order to consider a comparable grid of configurations in size as that of OLCTs and MARGOTs, we used the same α for each branch layer h letting C h vary in {10 −2 , 1, 10 2 } and α in {10 −1 , 1, 10}, so that a total number of 27 models is considered for the SFS variant.Finally for OCT-H we used the grid α ∈ {2 i , i = −8, −1, ..., 2} ∪ {0} to tune the α parameter that penalizes the number of features used at each branch node to make the decision.
Again, we initialize the training phase of each MIP model by injecting a warm start solution, obtained training logistic or SVM classifiers, depending on the particular tree classifier.For OCT models we used an analogous greedy strategy: we solve the MILP problem at each individual single node, i.e., setting the depth of the tree equal to 1 and only using the set of points reaching the considered branch node, with a time limit of 30s.As also mentioned in [13], this strategy can significantly speedup the optimization, providing a good initial upper bound of the loss that may help the branch and bound method.The results of the experiment are shown in Figure 6 in the form of cumulative distribution of absolute gap from the best balanced accuracy on the test set and performance profiles [55] of runtime and sparsity.Sparsity is measured by the average number of features used at each node of the tree, and constitutes a proxy measure for interpretability of oblique trees.
Observing Figure 6(a) we can infer, for each model, estimates of the probability to obtain a test balanced accuracy value distant at most t points from the best one attained by any model.In other words, for any t, we can observe the fraction of times a model reaches an accuracy level within t points from the best one.This kind of graph is very informative.For example, for t = 0 we obtain the fraction of times each model is the best one among all those considered.From this point of view, we observe that our proposed model attains the top accuracy in about 50% of cases; moreover, it also appears to be the most robust, consistently being the most likely to obtain an accuracy value close to the best one, as the gap parameter t increases.Our model is thus not only most frequently the best one, but when it is not, it is still the one with the lowest probability of falling shorter than any given threshold from the best result.
Another interesting observation is that both SFS-MARGOT and MARGOT exhibit similar performances.That is, the SFS variant is able to produce sparse trees without drastically reduce the out-of-sample prediction performance.In this regard, from Figure 6(c) we can observe that OLCT outperforms both MARGOT and SFS-MARGOT in terms of sparsity, i.e., interpretability.
The high average levels of sparsity in OLCTs branches support the effectiveness of the ℓ 1 -regularization approach to this aim.This is of course not particularly surprising, as the effects of LASSO regularization are well known.Yet, the ℓ 1 -regularization approach appears to also be able to somewhat outperform the SFS strategy based on ℓ 0 penalization: this result was not granted and is certainly worth to be remarked, even more so taking into account that it is coupled with the efficiency advantage due to the avoided use of binary variables.
On the other hand, the OCT-H approach proves to be the best one in terms of sparsity.This result is driven by two elements within the model as described in [13]: the penalty term in the objective aimed to encourage splits considering a low number of features, and a constraint in the formulation which forces the weights of each split to have a unitary ℓ 1 norm.However, as a consequence of the combined use of these strategies, the final model tends to be "overregularized", resulting the clearly worst one in terms of balanced accuracy as shown in Figure 6(a).
Finally, in Figure 6(b), performance profiles of the running times highlight that the vanilla MARGOT is, in general, the most likely model to close the optimality gap in less than 300s.However, OLCT appears to have a comparable cost, whereas the SFS-MARGOT and the OCT-H approaches are much more computationally demanding since these latter models make use of many more binary variables.By the way, it is worth to notice that no model was able to close the optimality gap in more than 70% of the instances with a time limit of 300s.
Summarizing, results highlight that our method is able to outperform other approaches in terms of balanced accuracy, to increase the interpretability, inducing sparser structures and exploiting the well known logistic properties discussed in Section 4.4, and finally that these improvements can be achieved in competitive running times with respect to the other approaches.

Performance Analysis with Larger Scale Problems
The mixed-integer optimization models associated with loss-optimal classification trees grow very fast in size as the number of data points or nodes layers increases.In particular, the number of integer variables (and of constraints) grows linearly with the number of data points and exponentially with the trees depth.
This increment is problematic from the perspective of solving the training optimization problem, as even the most efficient solvers from the state-of-the-art struggle when the number of integer variables becomes large.
We are thus interested in conducting, at least, a preliminary assessment of the scalability of the OLCT approach compared to the behavior of MARGOT and OCT-H models.In this larger scale setting, we extended the time limit to 3600 seconds for Gurobi during the final fitting of the models, after the cross-validation phase.
First, we considered 9 additional datasets1 , whose size is reported in Table 2. Since the computational burden to carry out the present experiment is significant, we only considered a single train-test split for each dataset.The comparison concerns the OLCT model with the V 0 tangent points set, the MARGOT model and the OCT-H model, setting trees depth to 2. The results of the experiment are reported in Table 2.We do not report Gurobi running times as for no instance optimality of the solution was certified.In fact, for all models, the optimality gap was consistently very close to 100% when the time limit of 3600s was reached.
We can observe that OCT-H heavily struggled on these problems, producing meaningful models only with the least large instances of the benchmark.In the most difficult problems, neither the warm start and the global steps were able to find, within their time limits, anything better than a trivial tree forwarding all the points up to the same leaf, leading to a value of B Acc of 0.5.On the other hand, even if unable to certify optimality, both OLCT and MARGOT were able to handle the datasets.The out-of-sample performance of the two approaches is close, with OLCT apparently having a slight advantage; yet, the size of the benchmark does not allow us to state that one model is better than the other.
We then proceeded to evaluate the models as the depth of the trees increased.Specifically, we conducted experiments for trees with a depth of 3. It is important to note that, in this scenario, the cost of the experiments dramatically rises.The addition of an extra layer of nodes necessitates tuning an additional hyperparameter.Consequently, the number of hyperparameters configurations to be considered in the cross-validation phase increases accordingly, rendering the entire process time-consuming.As a result, we only considered 5 datasets from Table 1 with a single train/test split.The results of the experiment are reported in Table 3.We can observe that the loss-optimal classification tree framework remains manageable when used to train more complex tree classifiers with a depth of 3. In fact, the produced models demonstrate generalization capabilities, underlining the effectiveness of the procedure.OLCT models appear to have a slight advantage over MARGOT and OCT-H, both in terms of accuracy and ease of training.However, it is important to note that the size of the benchmark considered is too small to draw definitive conclusions in this regard.

Conclusions
In this work we proposed a general framework for loss-optimal classification trees, showing how different losses can be handled through a proper definition of the slack variables.We reviewed the recent state-of-art MIP methods for

3. 1
Notation and Oblique Tree Structure by Mixed-Integer Programming Fig. 3: Confidence at each branch node of the logistic classification tree from Figure 2 for the sample x = [0.69,0.87, −1.23, −0.8, −0.4,1.01, −1.87, 1.4, −0.94, 0.63, 0.35, −0.89, −0.07]of the heart dataset.The true label y is equal to 1.In this case, the model predicts the correct class of the point and, given the sigmoid activation, we are able to get the confidence of the forwarding decision at each branch node.
Cumulative distribution of the absolute gap from the best test balanced accuracy value attained by any model.
Performance profiles of the running times for solving the MIP models.
Cumulative distribution of the absolute gap from the best (exact) loss attained on the training set by any model.
Cumulative distribution of the absolute gap from the best test balanced accuracy value attained by any model.
Cumulative distribution of the absolute gap from the best (exact) loss attained on the training set by any model.
Cumulative distribution of the absolute gap from the best test balanced accuracy value attained by any model.
Performance profiles of the running times for solving the MIP models.
Performance profiles of the average number of features considered at each split of the classification tree.

Fig. 6 :
Fig. 6: Comparison of the performance of OLCT, MARGOT and MARGOT-SFS models, on a benchmark of 50 problem instances (10 datasets, 5 random seeds for train/test split).

Table 2 :
Performance of depth-2 trees fit based on different OCT models on larger scale datasets (available at https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/).Balanced accuracy (B Acc ) is reported for an 80/20 train/test split.We set 300s and 3600s as time limits for the validation and the training phase respectively.

Table 3 :
Performance of depth-3 trees fit based on different OCT models.Balanced accuracy (B Acc ) on test set and running times are reported on an 80/20 train/test split.We set 300s and 3600s as time limits for the validation and the training phase respectively.