New Advances In Data Science Problems Through Hyperplanes Location
Abstract
This thesis dissertation focus on developing new approaches for different Data Science problems from a Location Theory perspective. In particular, we concentrate on locating hyperplanes by means of solving Mixed Integer Linear and Non Linear Problems. Chapter 1 introduces the baseline techniques involved in this work, which encompass Support Vector Machines, Decision Trees and Fitting Hyperplanes Theory. In Chapter 2 we study the problem of locating a set of hyperplanes for multiclass classification problems, extending the binary Support Vector Machines paradigm. We present four Mathematical Programming formulations which allow us to vary the error measures involved in the problems as well as the norms used to measure distances. We report an extensive battery of computational experiment over real and synthetic datasets which reveal the powerfulness of our approach. Moreover, we prove that the kernel trick can be applicable in our method. Chapter 3 also focus on locating a set of hyperplanes, in this case, aiming to minimize an objective function of the closest distances from a set of points. The problem is treated in a general framework in which norm-based distances between points and hyperplanes are aggregated by means of ordered median functions. We present a compact formulation and also a set partitioning one. A column generation procedure is developed in order to solve the set partitioning problem. We report the results of an extensive computational experience, as well as theoretical results over the scalability issues and geometrical analysis of the optimal solutions. Chapter 4 addresses the problem of finding a separating hyperplane for binary classification problems in which label noise is considered to occur over the training sample. We derive three methodologies, two of them based on clustering techniques, which incorporate the ability of relabeling observations, i.e., treating them as if they belong to their contrary class, during the training process. We report computational experiments that show how our methodologies obtain higher accuracies when training samples contain label noise. Chapters 5 and 6 consider the problem of locating a set of hyperplanes, following the Support Vector Machines classification principles, in the context of Classification Trees. The methodologies developed in both chapters inherit properties from Chapter 4, which play an important role in the problems formulations. On the one hand, Chapter 5 focuses on binary classification problems where label noise can occur in training samples. On the other hand, Chapter 6 focus on solving the multiclass classification problem. Both chapters present the results of our computational experiments which show how the methodologies derived outperform other Classification Trees methodologies. Finally, Chapter 7 presents the conclusions of this thesis.
Full text
Programa de doctorado “Matem´aticas” PhD Dissertation NEW ADVANCES IN DATA SCIENCE PROBLEMS THROUGH HYPERPLANES LOCATION Author Alberto Jap´on S´aez Supervisors Prof. Dr. Justo Puerto Albandoz Prof. Dr. V´ıctor Blanco Izquierdo March 1, 2022
Resumen El trabajo de esta tesis se centra en desarrollar nuevas metodolog´ıas para abordar problemas cl´asicos de Ciencia de Datos desde el punto de vista de la Teor´ıa de la Localizaci´on. En particular, nos centramos en problemas de localizaci´on de hiperplanos que obtenemos modelando y resolviendo problemas de programaci´on lineal, y no lineal, entera mixta. El Cap´ıtulo 1 presenta las t´ecnicas desarrolladas en la literatura desde las que iniciamos este trabajo, que comprenden M´aquinas de Vector Soporte, ´ Arboles de Clasificaci´on y Teor´ıa de Ajuste de Hiperplanos. El Cap´ıtulo 2, basado en el art´ıculo Blanco et al. (2020b), se centra en el estudio de localizar un conjunto de hiperplanos para resolver un problema de clasificaci´on multiclase, extendiendo las M´aquinas de Vector Soporte al escenario multiclase. Presentamos cuatro formulaciones resultado de combinar dos formas de medir los errores de clasificaci´on y las normas utilizadas para medir las distancias. Reportamos los resultados de experimentos en conjuntos de datos reales y sint´eticos donde se muestra la alta capacidad predictiva de nuestros modelos frente a otros existentes en la literatura. Adem´as, demostramos tambi´en que el truco del kernel es aplicable en nuestro m´etodo. En el Cap´ıtulo 3, inspirado en el art´ıculo Blanco et al. (2021a), tambi´en tratamos el problema de localizar un conjunto de hiperplanos, sin embargo, en este caso lo hacemos bajo la premisa de minimizar una funci´on objetivo basada en el ajuste de distancias del conjunto de puntos a los hiperplanos. El problema lo tratamos en un contexto general donde las distancias de los puntos a los hiperplanos se contabilizan mediante funciones de tipo mediana ordenada (Ordered Weighted Averaging-OWA). Para la resoluci´on del problema reportamos una formulaci´on compacta y otra basada en conjuntos de las observaciones, para la cual desarrollamos un m´etodo de generaci´on de columnas. Reportamos numerosos experimentos que permiten evaluar las dos metodolog´ıas presentadas. Discutimos resultados te´oricos sobre la escalabilidad de los problemas y reportamos tambi´en algunos an´alisis geom´etricos de las soluciones. II
III El Cap´ıtulo 4, que sigue el trabajo planteado en Blanco et al. (2020a), trata el problema de encontrar un hiperplano separador de margen m´aximo con la particularidad de que consideramos que puede haber cierto ruido en las etiquetas de la muestra de entrenamiento. Desarrollamos tres metodolog´ıas, dos de ellas basadas en t´ecnicas de cluster, que incorporan en el entrenamiento la capacidad de reetiquetar observaciones, es decir, de considerar observaciones como pertenecientes a su clase contraria. En los resultados computacionales mostramos como nuestros modelos obtienen mayores tasas de acierto en la clasificaci´on cuando hay ruido aleatorio en las etiquetas de las muestras de entrenamiento. Los cap´ıtulos 5 y 6, basados en los trabajos presentados en Blanco et al. (2021c) y Blanco et al. (2021b) respectivamente, se centran en el estudio de ´ Arboles de Clasificaci´on en la que las ramas se definen mediante hiperplanos basados en M´aquinas de V´ector Soporte. Las metodolo´ıas desarrolladas en estos cap´ıtulos heredan propiedades de los m´etodos presentados en el Cap´ıtulo 4, que forman un pilar esencial en los mismos. Por un lado, en el Cap´ıtulo 5 se centra en el estudio de problemas de clasificaci´on binaria donde adem´as se contempla la posibilidad de existencia de ruido aleatorio en las etiquetas de las observaciones en las muestras de entrenamiento. Por otro lado, el Cap´ıtulo 6 se centra en el estudio de problemas de clasificaci´on multiclase. En ambos cap´ıtulos se presentan experimentos computacionales que muestran mejoras en las predicciones sobre bases de datos reales frente a otros m´etodos de la literatura. Finalmente, en el Cap´ıtulo 7 se muestran las conclusiones de la tesis.
Abstract This thesis dissertation focus on developing new approaches for different Data Science problems from a Location Theory perspective. In particular, we concentrate on locating hyperplanes by means of solving Mixed Integer Linear and Non Linear Problems. Chapter 1 introduces the baseline techniques involved in this work, which encompass Support Vector Machines, Decision Trees and Fitting Hyperplanes Theory. In Chapter 2, which is based in Blanco et al. (2020b), we study the problem of locating a set of hyperplanes for multiclass classification problems, extending the binary Support Vector Machines paradigm. We present four Mathematical Programming formulations which allow us to vary the error measures involved in the problems as well as the norms used to measure distances. We report an extensive battery of computational experiment over real and synthetic datasets which reveal the powerfulness of our approach. Moreover, we prove that the kernel trick can be applicable in our method. Chapter 3, which is inspired in Blanco et al. (2021a), also focus on locating a set of hyperplanes, in this case, aiming to minimize an objective function of the closest distances from a set of points. The problem is treated in a general framework in which norm-based distances between points and hyperplanes are aggregated by means of ordered median functions. We present a compact formulation and also a set partitioning one. A column generation procedure is developed in order to solve the set partitioning problem. We report the results of an extensive computational experience, as well as theoretical results over the scalability issues and geometrical analysis of the optimal solutions. Chapter 4, which follows the work presented in Blanco et al. (2020a), addresses the problem of finding a separating hyperplane for binary classification problems in which label noise is considered to occur over the training sample. We derive three methodologies, two of them based on clustering techniques, which incorporate the ability of relabeling observations, i.e., treating them as if they belong to their contrary class, during the training process. We report computational experiments that show how our methodologies obtain higher accuracies when training samples VI
VII contain label noise. Chapters 5 and 6, which are based in the works Blanco et al. (2021c) and Blanco et al. (2021b), respectively, consider the problem of locating a set of hyperplanes, following the Support Vector Machines classification principles, in the context of Classification Trees. The methodologies developed in both chapters inherit properties from Chapter 4, which play an important role in the problems formulations. On the one hand, Chapter 5 focuses on binary classification problems where label noise can occur in training samples. On the other hand, Chapter 6 focus on solving the multiclass classification problem. Both chapters present the results of our computational experiments which show how the methodologies derived outperform other Classification Trees methodologies. Finally Chapter 7 presents the conclusions of this thesis.
4 Chapter 1. Introduction Historically, Mathematics has been considered an abstract field which has no, or at best little, relation to the sensitive world for a majority of society. This is no longer true. The collective conscience about the generation and use of data has implied a change in the perception of mathematical knowledge amongst people. Nowadays, words like algorithms, Machine Learning, Artificial Inteligence or Big Data are not unlikely to appear in a casual conversation. At the same time, proper data management has proven to be an essential requirement for decision making, and according to this, institutions spend a large amount of resources on such a purpose. Statistics and Operations Research provide an extensive pool of techniques so as to deal with data. When facing a data problem in real life, two scenarios are mainly presented: On the one hand there are cases where, guaranteeing to provide a good solution, computational time is desired to be as small as possible. On the other hand, there are problems in which the quality of the solution is the most important part, therefore the optimal solution must be found despite the fact that it can be computationally expensive. This optimal solution is usually obtained by means of solving an optimization problem. In this thesis we focus on deriving and solving optimization problems in which a set of hyperplanes are located in order to solve some Data Science problems. Although Data Science does not have a concrete definition, authors usually rely on the same ideas. Before Data Science was even named, Tukey stated (Tukey, 1962): By and large, the great innovations in statistics have not had correspondingly great effects upon data analysis. The extensive calculation of ”sum of squares” is the outstanding exception. [...] Is it not time to seek out novelty in data analysis? Tukey pointed out the need of a new way of analysing data without taking into account as much statistical constraints as usually done, encouraging scientists to draw hypotheses directly from data and not necessarily to reduce the problem to finding a statistical model fitting the data. Data Science is later on introduced, in the late 1990s, as a discipline encompassing the path from data colection to decision making passing through data modeling, combining Statistics and Computer Science with the main objective of helping in the decision making process. Data science is currently experiencing an outstanding time. Technology is allowing computers to go further every day. Decision making is becoming a data driven procces throughout industries and organizations around the whole world. Data itself is often a valuable asset ahead of other tangible assets. Accordingly, Data Science has attracted the interest of people and a lot of new techniques are constantly being developed.
5 Classification and regression problems represent some of the most challenging problems within Data Science. In classification problems we consider given a training sample in the form X={(x1, y1),...,(xn, yn)} ⊆ Rp×{1, . . . , k}, where a set of pfeatures has been measured for a set of nindividuals (x1, . . . , xn), as well as a label (y1, . . . , yn) identified with an element of a finite set of kclasses. The goal of supervised classification is to derive a decision rule DX:Rp→ {1, . . . , k}in order to accurately predict the labels of the forthcoming unobserved data. We call it a multiclass classification problem when the number of classes is greater than two, and we call it a binary classification problem when only two classes are involved. When facing a binary classification problem, the notation on the labels set {1,2}is often replace by {−1,+1}, and hence the classes are refered to as the negative and the positive ones. On its continuous counterpart, in regression problems the label is replaced by a continuous target variable. Therefore, training samples are given in the form X={(x1, y1),...,(xn, yn)} ⊆ Rp×R, where the decision rule DX:Rp→R returns a real number as a prediction. Many different classification and regression techniques can be found in the literature, as for instance, Logistic Regression and Linear Discriminant Analysis (Friedman (2017)), Classification and Regression Trees (Breiman et al. (2017)), Support Vector Machines (Cortes and Vapnik (1995)) and Support Vector Regression (Drucker et al. (1997)), k-nearest neighbor (Cover and Hart (1967)), Random Forest (Breiman (2001)) or Deep Learning methods (Goodfellow et al. (2016)). In addition to supervised problems, there are Data Science problems which lie in the perspective of an unsupervised paradigm, i.e., problems where no target variables are presented. Cluster problems are the most famous ones. The goal on these problems is to create a finite number of groups which present, in a certain sense, homogeonity within the elements of the groups and at the same time heterogenity when elements from different groups are compared. In such a context, the k-means algorithm or hierarchical clusters are the most popular baseline techniques, nevertheless, more complex methods have been studied in the literature as for instance the ordered p-median problem (Boland et al. (2006); Kalcsics et al. (2002); Labb´e et al. (2017)). Throughout the chapters of this thesis the focus is set on developing new strategies for classification and regression problems. However, cluster methodologies are not left out and they are present in most of the chapters. The main elements from where we start this work involve Support Vector Machines, Decision Trees and Fitting Hyperplanes Theory. The following sections of this introduction are devoted to introduce these elements as well as to detail the contributions of this dissertation.
6 Chapter 1. Introduction 1.1 Support Vector Machines 1.1.1 Original problem Support Vector Machines (SVM), introduced by Cortes and Vapnik (1995), is a technique originally proposed to approach binary classification problems by means of solving an optimization problem. Given a training sample {(x1, y1),...,(xn, yn)} ⊆ Rp×{−1,+1}, the goal of SVM is to obtain a hyperplane separating the data (x∈ Rp) into their two different classes (y∈ {−1,+1}). Among all possible hyperplanes that can obtain such a separation between the classes, SVM looks for the one with maximum margin (maximum distance from classes to the separating hyperplane) while minimizing the misclassification errors. Let us denote by Ha hyperplane in Rpin the form H={z∈Rp:ω0z+ω0= 0}for some ω∈Rpand ω0∈R(the vector v0is the result of the transpose operator applied to the vector v∈Rp). This hyperplane will induce a subdivision of the data space Rpinto three regions: the +1 (positive) half-space H+={z:ω0z+ω0>1}, the −1 (negative) halfspace H−={z:ω0z+ω0<−1}and the strip S={z:−1≤ω0z+ω0≤1}. In the SVM model, positive-class observations (y= +1) will be forced to lie on the positive half-space, and the same constraint will be imposed for the negative-class (y=−1) observations on the negative half-space. Moreover, observations will also be penalized if they lie inside the strip. When these constraints are violated for an observation, a classification error is accounted for in the optimization problem. The separation (margin) between classes is computed as the width of the strip S. The points defining the width of the strip, which are known as the support vectors, verify |ω0xi+ω0|= 1. Therefore, recalling that given a point z∈Rp, its Euclidean distance to the hyperplane His |ω0z+ω0| kωk2 , we obtain that the width of the strip is given by 2 kωk2 . As mentioned before, the SVM separating hyperplane will be obtained from an equilibrium of maximizing the separation between classes and minimizing these classification errors. Denoting by di∈R+the classification error of observation i, and by cthe penalty for these errors, the SVM can be formulated as the following Non Linear Problem (NLP): min 1 2kωk2 2+c n X i=1 di(SVM) s.t. yi(ω0xi+ω0)≥1−di,∀i= 1, . . . , n, di≥0,∀i= 1, . . . , n, ω∈Rp, ω0∈R.
1.1. Support Vector Machines 7 In agreement with this, the decision rule for out-of-sample observations, DX, arises naturally according to the region in which an observation is located: DX(z) = (−1 if ω0z+ω0<0, 1 if ω0z+ω0>0. We see in Figure 4.1 a graphical example of the elements involved in a SVM problem. Classes are being represented by the blue and green colors. The black line is the separating hyperplane, and the dashed parallels are the boundary of the strip. The two blue points and the two green points which lie on these lines are the support vectors of this solution. Moreover, the red lines are being used to remark the two kind of errors that can be accounted for in the problem. On the one hand, the small red line is representing a penalization of a blue observation because, despite the fact that it is being well classified, it lies inside the margin, which is desired to be an isolated area. On the other hand, the big red line is representing a penalization of a blue observation which is being wrong classified, since it lies in the side of the green observations. Figure 1.1: SVM example SVM is one of the most famous classification tools. Nevertheless, this is not only due to its original (linear) formulation, but also because of the so called kernel trick. When facing a problem which is not linearly separable, i.e., it is imposible to obtain a hyperplane separating the classes with no classification errors, or if the amount of classification errors obtained when using a separating hyperplane is too big, SVM provides a very useful (kernel) trick in order to deal with these situations. This trick relies on the idea of proyecting data onto a higher dimensional space in which data can be linearly separable, or where at least misclassification errors are not as big as in the original space. This higher dimensional space is called the feature space. In Figure 1.2 we see an example of some points that can not
8 Chapter 1. Introduction be linearly separated in their original data space (left), meanwhile a hyperplane is separating them in the feature space (right). The advantages of this kernel trick Figure 1.2: Original space (left) and feature space with a hyperplane (right) are that, for a transformation ϕ:X→ F from the original space Xto the feature space F, there is no need of knowing the exact transformation ϕ, it is enough with the knowledge of the form of the inner products on the feature space. Moreover, the computational cost of solving the problem in Fis similar to the cost of solving the original problem. These statements can be justified by formulating the dual (Lagrangean) of the problem. In order to obtain this dual problem, we must start from its primal, which can be written as follows, min ω0ω 2+c n X i=1 di(F −SVM) s.t. yi(ω0ϕ(xi) + ω0)≥1−di,∀i= 1, . . . , n, di≥0,∀i= 1, . . . , n, ω∈ F, ω0∈R, and then to apply the Karush-Kuhn-Tucker (KKT) conditions. Therefore, if we consider the following set of dual variables α= (α1, . . . , αn)≥0, λ= (λ1, . . . , λn)≥ 0, and the Lagrangean function L(ω, ω0, α, λ, d) defined as L(ω, ω0, α, λ, d) = ω0ω 2+c n X i=1 di− n X i=1 αiyi(ω0ϕ(xi) + ω0)−(1 −di)− n X i=1 λidi, where d= (d1, . . . , dn) is the vector of the error variables, we obtain by differentiating with respect to ω,ω0and di, and equalazing to zero, the following optimality
1.1. Support Vector Machines 9 conditions: ∂L(ω, ω0, α, λ, d) ∂ω = 0 =⇒ω= n X i=1 αiyiϕ(xi), ∂L(ω, ω0, α, λ, d) ∂ω0 = 0 =⇒ n X i=1 αiyi= 0, ∂L(ω, ω0, α, λ, d) ∂di = 0 =⇒αi=c−λi. According to this, substituting these conditions on L(ω, ω0, α, λ, d), the dual problem is derived as follows: max n X i=1 αi−1 2 n X i=1 n X j=1 αiαjyiyjϕ(xi)0ϕ(xj) s.t. n X i=1 αiyi= 0, αi=c−λi,∀i= 1, . . . , n, 0≤λi,∀i= 1, . . . , n, 0≤αi,∀i= 1, . . . , n, which can be rewritten as max n X i=1 αi−1 2 n X i=1 n X j=1 αiαjyiyjϕ(xi)0ϕ(xj) (Dual −F −SVM) s.t. n X i=1 αiyi= 0, 0≤αi≤c, ∀i= 1, . . . , n. Hence, as we can see, data dependency is shown only through the inner products in F. Therefore, the kernel trick is based on using a known non linear function, K:Rp×Rp→R (xi, xj)→ϕ(xi)0ϕ(xj), instead of computing the inner products of the transformed data, mantaining this way the computational cost of solving the original problem and avoiding us to know
10 Chapter 1. Introduction the explicit form of ϕ. Finally, we conclude this SVM presentation with a brief discussion about different error measures that can be used when solving the problem. The formulations presented before attend to an unbounded continuous meassure of the errors through the divariables, which are referred to as Hinge Loss (HL) measures. SVM can struggle dealing with outlier observations when minimizing HL error measures. This is due to the fact that, in extreme cases, the sum of loads of little classification errors can be smaller than the sum of a few big errors. With the aim of obtaining robust SVM models against outlier observations, Brooks introduced in (Brooks (2011)) two approaches which provide an upper bounded error meassure by means of solving a Mixed Integer Non Linear Programming Problem (MINLP). The first approach is the Ramp Loss SVM (RL-SVM), which maintains the same value as standard SVM on the divariables for observations within the margin, but considers these errors to be constant for wrong classified observations outside the margin. In order to model the problem, we need to introduce the following set of binary variables ξi∈ {0,1} for each of the observations defined as: ξi= 1 if xiis wrong classified outside the margin, 0 otherwise. According to this, the RL-SVM is formulated as follows: min 1 2kωk2 2+c n X i=1 di+ 2 n X i=1 ξi!(RL −SVM) s.t. yi(ω0xi+ω0)≥1−di−Mξi,∀i= 1, . . . , n, 0≤di≤2,∀i= 1, . . . , n, ξi∈ {0,1},∀i= 1, . . . , n, ω∈Rp, ω0∈R, where Mis a big enough constant (see Baldomero-Naranjo et al. (2020) for further details on estimating this constant), and cis the penalization cost. The second approach is the Hard-Loss SVM, which is rather similar to the RL-SVM with the difference that in this method, all classification errors are considered to be constant regardless of the distance from the observations to the separating hyperplane.
1.1. Support Vector Machines 11 Therefore, in this problem ξivariables are slightly modified and defined as: ξi= 1 if xiis wrong classified or is located inside the margin, 0 otherwise, and accordingly, being Ma big enough constant and cthe penalization cost, the problem is formulated as follows: min 1 2kωk2 2+c n X i=1 ξi(HardLoss −SVM) s.t. yi(ω0xi+ω0)≥1−Mξi,∀i= 1, . . . , n, ξi∈ {0,1},∀i= 1, . . . , n, ω∈Rp, ω0∈R. 1.1.2 Multiclass approaches The good perfomance of SVM in binary classification problems has motivated the study of its extension to the multiclass scenario from different points of view. There is no a unique extension to this problem because when dealing with more than two classes, some elements which were unequivocally defined for the biclass problem are left open to interpretation. When facing a multiclass problem, as for instance the one shown in Figure 1.3 of a three class problem in R2, some questions come up naturally: How many hyperplanes do we need to solve the problem? What is the margin in this context? How do we define classification errors? How do we define a decision rule? Figure 1.3: Multiclass example
12 Chapter 1. Introduction There are two main lines to approach the problem. On the one hand, there are sequential methods in which the multiclass problem is divided into a set of binary ones, which are afterwards combined in order to create a solution for the original problem. On the other hand, there are global methods which provide a solution by solving a single optimization problem. In the following paragraphs we describe some of the most known multiclass SVM methods. •Sequential methods The first method we introduce in this section is the One Versus All (OVA) SVM. In such a method, considering kpossible classes, the problem is divided into k binary classification subproblems. Each of these subproblems consists on separating one class from all the other classes. Therefore, positive class is defined as the class to be separated, and the negative one is a fictitious class which is formed gathering the remaining classes. For the multiclass example presented in Figure 1.3, we show in Figure 1.4 the first hyperplane of the OVA method, resulting from separating the yellow class from the other two classes (represented together by the black color). Figure 1.4: OVA - yellow class versus rest This process is repeated until all the hyperplanes Ht={z∈Rp:ω0 tz+ωt0= 0}, ωt∈Rp,ωt0∈R,t= 1, . . . , k, are located, and the decision rule is derived as follows: DX→ {1, . . . , k} x→arg max t=1,...,k ω0 tx+ωt0 The remaining subproblems are shown in Figure 1.5, and the three hyperplanes combined together in Figure 1.6. The main advantage of OVA is that it is a fast method to calculate due to the fact that only kbinary problems are involved. Nevertheless, this method can struggle when the number of classes is large. For a large number of classes, the resulting
1.2. Classification Trees 19 Given a training sample containing nobservations with pfeatures and a label over kpossible classes, OCT is formulated according to minimize an objective function formed by two terms. The first one manages the number of misclassified observations whereas the second one regulates the tree complexity. Complexity is understood as the number of non-empty nodes. On the one hand, in order to control complexity, some binary variables are introduced for each of the branch nodes, δt∈ {0,1}, t ∈τb. δttakes value one if tis a non-empty node, and zero otherwise. On the other hand, for the purpose of accounting for misclassification errors, some variables, Lt∈Z+, are introduced at each of the leaf nodes, t∈τl, where Ltis equal to the number of misclassified observations in leaf t. Therefore, the authors derive the objective function of the problem as follows: 1 ˆ LX t∈τl Lt+cX t∈τb δt, where cis the complexity cost parameter, and ˆ Lis the baseline accuracy (the accuracy obtained predicting the most represented class for the whole dataset). Regarding to OCT constraints, one may start by defining the univariate splits. The authors define the splitting hyperplane coefficients as binary variables, ω∈ {0,1}p, an they impose the sum over these coefficients to be equal to one. Nevertheless, complexity plays its role in this part as well, since δt= 0 means that no observations belong to node t, and therefore no splitting hyperplane should be built in such a node. In order to satisfy the above requirements, the following constraints are added to the problem: p X j=1 ωtj =δt,∀t∈τb, 0≤ω0t≤δt,∀t∈τb. These two constraints allow the model to build a univariate split in node tif and only if δt= 1 (note that the second constraint is valid for normalized data, xi∈[0,1]p, which can be done without loss of generality). Once splits and complexity variables are defined, the authors endow the tree with a hierarchical structure by means of the following constraints: δt≤δp(t),∀ ∈ τb\{1}. According to this, if a split is not applied in node t, further splits will not be allowed in its successors nodes. Another important part of the model is tracking the observations allocation
20 Chapter 1. Introduction within the tree nodes. In order to do this, a set of binary variables, zit, i = 1, . . . , n, t = 1, . . . , T, is added to the problem. These variables take value one if observation iis in node t, and zero otherwise. Moreover, to assure the proper functioning of the model, each observation must be assigned to exactly one leaf node (where classification errors are accounted for), which can be done by imposing the usual assignment constraints: X t∈τl zit = 1,∀i= 1, . . . , n. Furthermore, OCT fixes a minimum number of observations per leaf (Nmin), and so as to regulate this, some binary variables, lt∈ {0,1}, l ∈τl, need to be added to the problem to indicate whether leaf tcontains points (lt= 1), or not (lt= 0). This requirement can be modeled by adding to the problem the following constraints: zit ≤lt,∀i= 1, . . . , n, t ∈τl, n X i=1 zit ≥Nminlt,∀t∈τl. With the above constraints, leaf tis forced to contain more than Nmin observations. The allocation on the leaf nodes is not correctly defined yet. Observations need to satisfy the branch inequalities which lead indeed to leaf tso as to actually be located in t. For such a purpose, the authors start by pointing out the need of the following constraints: ω0 mxi< ω0m−M1(1 −zit),∀i= 1, . . . , n, t ∈τl, m ∈AL(t), ω0 mxi≥ω0m−M2(1 −zit),∀i= 1, . . . , n, t ∈τl, m ∈AR(t), where M1and M2are big enough constants. Note that the first set of constraints use a strict inequality, and therefore the authors use a vector of parameters to transform these into non-strict inequalities. The way of computing this vector as well as the discussion on the values M1= 1 + max(), and M2= 1 is detailed in the paper. In the following we introduce the reformulation of these constraints presented by the authors: ω0 m(xi+)≤ω0m−(1 + max())(1 −zit),∀i= 1, . . . , n, t ∈τl, m ∈AL(t), ω0 mxi≥ω0m−(1 −zit),∀i= 1, . . . , n, t ∈τl, m ∈AR(t). Finally, misclassification errors must be defined. In order to do this, the authors
1.2. Classification Trees 21 define first the variables Nt∈Z+, and Nst ∈Z+,t∈τl, s = 1, . . . , k, as follows: Nst =1 2 n X i=1 (1 + Yis)zit,∀s= 1, . . . , k, t ∈τl, Nt= n X i=1 zit,∀t∈τl, where Yis, which is derived from original data, takes value one if observation i belongs to class s, and minus one otherwise. According to this, Ntis the number of observations located in leaf t, and Nst is the number of observations belonging to class sinside leaf t. Hence, misclassified observations in leaf tcan be calculated as the difference between Ntand the maximum value amongst the Nst variables. Recalling that Ltvariables are being minimized at the objective function, the authors compute these differences by means of the following constraints: Lt≥Nt−Nst −M(1 −cst),∀s= 1, . . . , k, t ∈τl, Lt≤Nt−Nst +Mcst,∀s= 1, . . . , k, t ∈τl, where Mis a big enough constant (nis a valid value for this constant), and cst ∈ {0,1}is a binary variable that takes value one if class sis predicted in leaf t, and zero otherwise. For these last variables, one more set of constraints is required, since exactly one class has to be predicted in leaf t(in case tis a non-empty leaf), which is done imposing the assingment constraints: k X s=1 cst =lt,∀t∈τl. Gathering all the constraints together, OCT is formulated as the following Mixed Integer Problem (MIO): min 1 ˆ LX t∈τl Lt+cX t∈τb δt(OCT)
22 Chapter 1. Introduction s.t. p X j=1 ωtj =δt,∀t∈τb, 0≤ω0t≤δt,∀t∈τb, δt≤δp(t),∀ ∈ τb\{1}, X t∈τl zit = 1,∀i= 1 . . . , n, zit ≤lt,∀i= 1, . . . , n, t ∈τl, n X i=1 zit ≥Nminlt,∀t∈τl, ω0 m(xi+)≤ω0m−(1 + max())(1 −zit),∀i= 1, . . . , n, t ∈τl, m ∈AL(t), ω0 mxi≥ω0m−(1 −zit),∀i= 1, . . . , n, t ∈τl, m ∈AR(t), Nst =1 2 n X i=1 (1 + Yis)zit,∀s= 1, . . . , k, t ∈τl, Nt= n X i=1 zit,∀t∈τl, Lt≥Nt−Nst −n(1 −cst),∀s= 1, . . . , k, t ∈τl, Lt≤Nt−Nst +ncst,∀s= 1, . . . , k, t ∈τl, k X s=1 cst =lt,∀t∈τl, ωt∈ {0,1}p, ω0t∈R, δt∈ {0,1},∀t∈τb, lt∈ {0,1}, Nt, Lt∈R,∀t∈τl, Nst ∈R, cst ∈ {0,1},∀s= 1, . . . , k, t ∈τl, zit ∈ {0,1},∀i= 1, . . . , n. 1.3 Fitting Hyperplanes Theory In this section we summarize some of the elements presented in Blanco et al. (2018) about fitting hyperplanes problems. The problem of locating hyperplanes with respect to a given set of point is well known in Location Theory (LT) (Sch¨obel (2013)). This problem is closely related to another common question in Data Analysis: to study the behavior of a given dataset with respect to a fitting body expressed with an equation of the form f(x) = 0, with x= (X1, . . . , Xp)∈Rp. This last problem reduces to the estimation of the ‘best’ function fthat expresses the relationship between the data or, in the jargon of LT, to the location of the surface f(x) = 0 that minimizes some aggregation function of the distances to these points (see Amaldi et al. (2016); Drezner et al. (2002); Dıaz-B´anez et al. (2004)). In many cases the family
1.3. Fitting Hyperplanes Theory 23 of functions where fbelongs to is fixed and then, the parameters defining such an optimal function must be determined. The family of linear functions is the most widely used. This implies that the above equation is of the form f(x) = ω0+Pp l=1 ωlXl= 0 for ω0, ω1, . . . , ωp∈R. To perform such a fitting, we are given a training sample {x1, . . . , xn} ⊂ Rp, where the goal is to find ˆω0and ˆω= (ˆω1,...,ˆωp) that minimizes some measure of the deviation of the data with respect to the hyperplane they induce, H={z∈ Rp: ˆω0z+ ˆω0= 0}. For a given point x∈Rp, we define the residual with respect to a generic xas a mapping εx:Rp+1 →R+, that maps any set of coefficients (ω0, ω)=(ω0, ω1, . . . , ωp)∈Rp+1, into a measure εx(ω0, ω) that represents the deviation of the given point xfrom the hyperplane with those parameters. The problem of locating a hyperplane for a given set of points {x1, . . . , xn} ⊆ Rpconsists of finding the coefficients minimizing an aggregation function, Φ : Rn→R, of the residuals of all the points. Different choices for the residuals and the aggregation criteria will give, in general, different optimal values for the parameters and thus different properties for the resulting hyperplanes. This problem is not new and some of these criteria, as the minisum, minimax and some other alternatives, have been widely analyzed from a LT perspective (see Megiddo and Tamir (1983); Sch¨obel (1996, 1997, 1998, 2013), amongst others). For a further analysis on the problem of locating a hyperplane to fit our training sample, minimizing different forms of measuring the residuals and their aggregation, we define the Fitting Hyperplane Problem (FHP) as the problem of finding ˆω0and ˆωsuch that: ˆω0,ˆω∈arg min ω0∈R, ω∈RpΦ((ω0, ω)),(FHP) where (ω0, ω)=(ε1(ω0, ω), . . . , εn(ω0, ω)) is the vector of residuals. Note that the difficulty of solving the FHP depends on both the expressions for the residuals and the aggregation criterion Φ. If Φ and εxare linear, the above problem becomes a linear programming problem (LP). In the following we present a multiparametric family of functions, called ordered median functions, introduced in Nickel et al. (2005). In order to do this, let λ1, . . . , λn∈ Rand let ∈Rnbe the vector of residuals of all of the points in the given training sample. Thus, the ordered median functions are derived by means of the aggregation criteria Φ : Rn→R+defined as: Φ() = n X i=1 λiρ (i),1≤ρ < +∞, where (i)∈ {1, . . . , n}is such that (1) ≤ ··· ≤ (n).
24 Chapter 1. Introduction According to this, the choice of the λ-weights will lead to different problems, capturing many of the models proposed in the literature. Most classical models assume that the residuals are defined as the vertical distance (with respect to the last coordinate) from the points to the hyperplane: x(ω0, ω) = xp−ω0 ωp− p−1 X l=1 ωl ωp xl, (assuming that ωp6= 0). Therefore, the difference between them comes from the choice of the aggregation criterion Φ. We show below how some classical methods can be accommodated to this framework. 1. The Least Sum of Squares (LSS) method, credited to Gauss (1877), is the most widely used approach to estimate the coefficients of a linear model due to its simplicity (a closed form for the optimal coefficients is obtained) and its theoretical implications for the inference over the total population. However, somehow restricting hypotheses are required in order to be applied (see, e.g., Giloni and Padberg (2002)). The LSS criterion is defined as the sum of the squares of the residuals, that is: ΦLSS(1, . . . , n) = Pn i=1 2 i,where the residuals iare given by the above expression. The reader may observe that LSS corresponds to the FHP with λ= (1,...,1), ρ= 2 and the vertical distance. 2. The Least Absolute Deviation (LAD) method, introduced by Edgeworth (1887), consists of minimizing the sum of the absolute value of the vertical residuals. Therefore, ΦLAD(1, . . . , n) = Pn i=1 |i|. Note that LAD corresponds to the FHP for λt= (1,...,1) and ρ= 1. 3. The Least Quantile of Squares (LQS), introduced by Bertsimas and Mazumder (2014), is a generalization of the Least Median of Squares (LMS) introduced by Hampel (1975). It also considers vertical distances as residuals, but they are aggregated to minimize the r-quantile of its distribution (rranges in {1, . . . , n}). Hence, ΦLQS(1, . . . , n) = r−quantile(2 1, . . . , 2 n) := 2 (r). This method also fits to the general form of the aggregation criteria considered in this paper. In this case, the LQS hyperplane can be obtained for ρ= 2 and λ= ( (r−1) z }| { 0,...,0,1, (n−r) z }| { 0,...,0). (Observe that LMS hyperplane is also obtained within the same scheme when p= 2 and λ= ( bn 2c z }| { 0,...,0,1, bn 2c z }| { 0,...,0)).
1.3. Fitting Hyperplanes Theory 25 4. The Least Trimmed Sum of Squares (LTS) method was introduced by Rousseeuw Rousseeuw (1984) as a robust alternative to the LSS method, in that it has a high breakdown point. Recall that, intuitively, the breakdown point of an estimator is the proportion of incorrect observations (e.g., arbitrarily large observations) an estimator can handle before giving an incorrect (e.g., arbitrarily large) result. With our notation, it corresponds to choose again as residuals the vertical distance, ρ= 2, and the aggregation criterion ΦLT S(1, . . . , n) = Ph i=1 2 (i)where (i)∈ {1, . . . , n}with (i)≤(i+1) for i= 1, . . . , n −1, and h∈ {1, . . . , n}. The most common choice for his bn 2c, considering the best 50% square residuals. Going back to the general agregation and residual function definitions, on the one hand we have that the function Φ is invariant against permutations of its components and, for non negative lambda weights, a monotone function, ensuring that the ordering of the individual residuals do not affect the overall goodness of the fitting. Moreover, it also implies that a componentwise smaller vector of residuals gives rise to a more accurate fitting. On the other hand, the natural implication of the assumption made about the definition of residuals is that, as expected, the response (projection) of a point on a given hyperplane differs from the classical evaluation. In this setting the response is the closest point, with respect to the distance D, to the hyperplane H, and for a given point x∈Rpthis is, recalling from Mangasarian (1999), x=|ω0x+ω0| kωk∗ , where k·k∗is the dual norm of the one inducing distance D. We should remark at this point that the standard residual (vertical distance) is a distance measure that is not induced by a norm, however its expression can be written in a analogous form and so it fits to the shape of the distances that are considered in this FHP. Gathering all these elements, and considering λ∈Rn, such that λ1≥ ··· ≥ λn≥0, the FHP is equivalent to the following Mathematical Programming problem (Blanco et al. (2018)): min n X j=1 uk+ n X i=1 vi s.t. uk+vi≥λkei,∀i, k = 1, . . . , n, ei≥xi(ω0, ω),∀i= 1, . . . , n, ei, ui, vi∈R+,∀i= 1, . . . , n, ω0∈R, ω ∈Rp. According to this, the first set of constraints control the order over the residuals
26 Chapter 1. Introduction meanwhile the second one ensures proper residual definition. Hence, the distance involved in such a problem will determine whether we obtain a LP (as for instance when using vertical distances, `1or `∞norms) or a NLP (as for instance when using the Euclidean norm). To conlcude this section, observe that when the points in the dataset lie exactly on a hyperplane, H, this hyperplane is always optimal for all versions of the FHP, although for some specific choices of λthe solution may not be unique and different hyperplanes may be alternative optimal. 1.4 Contributions of this thesis In this thesis we analyze different classical Data Science problems through a Location Theory perspective. In particular, we focus on locating hyperplanes to address the problems. Chapter 2 is based on the paper Blanco et al. (2020b). In this chapter we present a multiclass extension of the SVM methodology. In contrast to other SVM multiclass approaches, we propose some MIP and MINLP formulations which maximize the minimum of the margins amongst a set of hyperplanes. Furthermore, our methodology is based on a novel idea consisting on dividing the pfeature space, by means of locating jhyperplanes, into a set of cells that are unequivocally identified with a pattern of signs over a j-tuple. Once the hyperplanes are located, each of the cells is assigned to the most represented class amongst its obervations. Therefore, the decision rule for out-of-sample observations is nothing but predicting the class of the cell in which observations lie in, or the class of the closest non-empty cell. We propose four different models depending on the norm used to measure the margins, `1 or `2, and depending on whether using a bounded or an unbounded error meassure. Moreover, we prove that the kernel trick can be used in all the models. Finally, we report some computational experiment results in which we show the effectiveness of our approaches. Afterwards studying a multiclass classification problem by means of locating mhyperplanes, we study in Chapter 3 the problem of locating mhyperplanes minimizing a distance based error function from observations to the hyperplanes. Although we study a very general problem, allowing one to use many different distances measures and error functions, particular cases result into usual Data Science problems. For instance, when fixing m= 1 and considering the vertical distance, our problem turns into classical linear regression. Therefore, locating mhyperplanes when using the vertical distance can be seen as a m-extension to the classical linear regression problem. Chapter 3 is based on the paper Blanco et al. (2021a). In this work, two approaches are presented to solve the above mentioned problems. On the one hand, exact prob-
1.4. Contributions of this thesis 27 lem formulations are presented by means of some MIP and MINLP. On the other hand, some set partitioning formulations are presented as well as a column generation procedure to solve them. Moreover, convergence results are proven under certain conditions, and some heuristics are also developed so as to speed up the solution processes. Lastly, a discussion of the presented methods performance over an extensive battery of computational experiments is reported. Whereas Chapters 2 and 3 consider a set of mhyperplanes, Chapter 4, which is based on the work Blanco et al. (2020a), focuses on locating a single hyperplane in a binary classification context. This hyperplane is located following a SVM classification principle where some noise might exist within the observation labels set. The motivation of this work comes from different real life situations in which the reliability of the label set in the training sample is not completely truthful. In order to approach this problem, we present three models which combine SVM and cluster ideas. The goal of these models is to create two clusters, as homogeneous as possible, which are separated by a hyperplane, and at the same time the margin betweem them is maximized. A SVM hyperplane separates observations with respect to their classes, nevertheless, we are able to separate observations with respect to the cluters since our models permit observations to be relabeled, i.e., to be treated as if they belong to their opposite class. These methods obtain much higher accuracies than standard SVM when some label noise is added to the training sample, as can be seen in the computational experiments of the work. Optimal Classification Trees have been studied when using splits that do not maximize the separating hyperplanes margin and do not take into account distance based errors. In Chapter 5, which is based on the paper Blanco et al. (2021c), we propose a formulation of an OCT with SVM-based splits, that we call OCTSVM, for binary classification problems. Moreover, in our solution we take advantage of the relabel formulation presented in Chapter 4, and hence our model is also robust against label noise on the training sample. The computational results presented in this chapter show how OCTSVM outperforms OCT, OCT-H and CART methodologies in both noisy label training samples and normal training samples. The multiclass extension of the OCTSVM, that we call MOCTSVM, is later on introduced in Chapter 6. This chapter is based on the work by Blanco et al. (2021b). The tools developed in Chapter 4 are again essential to allow us to formulate the MOCTSVM. This is due to the fact that SVM splits are designed to separate two classes, and therefore we permit the model to create SVM-splits by relabeling observations over two fictitious classes that are separated in the branch nodes, until observations reach the leaf nodes were misclassification errors amongst the original labels are measured. This method is rather similar to the multiclass SVM approach presented in Chapter 2, nevertheless, MOCTSVM follows a hierarchical structure
28 Chapter 1. Introduction when locating the mhyperplanes. According to this, when looking at the global configuration of the hyperplanes in the p-features space, some parts of themselves do not play any action in the decision rule, what can lead to a lower number of cells in a classification problem. This does not necessarily have an impact on the predictive power of the model, but the interpretability of the solution is often easier when the number of cells is low. Finally, we conclude this chapter with some computational experiment results which show that MOCTSVM reports more accurate results than CART, OCT and OCT-H in a battery of real life datasets.
2.2. Preliminaries 35 how an analogous to the kernel trick can be extended to be applied in this model. Section 2.4 describes some heuristic strategies, preprocessing and dimensionality reductions to obtain good quality initial solutions of the MINLP. Finally, in section 2.5 we report our computational results on different real and synthetic datasets, and compare our method with the standard ones for multiclass SVM. 2.2 Preliminaries In this section we introduce the problem under study. Given a training sample X={(x1, y1),...,(xn, yn)} ⊆ Rp× {1, . . . , k}the goal of supervised classification is to find a decision rule to assign labels to data, in order to be applied to out-of- sample data. We assume that a given number, m, of hyperplanes in Rphave to be built to obtain a subdivision of this space into full dimension polyhedral regions that we shall denote as cells. (Here, we would like to mention that the term cell stands for a nonempty intersection of the semispaces induced by the hyperplanes in the considered family). Let us denote by H1,...,Hmthe hyperplanes to be found, which are in the form Hr={z∈Rp:ω0 rz+ωr0= 0}for some ωr∈Rp, ωr0∈R, for r= 1, . . . , m. Each cell induced with such an arrangement of hyperplanes will be then assigned to a label in {1, . . . , k}. In Figure 2.1 we illustrate a subdivision of R2induced by 2 hyperplanes and the labels assigned to each cell. In the left figure, we represent the observations, highlighting the classes with different symbols (stars, circles and squares). In the right figure, two hyperplanes which induce 4 cells are constructed to separate the three classes. Each cell is assigned to a class (north →circles, south →stars, east →stars and west →squares). In this example the subdivision in cells and the assignment of labels reaches a perfect classification on the given observations. Class Stars Class Stars Class Squares Class Circles Figure 2.1: Illustration of a subdivision induced by 2 hyperplanes in R2. From the above, we would like to construct an arrangement of mhyperplanes, H={H1,...,Hm}, determined by ω1, . . . , ωm∈Rp(coefficients) and ω10, . . . , ωm0∈
36 Chapter 2. Multiclass Support Vector Machines R(intercepts) and a decision rule that assigns a single label to each one of the cells in the subdivision of the space induced by such an arrangement. We would like to point out that each cell in the subdivision can be univocally identified with a {−1,+1}-vector in Rm: the `-component of that vector represents the side (positive or negative) with respect to the hyperplane H`where that cell lies in. Definition 2.1 (Suitable Assignment).Given a subdivision Cof Rpinto cells induced by the arrangement of hyperplanes H={H1,...,Hm}in Rp, a function g:{−1,1}m→ {1, . . . , k}is said a suitable assignment, if gunivocally maps cells (equivalently, sign-patterns) to labels in {1, . . . , k}. Observe that a suitable assignment, g, allows us to classify any observation x∈Rpwithin the set of classes {1, . . . , k}, as follows: 1. Identify xwith a sign-pattern: s(x) = (s1(x),...,sm(x)) ∈ {−1,+1}m, where sr(x) = sign(ω0 rx+ωr0) for r= 1, . . . , m. 2. Apply the function gto the sign-patterns: by(x) = g(s(x)) ∈ {1, . . . , k}, is the predicted label of x. The quality of the decision rule is based, on comparing predictions and actual labels on a training sample, but also on maximally separating the classes in order to find good predictions and avoid undesired overfitting. SVM is a particular case of our approach for classifying two-class datasets if m= 1, i.e., a single hyperplane to subdivide the feature space is used. In such a case, signs are in {−1,1}and classes in {1,2}, so whenever there are observations in both classes, the assignment is one-to-one. However, even for biclass instances, if more than one hyperplane is used, one may find better classifiers (we illustrate this behavior with the dataset 2C4N of our computational experiments in Table 2.3). In Figure 2.2, left-and-right, we draw the same dataset of labeled (red and blue) observations and the result of applying a standard SVM (left) and our method with 2 hyperplanes. In that picture one may see that not only the misclassification errors are smaller with two hyperplanes, as expected, but also the separation between classes is larger, improving the predictive power of the classifier. The rationale of our approach is particularly adequate for datasets in which there are several separated “clouds” of observations that belong to the same class. In Figure 2.3, we show two different instances in which, again, the colors indicate the class of the observations. The classes in both instances cannot be appropriately separated using any of the available linear SVM-based methods in the literature since they are based on subdividing the space on class-connected regions. However, we are able to perfectly separate the classes using 5 hyperplanes.
2.2. Preliminaries 37 Figure 2.2: Standard SVM (left) and our approach with 2 hyperplanes (right). In Figure 2.4 we compare our approach and the OVO approach in an instance with 24 observations. In the left figure we show the result of separating the classes with four hyperplanes, reaching a perfect classification on the training sample. In the right figure we show the best linear OVO classifier, in which only 66% of the data were correctly classified. We would also like to highlight that, although nonlinear SVM-approaches may separate the data more conveniently, our approach may help to avoid using kernels and ease the interpretation of the results. Different alternatives could be admissible to justify the rationale of the multiclass classifiers in our framework. To simplify the presentation, we will concentrate on two different models which share the same paradigm but differ in the way they account for misclassification errors. Recall that in SVM-based methods, two criteria are simultaneously optimized when constructing a classifier. On the one hand, a measure of the quality of the decision rule on out-of-sample observations, based on finding a maximum separation between classes; and on the other hand a measure of the misclassification errors for the training set of observations. Both criteria are adequately weighted in order to find a good compromise between the two goals. In what follows we describe how similar measures can be defined in our multiclass classification framework and the way we account them for. 2.2.1 Separation between classes Separation between classes will be measured as it is usual in SVM-based methods. Let ω1, . . . , ωm∈Rpand ω10, . . . , ωm0∈Rbe the coefficients and intercepts of a set of hyperplanes. The distance induced by a norm k·kbetween the shifted hyperplanes H+ r={z∈Rp:ω0 rz+ωr0= 1}and H− r={z∈Rp:ω0 rz+ωr0=−1}is given by
38 Chapter 2. Multiclass Support Vector Machines Figure 2.3: A 5-classes instance classified with our approach using 5 hyperplanes (left) and a 6-classes instance classified with our approach using 5 hyperplanes (right). Figure 2.4: A 4-classes instance classified with our approach using 4 hypeplanes (left) and the same intance classified using the OVO SVM approach (right). 2 kωrk∗, where k·k is a given norm in Rpand k·k∗is its dual norm (see Mangasarian (1999)). Unless explicitly mentioned, we will consider that k · k is the Euclidean norm which dual is also the Euclidean norm. Hence, in order to find globally optimal hyperplanes with maximum separation, we maximize the minimum separation between classes, that is min n2 kω1k,..., 2 kωmko. This measure will conveniently keep the minimum separation between classes as large as possible. Observe that finding the maximum min-separation is equivalent to minimize max{1 2kω1k2,...,1 2kωmk2}. For a given arrangement of hyperplanes, H= {H1,...,Hm}, we will denote by hH(H1,...,Hm) = max{1 2kω1k2,...,1 2kωmk2}. We note in passing that different criteria could have been used to model the separation between classes. For instance, one may consider to maximize the summation of all separations namely m X r=1 2 kωrk2or the inverse of that summation, namely
2.2. Preliminaries 39 2 Pm r=1 kωrk2. However, although mathematically possible, these approaches do not capture the original concept in classical SVM and we have left them to be developed by the interested reader. 2.2.2 Misclassification errors The performance of a classifier on the training set is usually measured with some function of the misclassification errors. Classical SVM with hinge-loss errors use, for non well-classified observations, a penalty proportional to the distance to the side in which they would have been well-classified. Then the overall sum of these errors is minimized. We extend the notion of hinge-loss errors to the multiclass setting as follows. Let H={H1,...,Hm}be an arrangement of hyperplanes and (x, y) a pair observation (x), label (y), with s(x) = (s1(x),...,sm(x)) being the sign-pattern of xwith respect to the hyperplanes in H. Let g:{−1,1}m→ {1, . . . , k}be a suitable assignment. We denote by t(x) = (t1(x),...,tm(x)) the signs of the closest cell to xwhose class by gis y. We will say that (x, y) is wrong-classified with respect to Hrif sr(x)6= tr(x), otherwise it is said that (x, y) is well-classified. In what follows we describe the different error measures (misclassification errors due to different causes) that will be considered for xin order to construct an optimal decision rule. Definition 2.2 (Multiclass In-Margin Hinge-Loss).The multiclass in-margin hingeloss for (x, y)with respect to the hyperplane Hris given as: hIx, y, Hr=(max{0,1−sr(x)·(ω0 rx+ωr0)}if xis well classified through Hr, 0otherwise. Observe that hImodels the error due to observations that although adequately classified with respect to Hr, belong to the margin between the shifted hyperplanes H+ rand H− r. These errors will be zero if the observation is wrong-classified, or if it is well-classified and does not belong to the margin induced by the r-th hyperplane. Definition 2.3 (Multiclass Out-Margin Hinge-Loss).The multiclass out-margin hinge-loss for (x, y)with respect to the hyperplane Hris given as: hOx, y, Hr=(1−tr(x)·(ω0 rx+ωr0)if xis not well classified through Hr, 0otherwise. hOmeasures, for wrong-classified observations, how far they are from being wellclassified. This error is zero whenever an observation is well-classified. Note that if an observation, besides being wrong-classified, belongs to the margin between H+ r
40 Chapter 2. Multiclass Support Vector Machines and H− r, then only hOshould be accounted for. In Figure 2.5 we illustrate the differences between the two types of losses. hI>0 hO>0 Figure 2.5: Illustration of the error measures considered in our approach. 2.3 Mathematical Programming formulations In this section we describe the mathematical optimization models that we propose for the multiclass classification problem. Using the notation introduced in previous sections, the problem can be mathematically stated as follows: min hH(H1,...,Hm) + c1 n X i=1 m X r=1 hIxi, yi,Hr+c2 n X i=1 m X r=1 hOxi, yi,Hr(2.1) s.t. Hris a hyperplane in Rp, for r= 1, . . . , m, where c1and c2are parameters which model the cost of misclassified and striprelated errors. Usually these constants will be considered equal, nevertheless, in practice analyzing different values for them might lead to better results on predictions. A case of interest results considering c2=mc1, i.e., the unitary cost of misclassification errors caused by out-margin observations is mtimes the unitary cost caused by in-margin observations. This method gives a larger penalty to wrongly classified observations, avoiding the calibration of a larger number of parameters. Observe that the problem above consists of finding the arrangement of hyperplanes minimizing a combination of the three quality measures described in the previous section: 1) the maximum margin between classes, 2) the overall sums of the in-margin errors and 3) the out-margin misclassification errors. In what follows,
2.3. Mathematical Programming formulations 41 we describe how the above problem can be rewritten as a mixed integer non linear programming problem by means of adequate decision variables and constraints. Furthermore, the proposed model will consist of a set of continuous and binary variables, a linear objective function, and a set of linear and second order cone constraints. This form allows us to push the model to a commercial solver. First, we describe the variables and constraints needed to model the first term in the objective function. We consider the continuous variables ωr∈Rpand ωr0∈R to represent the coefficients and intercept of the hyperplane Hr, for r= 1, . . . , m. Since there is no distinction between hyperplanes, we can assume, without loss of generality that they are non-decreasingly sorted with respect to the norms of their coefficients, i.e., kω1k≥kω2k≥···≥kωmk. Then, it is straightforward to see that the term hH(H1,...,Hm) can be replaced in the objective function by 1 2kω1k2, once the following set of constraints is included in the model: 1 2kωr−1k2≥1 2kωrk2,∀r= 2, . . . , m. (2.2) For the second term, the in-margin misclassification error, hIxi, yi,Hr, corresponding to the observation (xi, yi) will be identified with the continuous variable eir ≥0, for i= 1, . . . , n,r= 1, . . . , m. Observe that to properly determine these errors, one has to determine whether the observation xiis well-classified or not with respect to the rth hyperplane. In order to do that we need to introduce some binary variables. First, we consider the two following sets of binary variables: tir =(1 if ω0 rxi+ωr0≥0, 0 otherwise, and zis =(1 if iis assigned to class s, 0 otherwise, for i= 1, . . . , n,r= 1, . . . , m,s= 1, . . . , k. The t-variables model the sign-pattern of the observations, while the z-variables give the allocation profile of observations to classes. As mentioned above, the classification rule is based on assigning signpatterns to classes. The adequate definition of the t-variables is assured with the following constraints: ω0 rxi+ωr0≥ −T(1 −tir),∀i= 1, . . . , n, r = 1, . . . , m, (2.3) ω0 rxi+ωr0≤Ttir,∀i= 1, . . . , n, r = 1, . . . , m, (2.4) where Tis a big enough constant. Observe that Tcan be accurately estimated based on the data set under consideration. Furthermore, the following constraints assure the adequate relationships between
42 Chapter 2. Multiclass Support Vector Machines the variables: k X s=1 zis = 1,∀i= 1, . . . , n, (2.5) kzi−zjk1≤2kti−tjk1,∀i= 1, . . . , n, =1, . . . , n. (2.6) Observe that (2.5) enforces that a single class is assigned to each observation while (2.6) assures that the assignments of two observations must coincide if their signpatterns are the same. Additionally, the set of z-variables determines whether an observation is well-classified. Indeed, let δi∈ {0,1}kbe defined as δis = 1 if yi=s and 0 otherwise. (Observe that δiis the binary encoding of the class of the ith observation.) Then, ξi=1 2kzi−δik1∈ {0,1}assumes the value zero if and only if the observation iis well-classified, i.e., ξi=(0 if iis well-classified, 1 otherwise. Now, we will model whether the ith observation is well-classified or not, with respect to the rth hyperplane. Observe that the measure of how far is a wrongclassified observation from being well-classified, needs a further analysis. One may have a wrong-classified observation and several training observations in its same class. We assume that the error for this observation is the misclassification error with respect to the closest cell for which there are well-classified observations in its class. Thus, we need to model the decision on the well-classified representative observation for a wrong-classified observation. In Figure 2.6, we illustrate this type of misclassification errors. The observation xiis wrong-classified but the misclassification error of xi, in case xjis chosen as its representative (well-classified) observation, is 0 with respect to hyperplane H1(note that both xiand xjare in the same side of H1), whereas the misclassification error with respect to H2is h. Observe that h is the distance between xiand the shifted hyperplane defining the halfspace where xjlies in. We consider the following set of binary variables: hij = 1 if xj, which is well classified and verifies yj=yi, is the representative of xi in its closest cell through hyperplanes, 0 otherwise.
2.3. Mathematical Programming formulations 43 These variables require to impose the following constraints: n X j=1: yi=yj hij = 1,∀i= 1, . . . , n, (2.7) ξj+hij ≤1,∀i, j = 1, . . . , n, (yi=yj),(2.8) hii = 1 −ξi,∀i= 1, . . . , n. (2.9) The first set of constraints, (2.7), imposes a single assignment between observations belonging to the same class. Constraints (2.8) avoid choosing wrong-classified representative observations. The set of constraints (2.9) enforces well-classified observations to be represented by themselves. h H1 H2 xj xi Figure 2.6: Illustration of the wrong-classification errors. With these variables, we can model the in-margin errors by means of the following constraints: ω0 rxi+ωr0≥1−eir −T(3 −tir −tjr −hij),∀r= 1, . . . , m, (2.10) ω0 rxi+ωr0≤ −1 + eir +T(1 + tir +tjr −hij),∀r= 1, . . . , m. (2.11) These constraints model, by using the sign-patterns given by t, that, eir = max{0,min{1,1− sr(x)(ω0 rxi+ωr0}}. Note that the constraints are active if either tir =tjr =hij = 1, i.e., if the well-classified observation xjis the representative observation for xiand both are in the positive side of the rth-hyperplane; or tir =tjr = 0 and hij = 1, i.e., if the well-classified observation xjis the representative observation for xiand both are in the negative side of the rth-hyperplane. Thus, constraints (2.10) and (2.11) adequately model the in-margin errors for all observations . Furthermore, because of (2.3) and (2.4), and those described above, the variables eir always take values
44 Chapter 2. Multiclass Support Vector Machines smaller than or equal to 1. Finally, the third addend, the out-margin errors, will be modeled through the continuous variables dir ≥0, for i= 1, . . . , n,r= 1, . . . , m. With the set of variables described above, the out-margin misclassification errors can be adequately modeled through the following constraints: dir ≥1−ω0 rxi−ωr0−T(2 + tir −tjr −hij),∀i, j = 1, . . . , n(yi=yj), r = 1, . . . , m, (2.12) dir ≥1 + ω0 rxi+ωr0−T(2 −tir +tjr −hij),∀i, j = 1, . . . , n(yi=yj), r = 1, . . . , m. (2.13) Constraints (2.12) are active only if tir = 0 and tjr =hij = 1, that is, if xjis a well-classified observation in the positive side of Hr, while xiis wrong-classified in the negative side of Hrbeing xjthe representative observation for xi(note that if xi is well-classified then hii = 1 by (2.9) and then, the constraint cannot be activated). The second set of constraints, namely (2.13), can be analogously justified in terms of the negative side of Hr. The main difference of these constraints with respect to (2.10) and (2.11) is that (2.12) and (2.13) are active only if xiis wrong-classified. According to the above constraints, a misclassified observation xiis penalized in two ways with respect to each hyperplane Hr. In case that xiis well-classified with respect to Hr, but it belongs to the margin, then eir = 1 −sign(ω0 rxi+ωr0)(ωt rxi+ ωr0)≤1 and dir = 0 (tir =tjr). Otherwise, if xiis wrong-classified with respect to Hr, then dir = 1 −sign(ω0 rxi+ωr0)(ωt rxj+ωr0)≥1 and eir = 0 (hij = 1 and tir 6=tjr). We illustrate the rationale of the proposed constraints on the data drawn in Figure 2.7. Observe that A is not correctly classified since it lies within a cell in which the blue-class is not assigned. Suppose that B, a well-classified observation, is the representative of A (hAB = 1), then the model would have to penalize two types of errors. The first one with respect to H2. If we suppose tB2= 1, then tA2= 0, leading to an activation on constraint (2.12) being dA2>0. On the other hand, even though A is well-classified with respect to H1, we also have to penalize its margin violation. Again, if we assume tB1= 1, then tA1= 1, what would activate constraint (2.10) being eA1>0.
2.3. Mathematical Programming formulations 51 Moreover, we need another problem that accounts for the first part of Φ. Φ1(h, z, t, ξ) = min ω1,ω10,e,d 1 2kω1k2+c1 n X i=1 ei1+c2 n X i=1 di1 s.t. −ω0 1xi−ω10 ≤κ1(ti1),∀i= 1, . . . , n, ω0 1xi+ω10 ≤κ2(ti1),∀i= 1, . . . , n, −ω0 1xi−ω10 −ei1≤κ3(u+ ij1),∀i, j = 1, . . . , n, ω0 1xi+ω10 −ei1≤κ3(u− ij1),∀i, j = 1, . . . , n, −ω0 1xi−ω10 −di1≤κ3(q+ ij1),∀i, j = 1, . . . , n (yi=yj),(SP1) ω0 1xi+ω10 −di1≤κ3(q− ij1),∀i, j = 1, . . . , n (yi=yj), −di1≤0,∀i= 1, . . . , n, −ei1≤0,∀i= 1, . . . , n, ω1∈Rp, ω10 ∈R. Thus, using the above notation, (MCSVM) is equivalent to the following problem: min h,z,t,ξ,ω1 Φ1(h, z, t, ξ) + m X r=2 Φr(h, z, t, ξ, ω1) s.t. (2.5) −(2.9) hij ∈ {0,1},∀i, j = 1, . . . , n, tir ∈ {0,1},∀i= 1, . . . , n, r = 1, . . . , m, zis ∈ {0,1},∀i= 1, . . . , n, s = 1, . . . , k, ξi∈ {0,1},∀= 1, . . . , n. Observe that the above problem only accounts for ω1and the binary variables. Once they are fixed can be plugged into the SPr,r= 1, . . . , m subproblems and then it allows one to find the optimal values of the continuous variables. The elements κ1(tir), κ2(tir), κ3(u+ ijr), κ3(u− ijr), κ3(q+ ijr), and κ3(q− ijr) are fixed constants once the binary variables are fixed. In order to solve the problem for a fixed set of binary variables we can proceed recursively, solving independently (SPr) for all r= 2, . . . , m, and then combining their solutions with (SP1). Therefore, to get that goal we apply Lagrangean duality to obtain and exact dual reformulation. Indeed, relaxing the constraints of (SPr) with dual multipliers µr≥0 (assuming that µ0 r6= 0) and denoting by αir =αir =µ1 ir −µ2 ir +X j:yi=yj (µ3 ijr − µ4 ijr +µ5 ijr −µ6 ijr), for all i= 1, . . . , n and r= 1, . . . , m; after some derivations that
52 Chapter 2. Multiclass Support Vector Machines can be followed in Lemma 2.1, one can check that evaluating Φ for given values of the binary variables, can be obtained solving the continuous optimization problem (JSP). (All details can be found in Lemma 2.1.) Next, we analyze how the evaluation of the function φ, given in (φval), depends on the original data. In order to do that we find the optimal solutions of (JSP). Dualizing the constraints that correspond to Φ1with multipliers µ1≥0 and those where the variables Γrappear, with dual multipliers γr, the Lagrangean function of the problem (JSP) results in: L(ω1, ω01, d, e;µ) = 1 2kω1k2(1 − m X r=2 γrµ0 r) + c1 n X i=1 ei1+c2 n X i=1 di1+ m X r=2 Γr(1 −γr) + m X r=2 γr 2µ0 r n X i,j=1 αirαjrx0 ixj+X i,j:u+ ijr=0 µ3 ijr +X i,j:u− ijr=0 µ4 ijr +X i,j:q+ ijr=0 µ5 ijr +X i,j:q− ijr=0 µ6 ijr + n X i=1 µ1 i1−ω0 1xi−ω10 −κ1(ti1) + n X i=1 µ2 i1ω0 1xi+ω10 −κ2(ti1) + n X i=1 µ3 ij1−ω0 1xi−ω10 −ei1−κ3(¯u+ ij1) + n X i=1 µ4 ij1ω0 1xi+ω10 −ei1−κ3(¯u− ij1) + n X i,j=1 µ5 ij1−dir −ω0 1xi−ω10 −κ3(¯q+ ij1) + n X i,j=1 µ6 ij1−dir +ω0 1xi+ω10 −κ3(¯q− ij1) + n X i=1 µ7 i1(−di1) + n X i=1 µ8 i1(−ei1), and its KKT optimality conditions reduce to: •(1 −Pm r=2 γrµ0 r)ω1= n X i=1 µ1 i1−µ2 i1+X j:yi=yj (µ3 ij1−µ4 ij1+µ5 ij1−µ6 ij1)xi.
2.3. Mathematical Programming formulations 53 • n X i=1 µ1 i1−µ2 i1+X j:yi=yj (µ3 ij1−µ4 ij1+µ5 ij1−µ6 ij1)= 0. •X j:yi=yj (µ3 ij1+µ4 ij1) + µ8 i1=c1,∀i= 1, . . . , n. •X j:yi=yj (µ5 ijr +µ6 ijr) + µ7 ir =c2,∀i= 1, . . . , n. •1−γr= 0, ∀r= 2, . . . , m. •µ≥0, γ≥0. Using the same notation as before, αi1=µ1 i1−µ2 i1+X j:yi=yj (µ3 ij1−µ4 ij1+µ5 ij1−µ6 ij1), for all i= 1, . . . , n, and simplifying the expressions using the KKT conditions and the complementary slackness conditions we get that the strong dual of (JSP) is:
54 Chapter 2. Multiclass Support Vector Machines max ω,ω.0,d,e;µ 1 2(1 −Pm r=2 µ0 r) n X i,j=1 αi1αj1x0 ixj+X i,j:u+ ij1=0 µ3 ijr +X i,j:u− ij1=0 µ4 ij1 +X i,j:q+ ij1=0 µ5 ijr +X i,j:q− ij1=0 µ6 ij1+ m X r=2 1 2µ0 r n X i,j=1 αirαjrx0 ixj +X i,j:u+ ijr=0 µ3 ijr +X i,j:u− ijr=0 µ4 ijr +X i,j:q+ ijr=0 µ5 ijr +X i,j:q− ijr=0 µ6 ijr s.t. n X i=1 αir = 0,∀r= 1, . . . , m, (DJSP) n X i=1 αi1xi= (1 − m X r=2 µ0 r)ω1, n X i=1 αirxi=µ0 rωr,∀r≥2, X j:u+ ijr=0 µ3 ijr +X j:u− ijr=0 µ4 ijr ≤c1,∀i= 1, . . . , n, X j:q+ ijr=0 µ5 ijr +X j:q− ijr=0 µ6 ijr ≤c2,∀i= 1, . . . , n, r = 1, . . . , m, µ1 ir = 0,∀i= 1, . . . , n, r = 1, . . . , m, with ¯ tir = 0, µ2 ir = 0,∀i= 1, . . . , n, r = 1, . . . , m, with ¯ tir = 1, αir =µ1 ir −µ2 ir +X i,j:u+ ij1=0 µ3 ijr −X i,j:u− ij1=0 µ4 ijr +X i,j:q−+ij1=0 µ5 ijr −X i,j:q− ij1=0 µ6 ijr, ∀i= 1, . . . , n, µr≥0,∈R,∀r= 1, . . . , m. Note that the objective function of (DJSP) only depends of the x-input data through the inner products x0 ixjfor i, j = 1, . . . , n, and also the expressions of ω1 from the dual variables is given as: (1 − m X r=2 µ0 r)ω1= n X i=1 αi1xi. The dependence of ω1with other ωrin the primal formulation is only through the nonincreasing sorted values of kω1k,...,kωmkin which we now that the largest value is kω1k. Thus, solving the dual problem (DJSP) allows us to determine all the optimal hyperplanes ωr, for all r= 1, . . . , m. In case a transformation ϕis
2.3. Mathematical Programming formulations 55 performed to the input data, the dependence of the data in the problem will be through the inner products ϕ(xi)0ϕ(xj) for all i, j = 1, . . . , n, and the kernel Theory can be applied. Lemma 2.1. The evaluation of Φfor given values of the binary variables ¯ t,¯ ξ,¯ h and ¯zcan be obtained solving the following continuous optimization problem: ˆ Φ(h, z, t, ξ) = min (1 2kω1k2+c1 n X i=1 ei1+c2 n X i=1 di1+ m X r=2 Γr) s.t. −µ0 r 2kω1k2+1 2µ0 r n X i,j=1 αirαjrx0 ixj +X i,j:u+ ijr=0 µ3 ijr +X i,j:u− ijr=0 µ4 ijr +X i,j:q+ ijr=0 µ5 ijr +X i,j:q− ijr=0 µ6 ijr +c1 ≤Γr,∀r≥2, −ω0 1xi−ω10 ≤κ1(ti1),∀i= 1, . . . , n, ω0 1xi+ω10 ≤κ2(ti1),∀i= 1, . . . , n, −ω0 1xi−ω10 −ei1≤κ3(u+ ij1),∀i, j = 1, . . . , n, ω0 1xi+ω10 −ei1≤κ3(u− ij1),∀i, j = 1, . . . , n, −ω0 1xi−ω10 −di1≤κ3(q+ ij1),∀i, j = 1, . . . , n, (yi=yj),(JSP) ω0 1xi+ω10 −di1≤κ3(q− ij1),∀i, j = 1, . . . , n, (yi=yj), −di1≤0, i = 1, . . . , n, −ei1≤0, i = 1, . . . , n, ω1∈Rp, ω10 ∈R, X j:u+ ijr=0 µ3 ijr +X j:u− ijr=0 µ4 ijr +µ8 ir ≥c1, i = 1, . . . , n, X j:q+ ijr=0 µ5 ijr +X j:q− ijr=0 µ6 ijr ≤c2,∀i= 1, . . . , n, r = 1, . . . , m, µ1 ir = 0,∀i= 1, . . . , n, r = 1, . . . , m, with ¯ tir = 0, µ2 ir = 0,∀i= 1, . . . , n, r = 1, . . . , m, with ¯ tir = 1, αir =µ1 ir −µ2 ir +X i,j:u+ ij1=0 µ3 ijr −X i,j:u− ij1=0 µ4 ijr +X i,j:q+ ij1=0 µ5 ijr −X i,j:q− ij1=0 µ6 ijr, ∀i= 1, . . . , n, µ≥0.
56 Chapter 2. Multiclass Support Vector Machines where µr≥0are the dual multipliers of the constraints in Problem (SPr). Proof. In order to prove the result, in what follows we derive the optimality conditions of the Lagrangean dual of (SPr) for r= 2, . . . , m. Relaxing the constraints of (SPr) with dual multipliers µr≥0, the Lagrangean function for given values of the binary variables ¯ t,¯ ξ,¯ hand ¯z(and consequently for values ¯u+, ¯u−, ¯q+and ¯q−) is: Lr(ωr, ω0r, d, e;µ) = c1 n X i=1 eir +c2 n X i=1 dir +µ0 r1 2kωrk2−1 2kω1k2 + n X i=1 µ1 ir−ω0 rxi−ωr0−κ1(tir)+X i=1 µ2 irω0 rxi+ωr0−κ2(tir) + n X i,j=1 µ3 ijr−ω0 rxi−ωr0−eir −κ3(¯u+ ijr)+ n X i,j=1 µ4 ijrω0 rxi+ωr0−eir −κ3(¯u− ijr) + n X i,j=1 µ5 ijr−dir −ω0 rxi−ωr0−κ3(¯q+ ijr)+ n X i,j=1 µ6 ijr−dir +ω0 rxi+ωr0−κ4(¯q− ijr) + n X i=1 µ7 ir(−dir) + n X i=1 µ8 ir(−eir). Therefore, the KKT optimality conditions read as: •µ0 rωr= n X i=1 µ1 ir −µ2 ir +X j:yi=yj (µ3 ijr −µ4 ijr +µ5 ijr −µ6 ijr)xi. • n X i=1 µ1 ir −µ2 ir +X j:yi=yj (µ3 ijr −µ4 ijr +µ5 ijr −µ6 ijr)= 0. •X j:yi=yj (µ3 ijr +µ4 ijr) + µ8 ir =c1,∀i= 1, . . . , n. •X j:yi=yj (µ5 ijr +µ6 ijr) + µ7 ir =c2,∀i= 1, . . . , n. •µr≥0. First of all, we observe that if kωrk<kω1kat optimality then µ0 r= 0 and actually, we do not have to consider the corresponding constraint nor the addend µ0 r 2(kωrk2−kω1k2) in the Lagrangean function. Hence, denoting by αir =µ1 ir −µ2 ir + X j:yi=yj (µ3 ijr −µ4 ijr +µ5 ijr −µ6 ijr), for all i= 1, . . . , n, and assuming that µ0 r6= 0, the
2.3. Mathematical Programming formulations 57 dual of (SPr) reads as: max ωr,ωr0,d,e;µ−µ0 r 2kω1k2+1 2µ0 r n X i,j=1 αirαjrx0 ixj− n X i=1 µ1 irκ1(¯ tir)− n X i=1 µ2 irκ2(¯ tir)− n X i=1 µ3 irκ3(¯u+ ijr) − n X i=1 µ4 irκ4(¯u− ijr)− n X i,j=1 µ5 ijrκ5(¯q+ ijr)− n X i,j=1 µ6 ijrκ6(¯q− ijr) s.t. n X i=1 αirxi=µ0 rωr, n X i=1 αir = 0, X j:yi=yj (µ3 ijr +µ4 ijr)≤c1,∀i= 1, . . . , n, X j:yi=yj (µ5 ijr +µ6 ijr)≤c2,∀i= 1, . . . , n, αir =µ1 ir −µ2 ir +X i,j:u+ ij1=0 µ3 ijr −X i,j:u− ij1=0 µ4 ijr +X i,j:q+ ij1=0 µ5 ijr −X i,j:q− ij1=0 µ6 ijr, ∀i= 1, . . . , n, µr≥0. Let us simplify further the expressions above. We observe that: X i µ1 irκ1(¯ tir) + X i µ2 irκ2(¯ tir)=0, − n X i,j=1 µ3 ijrκ3(¯u+ ijr)− n X i,j=1 µ4 ijrκ3(¯u− ijr) = X i,j: ¯u+ ijr=0 µ3 ijr +X i,j: ¯u− ijr=0 µ4 ijr, and −X i,j µ5 ijrκ3(¯q+ ijr)−X i,j µ6 ijrκ4(¯q− ijr) = X i,j: ¯q+ ijr=0 µ5 ijr +X i,j: ¯q− ijr=0 µ6 ijr.
58 Chapter 2. Multiclass Support Vector Machines Using those equations and substituting the problem becomes: max ωr,ω0,d,e;µ−µ0 r 2kω1k2+1 2µ0 r n X i,j=1 αirαjrx0 ixj+X i,j:u+ ijr=0 µ3 ijr +X i,j:u− ijr=0 µ4 ijr +X i,j:q+ ijr=0 µ5 ijr +X i,j:q− ijr=0 µ6 ijr s.t. n X i=1 αirxi=µ0 rωr, n X i=1 αir = 0,∀r= 1, . . . , m, X j:u+ ijr=0 µ3 ijr +X j:u− ijr=0 µ4 ijr ≤c1,∀i= 1, . . . , n, X j:q+ ijr=0 µ5 ijr +X j:q− ijr=0 µ6 ijr ≤c2,∀i= 1, . . . , n, r = 1, . . . , m, (DSPr) µ1 ir = 0,∀i= 1, . . . , n, r = 1, . . . , m with ¯ tir = 0, µ2 ir = 0,∀i= 1, . . . , n, r = 1, . . . , m with ¯ tir = 1, αir =µ1 ir −µ2 ir +X i,j:u+ ij1=0 µ3 ijr −X i,j:u− ij1=0 µ4 ijr +X i,j:q+ ij1=0 µ5 ijr −X i,j:q− ij1=0 µ6 ijr, ∀i= 1, . . . , n, µ≥0. Using the strong duality in all the subproblems (SPr) for r= 2, . . . , m, we can obtain the following expansion of the join subproblem (SPr) that allows one the evaluation of φdefined by φval.
2.4. Math-heuristic approach 59 ˆ Φ(h, z, t, ξ) = min (1 2kω1k2+c1 n X i=1 ei1+c2 n X i=1 di1) + max ωr,ω0,d,e;µ m X r=2 −µ0 r 2kω1k2+1 2µ0 r n X i,j=1 αirαjrx0 ixj+X i,j:u+ ijr=0 µ3 ijr +X i,j:u− ijr=0 µ4 ijr +X i,j:q+ ijr=0 µ5 ijr +X i,j:q− ijr=0 µ6 ijr s.t. −ωt 1xi−ω10 ≤κ1(ti1),∀i= 1, . . . , n, ω0 1xi+ω10 ≤κ2(ti1),∀i= 1, . . . , n, −ω0 1xi−ω10 −ei1≤κ3(u+ ij1),∀i, j = 1, . . . , n, ω0 1xi+ω10 −ei1≤κ3(u− ij1),∀i, j = 1, . . . , n, −ω0 1xi−ω10 −di1≤κ3(q+ ij1),∀i, j = 1, . . . , n (yi=yj), ω0 1xi+ω10 −di1≤κ3(q− ij1),∀i, j = 1, . . . , n (yi=yj), −di1≤0,∀i= 1, . . . , n, −ei1≤0,∀i= 1, . . . , n, ω1∈Rp, ω10 ∈R, X j:u+ ijr=0 µ3 ijr +X j:u− ijr=0 µ4 ijr ≤c1,∀i= 1, . . . , n, X j:q+ ijr=0 µ5 ijr +X j:q− ijr=0 µ6 ijr ≤c2,∀i= 1, . . . , n, r = 1, . . . , m µ1 ir = 0,∀i= 1,...,1, r = 1, . . . , m, with ¯ tir = 0, µ2 ir = 0,∀i= 1,...,1, r = 1, . . . , m, with ¯ tir = 1, αir =µ1 ir −µ2 ir +X i,j:u+ ij1=0 µ3 ijr −X i,j:u− ij1=0 µ4 ijr +X i,j:q+ ij1=0 µ5 ijr −X i,j:q− ij1=0 µ6 ijr, ∀i= 1, . . . , n, µ≥0. The usual transformation of the maximum in the objective function gives rise to the equivalent reformulation of the above problem as (JSP). 2.4 Math-heuristic approach As mentioned in previous section, the computational burden for solving (MCSVM), that is a mixed integer non linear programming problem (in which the nonlinearities
60 Chapter 2. Multiclass Support Vector Machines come from the norm minimization in the objective function), is the combination of the discrete aspects and the non-linearities in the model. In this section we provide some heuristic strategies that allow us to cut down the computational effort by fixing some of the variables. It will also provide good-quality initial feasible solutions when solving, exactly, (MCSVM) using a commercial solver. Two different strategies are provided. The first one consists of applying a variable fixing strategy to reduce the number of h-variables in the model (originally n2). The second approach consists of fixing to zero some of the z-variables. These nk variables allow us to model assignments between observations and classes. The proposed method is a math-heuristic approach, since after applying the adequate dimensionality reductions, Problem (MCSVM) (or (MCSVMRL)) has to be solved. Also, although our strategies do not ensure any kind of optimality certificate, they produce a very good performance as will be shown in our computational experiments. Observe that when classifying datasets, the measure of the efficiency of a decision rule, as ours, is usually assessed by means of the accuracy of the classification on out-of-sample data, whereas the objective value of the proposed model is just an approximated measure of such an accuracy which cannot be computed only with the training data. Algorithm 1 A math-heuristic approach. 1. Apply dimensionality reductions test based on algorithms 2 and 3. 2. Find an initial solution generating kseparating hyperplanes. 3. Solve problem (MCSVM) (or (MCSVMRL)) up to a prescribed accuracy for the train data. In what follows we describe two strategies to reduce the dimensionality of the problem. These approaches are based on applying clustering techniques to the data. The methods are sensible to the number of clusters. For determining this parameter, we run a hierarchical clustering method, using as termination criterion a given squared Euclidean distance between the observations and their centroids. 2.4.1 Reducing the h-variables Our first strategy comes from the fact that for a given observation xi, there may be several possible choices for hij to assume the value one with the same final result. Recall that hij could be equal to one whenever xjis a well-classified observation in the same class as xi. The errors eir and dir are then computed by using the class of xjbut not the observation xjitself. Thus, if a set of well-classified observations of the same class is close enough, only one of them can be the representative element of the group. In order to illustrate the procedure, we show in Figure 2.8 (left) a
Chapter 3 Multisource hyperplanes location problem to fitting set of points 68
70 Chapter 3. Multisource hyperplanes location problem to fitting set of points In this chapter we study the problem of locating a given number of hyperplanes minimizing an objective function of the closest distances from a set of points. We propose a general framework for the problem in which norm-based distances between points and hyperplanes are aggregated by means of ordered median functions. A compact Mixed Integer Linear (or Non Linear) programming formulation is presented for the problem and also an extended set partitioning formulation with an exponential number of variables is derived. We develop a column generation procedure embedded within a branch-and-price algorithm for solving the problem by adequately performing its preprocessing, pricing and branching. We also analyze geometrically the optimal solutions of the problem, deriving properties which are exploited to generate initial solutions for the proposed algorithms. Finally, the results of an extensive computational experience are reported. The issue of scalability is also addressed showing theoretical upper bounds on the errors assumed by replacing the original datasets by aggregated versions. 3.1 Introduction Location Analysis deals with the determination of the optimal positions of facilities to satisfy the demand of a set of customers. The problems analyzed in the field are diverse but can be usually classified as: Discrete Location problems (DLP) and Continuous Location problems (CLP). In the first family, a set of potential facilities is previously given and the goal is to select, among them, the optimal ones under one or more criteria. The main tools for solving these problems come from Discrete Optimization, or more precisely, from Integer Linear Programming. In the second family of problems, the facilities have to be located in a continuous space and then, convex analysis and global optimization tools are needed to solve the problems. The most popular problem in the latter family is the Weber problem (Weber, 1909) in which a single point-facility has to be positioned on the plane so as to minimize the overall sum of the (Euclidean) distances to a set of (planar) demand points. The applications of both types of location problems are vast. DLP are more common in the location of physical facilities (as ATMs, supermarkets, stations, etc), while CLP are more useful when locating facilities in telecommunication networks (as wifi routers, servers, etc) or even to provide the sets of potential facilities for a DLP. In this chapter we study a problem that falls into the family of CLP. More specifically, we focus on the determination of optimal hyperplanes fitting a given finite set of demand points. The location of a single hyperplane is a classical problem that has been addressed in different fields. On the one hand, this problem clearly extends the classical Weber problem, but where instead of locating zero-dimensional facilities one looks for locating higher dimensional structures. On the other hand, in
3.1. Introduction 71 Statistics and Data Analysis, the determination of a hyperplane minimizing the sum of squares of vertical residuals is key for estimating a multivariate linear regression model using the Least Sum of Squares (Gauss (1877)). One can also find useful applications, both in Location Science and Data Analysis, for the problem of finding optimal hyperplanes fitting a set of points. For instance, Espejo and Rodr´ıguez-Ch´ıa (2011) deals with the location of a rapid transit line on the plane to be used as an alternative transportation mean. Analogously, SVM is also based on constructing a hyperplane minimizing certain loss functions of the distances to a given set of points. Scanning the literature one can find that most of the attention has been devoted to finding hyperplanes with any of the following assumptions (see e.g., (Martini and Sch¨obel, 1998; Sch¨obel, 2013, 2003, 2015; Martini and Sch¨obel, 2001; Brimberg et al., 2002, 2003; Blanco et al., 2018)): (a) the problem is embedded on the plane; (b) a single hyperplane has to be located; (c) the vertical distance between each point and the hyperplane is considered; or (d) the residuals are aggregated by the sum or the maximum operators. Our goal here is to study a generalization of this problem in which, we construct simultaneously a given number, m, of hyperplanes in any finite dimensional space, Rp, by minimizing a rather general globalizing function, an ordered median function, of the residuals from the points to the fitting bodies. Ordered median functions aggregate the set of distances from the demand points to their closest hyperplanes (residuals) by means of a sorting, weighting averaging operation: distances are sorted and then their weighted sum is performed. The sum and maximum functions can be easily represented as ordered median functions with adequate choices of the weights inducing the median and center objective functions. Also, the k-centrum (sum of the k-th largest distances) or the cent-dian (convex combination of the sum and the max criteria) can be cast within this family of functions. In addition, different point-to-hyperplane norm-based distances are considered as a measure of the residuals of the fitting. Thus, this chapter naturally extends the analysis performed in Blanco et al. (2018) where the location of a single ordered median hyperplane was studied. As in the classical Weber problem (Weber, 1909), the extension from the location of one to several facilities (the so-called multisource problem) is not trivial (Blanco et al., 2016). Actually, while the classical single-facility point location problem with standard distances (`p, polyhedral, etc) can be formulated as a Second Order Cone programming problem (Blanco et al., 2014) (being then polynomially solvable), its multisource version becomes a non-convex NP-hard problem (Blanco et al., 2016). In the case of locating hyperplanes, the situation is even harder, since the location of a single hyperplane is, in general, an NP-hard problem (see Blanco et al. (2018)) whose exact solution can only be obtained for vertical and polyhedral norm based residuals.
72 Chapter 3. Multisource hyperplanes location problem to fitting set of points The problem considered in this chapter is not fully new although, in our opinion, it has not been fully analyzed and there is room for further improvement. In particular, similar problems have been analyzed from the Data Analysis field, and different names have been adopted. In the so-called Clusterwise Linear Regression (CLR) problem, a set of observations is provided and the goal is to cluster them by means of the sum of the squared residuals of several multivariate regression models (Sp¨ath, 1982; Hennig, 1999; Carbonneau et al., 2014; Park et al., 2017; Gitman et al., 2018). In (Bertsimas and Shioda, 2007), classification and regression are simultaneously performed, and also clustering by classical linear regression approaches. Finally, in (Bradley and Mangasarian, 2000), the clusters are constructed based on the closest distances to optimal hyperplanes in a given p-dimensional space. In the so-called Piecewise Linear Regression problem, a dependent variable is partitioned into jintervals and it adjusts linear bodies to each of them (see (McGee and Carleton, 1970)). However, only local search heuristic algorithms have been proposed for these problems, alternating clustering and regression techniques sequentially. Carbonneau et al. (2014) present a column generation algorithm for the (planar) clusterwise regression problem with sum of squared residuals which combined with some heuristic strategies outperforms previous results in the literature. Moreover, Park et al. (2017) generalized the clusterwise regression problem by allowing each entity to have more than one observation and propose an exact Mathematical Programming-based approach relying on column generation, and several heuristics. In this chapter we provide a general framework for the simultaneous location of several hyperplanes to fit a data set using Mathematical Programming tools. We formulate the problem by using general norm-based error measures of the distance from points to hyperplanes and ordered median functions to aggregate the residuals. This approach generalizes both the standard multisource regression (Carbonneau et al., 2014; Park et al., 2017) and also the more recent proposal for the m= 1 case (Blanco et al., 2018). The rest of the chapter is organized as follows. In section 3.2 we introduce the problem and fix the notation for the rest of the sections. This section also contains two illustrative examples taken from the literature. Section 3.3 is devoted to a first compact formulation for the problem. This formulation has a polynomial number of variables and constraints but its performance is not always good since it has a large integrality gap. For that reason, in section 3.4 we develop an alternative formulation with an exponential number of variables that is solved (exactly, for vertical and polyhedral-norm based residuals) within a branch-and-price (B&P) algorithm using column generation at each node of the branching tree. This section describes all the elements of this B&P: initialization, pricing (exact and heuristic) and branching. Section 3.5 reports our computational results based on two different datasets: the
3.2. Multisource location of hyperplanes 73 classical 50 points dataset by Eilon et al. (1971) and another synthetic dataset randomly generated. Section 3.6 is devoted to explore scalability issues and finally section 3.7 draws some conclusions and future extensions. 3.2 Multisource location of hyperplanes In this section we describe the problem under study and fix some set-related notation used in the chapter. We are given a set of nobservations/demand points (denoted as points from the rest of the chapter) in Rp,{x1, . . . , xn} ⊂ Rpand m∈Z+(p > 0). Our goal is to find mhyperplanes in Rpthat minimize an objective function of the closest distances from points to hyperplanes. We denote the index sets of demand points and hyperplanes by I={1, . . . , n}and J={1, . . . , m}, respectively. Given ω∈Rp and ω0∈R, we denote by H(ω, ω0) = {z∈Rp:ω0z+ω0= 0},the hyperplane in Rpwith coefficients ωand intercept ω0. Several elements are involved when finding the best mhyperplanes to fit a set of demand points. In what follows we describe them: •Residuals: The point-to-hyperplane measure of closeness. Given a demand point x= (x1,. . .,xp)∈Rpand a hyperplane H(ω, ω0), how far/close is the point from the hyperplane? The classical fitting methods use the so-called vertical-distance measure, which given a reference coordinate, say the p-th, computes the deviation xp+ω0 ωp+Pp−1 `=1 ω` ωpx`for all i∈I, whenever ωp6= 0. However, it has been already proposed that the use of more general distance measures based on norms may be advisable. In particular, some authors (see Blanco et al. (2020c)) have shown the usefulness of norm-based distances, such as polyhedral, or `p-distances (p≥1). Amongst them, we mention, for their importance, the Manhattan (`1-norm), the Tchebyshev (`∞-norm) or the Euclidean (`2-norm) distances. Thus, for a point x∈Rpand a hyperplane H(ω, ω0), we consider the residual from xto H(ω, ω0) as: εx(ω, ω0)=Dx, H(ω, ω0):= min{D(x, y) : y∈ H(ω, ω0)}, where D is a norm-based distance or the vertical distance in Rp(see Mangasarian (1999); Blanco et al. (2018) for further details on this projection). •Allocation Rule: Given a set of hyperplanes and a point, once the residuals to each of the hyperplanes are calculated, one has to allocate the point to a single hyperplane. Different alternatives can be considered, as the allocation to the
74 Chapter 3. Multisource hyperplanes location problem to fitting set of points closest or the furthest hyperplane. In our framework we assume, as usual in Location Analysis, that each point is allocated to the hyperplane with the smallest residual, i.e., for a point x∈Rpand an arrangement of hyperplanes H=nH(ωj, ω0j) : j∈Jo, the final residual point-to-hyperplanes is computed as: εxH= min j∈Jεx(ωj, ω0j), and the hyperplane, H(ω, ω0j), reaching such a minimum is the one where xis allocated (in case of ties among more than one hyperplane, a random assignment is performed). •Aggregation of Residuals: Given a set of points and an arrangement of hyperplanes, once the residuals are computed with respect to the arrangement, and in order to find the mhyperplanes that best fit the ndata points, a global measure of goodness must be chosen for aggregating the residuals. The classical aggregation functions are the sum or maximum of squared residuals. Most of these criteria can be cast within the framework of the family of ordered median aggregation criteria. More explicitly, given x1, . . . , xn∈Rp, an arrangement of hyperplanes H=nH(ωj, ω0j) : j∈Jo, and λ∈Rn +(with λ1≥ ··· ≥ λn≥0), we can extend the definition presented in Chapter 1 and consider λ-ordered median function defined as: OMλ(ε1, . . . , εn) = n X i=1 λieρ (i),1≤ρ < ∞(OMF) where e(1), . . . , e(n)are defined such that e(i)∈ {εx1(H), . . . , εxn(H)}for all i∈Iand e(1) ≥ ··· ≥ e(n). Observe that particular cases of ordered median problem for ρ= 1 (we consider ρ= 1 for the rest of the chapter) are the sum (λi= 1, i= 1, . . . , n), the maximum (λ1= 1, λi= 0, i6= 1), the k-centrum (λi= 1, i= 1, . . . , k,λj= 0, j > k) or the ν-centdian, a convex combination of sum and max criterion (λ1= 1, λi=ν, i = 2, . . . , n), for 0 < ν < 1. Summarizing all the above considerations, the Multisource Ordered Median Fitting Hyperplanes Problem (MOMFHP) can be stated as the problem of finding ω1, . . . , ωm∈Rpand ω01, . . . , ω0m∈Rsolving the following optimization problem: min n X i=1 λie(i)(MOMFHP0) s.t. ei≥min j∈Jεxi(ωj, ω0j),∀i= 1 . . . , n, ω∈Rp, ω0j∈R,∀j= 1, . . . , m, ei≥0,∀i= 1 . . . , n.
3.2. Multisource location of hyperplanes 75 where eirepresents the residual for the i-th point in the data set, for all i= 1 ...,n. MOMFHP appears when different trends or clouds have to be differentiated on the demand points, and then, different hyperplanes want to be use to fitting to the points, such that the global error assumed, when the points are allocated to their closest hyperplanes, is as small as possible. In Figure 3.1 we illustrate a set of demand points in the plane which could be clustered into three groups according to different linear trends which are drawn in gray color. Figure 3.1: Illustration of a feasible solution of our problem for a set of demand points. In the following example we illustrate the problem under analysis in two classical instances. Example 3.1. In the seminal paper by McGee and Carleton (1970), the authors illustrate the Clusterwise Linear Regression method with two instances. The first instance, (Quandt, 1958), consists of 20 points on the plane, {x1, . . . , x20}generated as follows: xi2= 2.5+0.7xi1+i, i = 1,...,12, and xi2= 5 + 0.5xi1+i, i = 13,...,20, where is randomly generated as a univariate normal distribution with mean 0and standard deviation 1. We run our model with this dataset choosing as residuals the `1-norm projection of the data onto the hyperplanes, and four different ordered median crite-
76 Chapter 3. Multisource hyperplanes location problem to fitting set of points ria: Weber, Center, dn 2e-Centrum (λ= ( dn 2e z }| { 1,...,1,0,...,0)) and 0.9-centdian (λ= (1,0.9,...,0.9)). The results are shown in Figure 3.2. Figure 3.2: Lines obtained for Quandt dataset for `1-norm residuals and different criteria. McGee and Carleton (1970) also analyzed a real instance, the Boston dataset. It was motivated by the fact that regional stock exchanges were hurt by the abolition of give-ups in 1968. The model tries to analyze the dollar volume of sales on the Boston Stock Exchange with respect to dollar volumes for the New York and American Stock Exchanges, based on a dataset with 35 monthly observations from January 1967 to November 1969. One can observe, in the results shown in (Figure 3.3), that our models are able to adequately cast the trends of these observations. 3.3 A compact formulation for (MOMFHP0) In this section we provide a Mathematical Programming formulation for (MOMFHP0). The main components which involve decisions in this problem, and that have to be adequately included in a suitable formulation, are the representation of general norm-based residuals and the aggregation of residuals using an ordered median function. We describe here how to incorporate all these elements into a Mathematical Programming formulation which in many cases is suitable to be solved with any of the available MIP/MINLP solvers.
3.4. Set partitioning formulation 83 The set partitioning formulation for (MOMFHP0) is: min X i∈IX S⊆I:S3i λie(i) SyS(3.21) s.t. X S⊆I yS=m, (3.22) X S⊂I: S3i yS= 1,∀i∈I, (3.23) yS∈ {0,1},∀S⊆I. (3.24) where e(i) Sis the i-th element in the sorted sequence of (active) residuals. In the above formulation the objective function computes the ordered median aggregation of the residuals (each demand point iallocated to its cluster S). Constraint (3.22) assures that mclusters have to be computed and constraints (3.23) that each observation belongs to a single cluster. In the same manner that we formulate the ordered median objective function in the compact formulation we can equivalently reformulate the problem above as follows: min n X k=1 uk+ n X i=1 vi(3.25) s.t. uk+vi≥λkX S3k ei SyS,∀i, k ∈I, (3.26) X S⊆I yS=m, (3.27) X S⊂I: S3i yS= 1,∀i∈I, (3.28) yS∈ {0,1},∀S⊂I, uk, vi∈R,∀i, k ∈I. This problem will be referred to as the Master Problem. The problem above, although adequately solves the problem of finding the m hyperplanes once the optimal clusters are computed, has an exponential number of variables (and coefficients to incorporate to constraints (3.26)), and then it is hard to solve unless the number of points is very small. Thus, we propose a column generation (CG) approach for solving, efficiently, the problem above by adding new variables to the model as needed and not considering all of them at the same time. A pseudocode indicating the general procedure is shown in Algorithm 4. Initially, a (small) subset of the y-variables is considered (those indexed by the sets in S0) and a relaxed version of the problem is solved with only these variables.
84 Chapter 3. Multisource hyperplanes location problem to fitting set of points It implies to compute the amounts ei Sfor all S∈ S0and i∈S. Next, it has to be checked whether the optimality condition is satisfied. If it is not the case, a new set of variables is found and added to the relaxed problem and the procedure is repeated. Algorithm 4 General Scheme for the CG approach. Data: {x1, . . . , xn} ⊆ Rp,m∈Z+(m > 0), λ1≥ ··· ≥ λn≥0. 1. Preprocessing: Compute a set of initial solutions for the problem S0= {S1, . . . , SK}with Sk⊆Ifor all k= 1, . . . , K. 2. Relaxed Master: Solve the relaxed master problem: min n X k=1 uk+ n X i=1 vi s.t. uk+vi≥λkX S3k ei SyS,∀i, k ∈I, X S∈S0 yS=m, (RMP) X S∈S0: S3i yS= 1,∀i∈I, 0≤yS≤1,∀S∈ S0. 3. New Columns : Check if new columns have to be added to (RMP). If: Optimality is satisfied C∗={S⊆I:y∗ S= 1}. Else: Update S0with the new columns and go to 2. End Result:{H(ωS, ω0S) : S∈ C∗}. The crucial steps in the implementation of the CG approach are the following: 1. Preprocessing: Generation of initial feasible solutions in the form of initial clusters (and their costs). This step may be improved by the theoretical properties verified by the corresponding optimal hyperplanes. We have implemented different initial solutions based on properties of the optimal solution of median and center hyperplanes (see Section 3.4.2). 2. Pricing: As already mentioned, in set partitioning problems, instead of solving initially the problem with the whole set of exponentially many variables, the variables have to be incorporated on-the-fly by solving adequate pricing
3.4. Set partitioning formulation 85 subproblems derived from previously computed solutions until the optimality of the solution is guaranteed. 3. Branching: The rule that creates new nodes of the branch and bound tree when a fractional solution is found at a node of the search tree. In this problem, we have adapted the Ryan-and-Foster branching complemented by a secondary ad-hoc branching in some special situations. 3.4.1 Preprocessing In the preprocesing phase, we generate different types of initial solutions, which implies the initialization of the CG algorithm with a given set of variables. We consider different types of initial solutions derived from the construction of hyperplanes fitting the sets of points. First, to initialize the pool of columns, S0, we randomly generate hyperplanes passing through poriginal points. Among the various strategies compared, we have eventually implemented one that performs completions with p−2 points of all possible couples of original points. This strategy augment n(n−1) 2new variables into the pool. In addition, we also augment to the pool the best single hyperplane that fits all the points, assuring that the problem is feasible at the root node of the branch-and-price tree. Finally, apart from the above initial columns, we also charge an initial heuristic solution (in the y-variables) so as to have a good upper bound in our branch-and-price algorithm. Our algorithm chooses at random mmutually disjoint subsets of ppoints and finds the hyperplanes determined by those msets of ppoints. Next, we perform a 1-interchange heuristic generating a new hyperplane that replaces, one at a time, one of those currently considered in the configuration until the first iteration where no improvement is possible. The incumbent set of hyperplanes and their corresponding allocations is considered an initial solution that is loaded into the solver. 3.4.2 Median and center optimal hyperplanes We have used the following properties to build the initial solutions of our CG approach since they determine optimal hyperplanes for specific objective functions, see e.g., Sch¨obel (2003). Lemma 3.1. The following properties are verified: 1. Weak incidence property: There exists an optimal median hyperplane passing through paffinely independent points. 2. Pseudo-halving property: Every optimal median hyperplane, H(ω∗, ω∗ 0)verifies #ni∈I:xi∈ H−(ω∗, ω∗ 0)o≤n 2and #ni∈I:xi∈ H+(ω∗, ω∗ 0)o≤n 2.
86 Chapter 3. Multisource hyperplanes location problem to fitting set of points 3. Weak blockedness property: There exists an optimal center hyperplane that is at maximum distance from p+ 1 of the points. 4. Parallel facets property: There exists an optimal center hyperplane that is parallel to a facet of the convex hull of the given points. For the more general ordered median objective function, we have proved the following result that characterizes the ordered median hyperplanes. In what follows, we derive a novel result for these hyperplanes that will be useful in the preprocessing phase of our CG approach. Let us introduce the following notation: •Let Bbe the subdivision of Rp+1 induced by the following arrangement of hyperplanes: Bab ij =n(ω, ω0)∈Rp+1 :a(ω0xi+ω0) = b(ω0xj+ω0)o,∀i, j ∈I, a, b ∈ {−1,1}. •Let Sbe the subdivision of of Rp+1 induced by the following arrangement of hyperplanes Si=n(ω, ω0)∈Rp+1 :ω0xi+ω0= 0o,∀i∈I. Lemma 3.2. If H(ω, ω0)is an optimal ordered median hyperplane then (ω, ω0)is an extreme point of a cell in the subdivision of Rp+1 induced by the intersection B∩S. Proof. For a given hyperplane H(ω, ω0), let us consider the objective function of the problem, namely Pi∈Iλie(i), where ei=D(H(ω, ω0), xi). It is clear that within each cell of the subdivision Bthe sorting of the residuals does not change. In addition, in each cell of the subdivision Sthe sign of ω0xi+ω0 is either positive or negative (but does not change) for each i∈I. Therefore, if C∈ B ∩ S is a cell in the subdivision induced by B ∩ S, there is a permutation σthat fixes the sorting of the residuals and also a constant vector (sign(ω0x1+ ω0), . . . , sign(ω0xn+ω0)) ∈ {−1,1}nsuch that X i∈I λie(i)=X i∈I λi sign(ω0xi+ω0)(ω0xi+ω0) kωk∗ =Pi∈Iλisign(ω0xi+ω0)(ω0xi+ω0) kωk∗ . The above function is the ratio of a non-negative linear function and a convex function, then it is quasiconcave provided that (ω, ω0)∈C. Therefore, it attains its minima at the extreme points of this region. Hence, if H(ω, ω0) is an optimal ordered median hyperplane (ω, ω0) must be an extreme point of some of those cells.
3.4. Set partitioning formulation 87 The above result allows us to interpret optimal ordered median hyperplanes also in terms of a geometrical description as those that meet pconditions between the following cases: i) passing through points xi, i ∈I, and ii) being at the same distance of two points xi, xj, i, j ∈I. Optimal ordered median hyperplanes must also satisfy, for some k= 1, . . . , p, the following property: it contains kpoints xi, i ∈Iand it is at the same distance from p−kpairs xi, xji, j ∈I. In our computational results we have computed the initial solutions and the initial pool of variables for the objective functions of type median, k-centrum and centdian, using the weak incidence property, whereas for the center objective function we use the weak blockedness property. 3.4.3 Pricing problem Certifying optimality in a CG approach avoiding the inclusion of all the columns into the relaxed master problem, (RMP), requires testing whether a new tentative column must be added to the problem. In case no new candidates are added to the master problem, the optimality is guaranteed, otherwise, one should add the new columns and repeat the process (Step 3 in Algorithm 4). Searching for new columns to be added to the model will be performed by looking at the dual formulation of the set partitioning formulation. Let γbe the dual variable for constraint (3.27), φithe dual variables for constraints (3.28) and δik the dual variables for constraints (3.26). Then, the dual of the Master Problem is the following: max −mγ + n X i=1 φi s.t. n X k=1 δik = 1,∀i∈I, n X i=1 δik = 1,∀k∈I, −X i∈S n X k=1 λkei Sδik −γ+X i∈S φi≤0,∀S⊆ S0, δik, γ, φi≥0. Hence, for any S⊂ S0, since ySdoes not appear in the objective function, the reduced cost for variable ySis: ¯eS=γ−X i∈S φi+X i∈S n X k=1 λkδikei S.
88 Chapter 3. Multisource hyperplanes location problem to fitting set of points Then, given an optimal dual solution (γ∗, φ∗, δ∗), and considering the binary variables wi=(1 if the ith observation is chosen for the set Sindexing the new column, 0 otherwise, the pricing problem is to choose the subset Swith minimum reduced cost, i.e., to solve: min − n X i=1 φ∗ iwi+γ∗+ n X i=1 c∗ iri s.t. zi≥εxi(ω, ω0)),∀i∈I, ri≥zi−M(1 −wi),∀i∈I, wi∈ {0,1},∀i∈I, zi, ri≥0,∀i∈I, ω∈Rp, ω0∈R, where c∗ i=Pn k=1 λkδ∗ ik,∀= 1, . . . , n. If the optimal value of this problem is negative, the new column yˆ Sis added to the pool, where ˆ S={i:wi= 1}, since its reduced cost in the (RMP) is negative, and thus, it improves the objective function of the master problem. Otherwise, optimality is certified and we are finished. Heuristic pricing The exact pricing routine described above is an NP-hard problem and thus in general, it takes time finding new columns to be added to the pool or to certify optimality of the reduced master problem. This last task cannot be avoided, provided that we design an exact solution algorithm. Nevertheless, in many occasions finding promising new variables can be done at very low computational time resorting to heuristic schemes. In our problem, we propose to test hyperplanes chosen from a discrete set of potential candidates. To do so, we set a p-dimensional grid on the normalized space of ω0and ωcoefficients. Each point represents a hyperplane to be tested. Once the candidate (ω, ω0) is chosen we determine which set of points Sis going to be added to the new variable yS. This is done with a simple greedy rule: choose those points with negative reduced cost with respect to H(ω, ω0). If after this process we find a hyperplane that produces a negative reduced cost, we add this new column to the pool. Otherwise, we proceed with the exact pricer. This scheme speeds up the search for new columns without loosing the exactness of
3.4. Set partitioning formulation 89 the whole algorithm. 3.4.4 Branching The set partitioning formulation of the MOMFHP0is often not solved at the root node, in contrast with what is stated in Park et al. (2017). Thus, some branching strategy must be implemented to cope with the branch and bound search. Ryan- Foster (R-F) is one of the most popular techniques for branching in set partitioning problems (see Ryan and Foster (1981)). If a fractional solution is reached at a node, R-F creates two new branches as follows: Given to elements i1, i2∈I, they may never go together on a set in the whole branch, or they may always go together, i.e., if one of them belongs to a set S, the other one must also be included in S. To implement this branching, we can take advantage of the wivariables defined on the previous section for the pricing subproblem, to easily adapt this way of branching in our problem, by means of the following constraints: A) wi1+wi2= 1 ensuring that elements i1and i2are not assigned to the same hyperplane. B) wi1=wi2ensuring that elements i1and i2are assigned to the same hyperplane. Moreover, in our formulation there is a new case in which, despite the fact of having fractional solutions on a node, we will not create new branches following the R-F rule. This fact is motivated because in our problem may appear different columns (different y-variables) but being associated to the same set S, although possibly with different hyperplanes. Let S⊆Ibe a subset of points and let y1 S, ..., yq Sbe fractional variables for the same set Salthough with different hyperplanes H(ωi, ωi 0), i= 1 . . . q, with q > 1, such that q X i=1 yi S= 1. If there are no more fractional variables, or the rest of the fractional variables of the node satisfy the same conditions for some other subsets of points, we cannot apply R-F rule and either the node need not be branched (see Theorem 3.2 and Remark 3.4) or a different branching strategy must be implemented in these cases. Without loss of generality, we will describe the new branching for the case in which two fractional variables, y1 Sand y2 S, with hyperplanes H(ω1, ω1 0) and H(ω2, ω2 0), are obtained in a node for the same subset S. In this situation, the new branching rule that we propose creates three new branches as follows: 1. A branch where y1 S= 1 meaning that this variable will be in the solution in this branch. This is easily implemented in the pricing routine since it amounts to avoid considering the elements in Sin any further column in that branch
90 Chapter 3. Multisource hyperplanes location problem to fitting set of points because they are already in the set Swhich is part of the solution. Therefore, it suffices to fix the variables wi= 0,∀i∈Sin all the subproblems in the branch. 2. Analogously, it creates another branch where y2 S= 1. 3. The third branch sets y1 S=y2 S= 0. This branch represents the case in which none of the original fractional solutions are part of the integer solution. Once again, this can be enforced by adding the following constraints to the pricing subproblems of the branch: |S|−X i∈S wi!+|S|− n X i=1 wi+ d X `=1 |ωj `−ω∗ `|+|ωj 0−ω∗ 0|!·M≥1, for j=1,2, and Ma big enough constant, where ω∗and ω∗ 0define the new hyperplane H(ω∗, ω∗ 0). These constraints will make the problem infeasible if and only if all the individuals in S, and only the individuals of S, belong to the new solution, and moreover, the new solution provides a hyperplane H(ω∗, ω∗ 0) that is equal to H(ω1, ω1 0) or H(ω2, ω2 0). The alternative branching may be necessary in case of using general norm based residuals. Nevertheless, as we show below, the situation is simpler using vertical distance residuals. Theorem 3.2. Ryan and Foster branching is enough in the set partitioning formulation of (MOMFHP0)for the vertical distance residuals: If for a subset of points S⊆I, there exists a fractional solution 0< yS<1with #{yS6= 0}>1, at a node of the branch and bound tree then there exists an explicit solution that combines these variables to obtain a single one satisfying yS= 1 (#{yS6= 0}= 1). Proof. Let us consider a subset of points S⊆I. At a fractional node, we can have two possible scenarios: 1) # {yS6= 0}= 1, hence, it would exist a single hyperplane (a facility) H(ω, ω0) that would serve the points of S, and hence, R-F branching is enough, and 2) # {yS6= 0}>1. This latter case needs a further analysis since it may seem as if more than one facility would need to be involved to optimally serve the points in S. Without loss of generality we can assume # {yS6= 0}= 2 (a case with # {yS6= 0}> 2 can be treated sequentially by smaller problems with two solutions). In this situation there are two variables, y1 Sand y2 S, with values σand 1 −σ,σ∈(0,1), so that y1 S+y2 S= 1. These variables are represented by two hyperplanes H(ω1, ω1 0) and H(ω2, ω2 0), where the cost of a point i∈Swith coordinates x∈Rp,ei, is given by ei=σD(x, H(ω1, ω1 0)) + (1 −σ)D(x, H(ω2, ω2 0)). We prove that the hyperplane H(ω∗, ω∗ 0) defined as
3.4. Set partitioning formulation 91 H(ω∗, ω∗ 0) = z∈Rp:σ(ω1z+ω1 0) + (1 −σ)(ω2z+ω2 0)=0, satisfies that D(x, H(ω∗, ω∗ 0)) ≤σD(x, H(ω1, ω1 0)) + (1 −σ)D(x, H(ω2, ω2 0)) for vertical distance residuals, and this would mean that there exists a unique hyperplane that optimally serves all the points in S. Therefore, considering y∗ S, no further branching is required. If we consider the normalized hyperplanes H(ω1, ω1 0), and H(ω2, ω2 0), such that ω1 p=ω2 p=−1, then ω∗ p=−1, the vertical distance from x to H(ω∗, ω∗ 0) is Dv(x, H(ω∗, ω∗ 0)) = xp−ω∗ 0− p−1 X `=1 ω∗ `x`. Hence, Dv(x, H(ω∗, ω∗ 0)) = xp−ω∗ 0− p−1 X `=1 ω∗ `x` =σxp+ (1 −σ)xp−(σω1 0+ (1 −σ)ω2 0)− p−1 X `=1 (σω1 `+ (1 −σ)ω2 `)x` ≤σxp−ω1 0− p−1 X `=1 ω1 `x`+ (1 −σ)xp−ω2 0− p−1 X `=1 ω2 `x` =σDv(x, H(ω1, ω1 0)) + (1 −σ)Dv(x, H(ω2, ω2 0)). Remark 3.4. We prove that under mild conditions, R-F branching is also enough for the `1-norm based residuals. Without loss of generality, assume that there is a solution with two fractional variables y1 Sand y2 S, with values σand (1 −σ), and corresponding hyperplanes H(ω1, ω1 0)and H(ω2, ω2 0). Let us define the set SP =nj:|ω1 j|=|ω2 j|= 1, j = 1, . . . , po. Hence, if SP 6=∅, RF-branching is enough for the `1-norm residuals. Indeed, let ˆ be an index such that |ω1 ˆ |=|ω2 ˆ |= 1 and define sign(ˆ) = (1if ω1 ˆ ·ω2 ˆ = 1 −1if ω1 ˆ ·ω2 ˆ =−1. It is clear that for any ˆ ∈SP,ω∗ ˆ =σω1+ sign(ˆ)(1 −σ)ω2satisfies kω∗ ˆ k∞= 1. Consider for any ˆ ∈SP the hyperplane H(ω∗ ˆ , ω∗ 0) = z∈Rp:σ(ω1z+ω1 0) + sign(ˆ)(1 −σ)(ω2z+ω2 0)=0, then for any individual i∈Swith coordinates xi∈Rp, taking into account that
92 Chapter 3. Multisource hyperplanes location problem to fitting set of points ||ω1||∞=||ω2||∞= 1, we obtain that D`1(xi,H(ω∗, ω∗ 0)) = |ω∗ 0+Pp `=1 ω∗ `xi`| ||ω∗||∞ ≤σ|ω1 0+Pp `=1 ω1 `xi`| ||ω1||∞ + (1 −σ)|ω2 0+Pp `=1 ω2 `xi`| ||ω2||∞ =σD`1(xi,H(ω1, ω1 0)) + (1 −σ)D`1(xi,H(ω2, ω2 0)). and hence, y∗ S= 1 is an optimal solution for the problem. 3.5 Computational results A series of computational experiments has been performed in order to test the two proposed methodologies. We consider two different sets of instances, one based on Eilon et al. (1971) dataset and another on synthetic data. For all of them we solve (MOMFHP0) for four different objective functions: Weber (W), Center (C), dn 2e- Centrum (K) (λ= ( dn 2e z }| { 1,...,1,0,...,0)) and 0.9-Centdian (D) (λ= (1,0.9,...,0.9)) and with the two proposed approaches: the compact approach based on formulation (MOMFHP) and with the branch-and-price methodology. We test the performance of the algorithms on two different types of residuals: `1-norm based residuals and absolute value vertical distance residuals. The models were coded in C and solved with SCIP v.6.0.1 using as optimization solver CPLEX 12.8 in a Mac OS El Capitan with a Core i7 CPU clocked at 2.8 GHz and 16GB of RAM memory. A time limit of 5 hours was fixed for all the instances. It is well-known in the field of location analysis that continuous multifacility ordered median problems are very difficult to solve and already problems of moderate sizes (n= 50 demand points) can not often be solved to optimality (see e.g., Blanco et al. (2016)). The same or even a harder behavior should be expected here since these problems introduce a new degree of difficulty in the representation of general distance based residuals. 3.5.1 Eilon et al. (1971) dataset First, we tested our approach on instances based on the classical planar 50-points dataset provided by Eilon et al. (1971). We randomly generate five instances from such a dataset with sizes n∈ {20,30,40,45}and the entire complete original instance with n= 50. We run the models for m∈ {2,5}hyperplanes. The average results obtained for these instances are shown in tables 3.1 and 3.2. There, for each combination of n(size of the instance), m(number of hyperplanes to be located) and type (ordered median objective function to be minimized), we provide both for
3.6. Scalability: bounding the error in aggregation procedures 99 possible choice can be substituting the set of original demand points by the centroids obtained by any of the available clustering techniques. In any case, when solving (MOMFHP0) for X0instead of using Xone incurs in aggregation errors. Let Hbe the optimal arrangement of mhyperplanes for the problem and e= (e1, . . . , en) with ei=εxiH, for i∈I, the residuals with respect to H. Analogously, let H0be the optimal arrangement for the demand points in X0and e0the vector of residuals. Theorem 3.3. Let T= max i=1,...,n D(xi, x0 i). Then, the following relation holds: |OMλ(e0)−OMλ(e)| ≤ 2OMλ(T,...,T).(3.29) Proof. First of all, observe that, based on the triangular inequality, for any H εxi(H)≤εx0 i(H) + D(xi, x0 i),∀i= 1, . . . , n. Let us also consider the vector t= (D(x1, x0 1),...,D(xn, x0 n)) of distances from the original points in Xto their corresponding points in X0and denote by ˜ e= (εx0 1(H), . . . , εx0 n(H)). Since the function OM is non-decreasing monotone and sublinear, it follows that: OMλ(e)≤OMλ(˜ e+t)≤OMλ(˜ e) + OMλ(t). Hence, since T≥D(xi, x0 i) for all i= 1, . . . , n, we get that: |OMλ(e)−OMλ(˜ e)| ≤ OMλ(T, . . . , T). From the above inequality we can apply (Geoffrion, 1977, Theorem 5) to conclude that |OMλ(e0)−OMλ(e)| ≤ 2OMλ(T,...,T). The difference considered in the above theorem is the excess due to the implementation of an approximate solution based on the reduced model with data set X0rather than the correct optimal solution for the original data in the larger set X. This result allows us to scale our CG algorithm to problems of any size using aggregation techniques and providing estimates on the deviation from the optimal value. We illustrate the application of the above result including the percent error obtained aggregating to 20 points some of our random problems with 50 points by the 20-mean clustering technique. As one can see in Table 3.5 the percent errors
100 Chapter 3. Multisource hyperplanes location problem to fitting set of points error (%) mtype p= 2 p= 3 2 W 3.48 2.78 K -17.84 -16.52 C 4.40 5.90 D 3.59 2.87 5 W -0.16 1.21 K -9.17 -2.99 C 6.60 20.86 D 0.16 -11.16 Table 3.5: % aggregation errors for 50 points problems and vertical distance. are small. Observe that in some cases they are even negative, for problems that were not solved to optimality, and where the hyperplanes obtained by aggregating points, once evaluated on the actual 50 points, produce a smaller error than the upper bound found by the algorithm on the original dataset. 3.7 Conclusions This chapter considers the problem of locating a given number of hyperplanes in order to minimize an objective function of the distances from a set of points. Each point is assigned to its closest hyperplane, thus inducing as many clusters as the number of fitting hyperplanes. The distance from each point to its corresponding fitting hyperplane can be seen as a residual and these residuals are aggregated using ordered median functions that are ordered weighted averages representing different types of utilities. Two exact approaches are presented to solve the problem. The first one is based on a compact mixed integer formulation whereas the second one is an extended set partitioning formulation with an exponential number of variables that is handled by a branch-and-price approach. To enhance the performance of this last method we have developed a generator of initial feasible solutions based on geometrical properties of the optimal solutions of the hyperplane location problem that we have also derived in this chapter, and that are used to initialize the column generation routine of this branch-and-price. We have also presented a heuristic pricing strategy that is used in combination with the exact one to speed up some pricing iterations. We report the comparison of both method to solve the problem in two different datasets on an extensive battery of computational experiments. The issue of scalability of the exact methods is also analyzed obtaining theoretical upper bounds of the error induced by some aggregated versions of the original dataset.
Chapter 4 SVM-based classification with label noise 102
104 Chapter 4. SVM-based classification with label noise In this chapter we propose novel methodologies to optimally construct SVM- based classifiers that takes into account that label noise occur in the training sample. We propose different alternatives based on solving Mixed Integer Linear and Non Linear models by incorporating decisions on relabeling some of the observations in the training dataset. The first method incorporates relabeling directly in the SVM model while a second family of methods combines clustering with classification at the same time, giving rise to a model that applies simultaneously similarity measures and SVM. Extensive computational experiments are reported based on a battery of standard datasets taken from UCI Machine Learning repository, showing the effectiveness of the proposed approaches. 4.1 Introduction Amongst the most relevant applications of classification methods there are those related with security, as in spam filtering or intrusion detection. The main difference of these applications with respect to other uses of classification approaches is that malicious adversaries can adaptively manipulate their data to mislead the outcome of an automatic analysis. For instance, spammers often modify their emails by obfuscating words which typically appear in known spam or by adding words which are likely to appear in legitimate emails. Also, as stated in Weerasinghe et al. (2019), when Machine Learning algorithms utilized in safety-critical environments are compromised by adversaries, it could even result in loss of human lives. Note that, doubting on the reliability of the labels on the target variable is usual when having suspicions about the possibility of an intentional flip amongst these labels. However, it is not by far the only case in which one must think about this possibility. Nowadays, it is commonly said that data scientists spend a large percentage of their time dealing with collecting and preprocessing data, meanwhile the remainder is used to model and extract information from databases. Mistakes converted into wrong label assignments are very likely to happen. For instance, data can be wrongly identified at the very beginning of the data collection phase, or code errors can occur when preprocessing a database, leading to a dataset with label noise. Then, one has to, not only derive a classification rule from a training sample, able to adequately classify out-of-sample data, but also to take into account that some of the labels might be incorrect. The goal of this chapter is to analyze the power of using Mathematical Programming tools for the label noise detection when constructing a SVM classifier. As pointed out in Ganapathiraju et al. (2000), amongst all the available optimization-based classifiers, SVM particularly suffer the effect of noisy labels because their reliance on support vectors and the feature interdependence assumption.
4.1. Introduction 105 Analyzing the vulnerabilities of classifiers and their robustness against attacks, to better understand how their security may be improved, has recently received growing interest from the scientific community. Bi and Zhang (2005) provided robust alternatives when the features of the training sample observations are corrupted. On the other hand, Biggio et al. (2011) proposed an algorithmic approach to handle adversarial modifications of the labels, in case the labels are independently flipped with the same probability, by correcting the kernel matrix. According to Nalepa and Kawulok (2019), three main groups of approaches for dealing with noisy datasets have been already proposed in the literature: (1) Design of algorithms which filter noisy and/or mislabeled vectors from the input data (as in Ekambaram et al. (2016), or Han and Chang (2013)); (2) Construction of robust classifiers against noisy labelling (see Duan and Wu (2016); Natarajan et al. (2017)); and (3) Use of noise models in parallel with the obtention of the classifier, which are finally coupled for a higher-quality classification (see Ganapathiraju et al. (2000); Weerasinghe et al. (2019)). Further details on the different approaches to deal with datasets containing mislabeled observations can be found in the survey in Fr´enay and Verleysen (2013). Most recent methodologies to deal with noisy datasets are sequential. Thus, loosing the optimal performance obtained by one shot methods based on Mathematical Programming approaches. For instance, in the recent method presented in Northcutt et al. (2021), based in SVM with Confident Learning (SVM-CL) approach, the authors propose a probabilistic method in three sequential phases: 1) estimate the transition matrix of class-conditional label noise, 2) filter out noisy examples, and 3) train the dataset once noisy data are removed via Co-Teaching. Analogously, in de Fran¸ca and Coelho (2015) it is proposed a novel method in which first the training sample is biclustered (see e.g., (Yang et al., 2003)) trying to capture correlation between features and observations, next the training sample is modified according to the biclusters, and then the classification is performed on the modified dataset. Furthermore, there are some globally optimal methods that have been proposed in the literature. In particular, in Bertsimas et al. (2019), the authors present different robust adaptations of classical classification methods to deal with uncertainty in labels and/or features in the training sample. In contrast to those methods that have been already proposed to deal with classification and noisy labels, our approach simultaneously construct a SVM-based classifier and re-label observations, leading to an optimal method. In addition, this approach allows one to get separating hyperplanes that would had been impossible to obtain throughout standard SVM and that report better results for many different problems. Although the method proposed in Bertsimas et al. (2019) also optimally constructs the classifier under the presence of noisy labels, it is thought to be robust
106 Chapter 4. SVM-based classification with label noise against the worse possible situation. On the contrary, our method builds the hyperplanes always on the convenience of finding good classifiers and not to be protected against the worst possible flip of labels which results in better classifiers in most scenarios. The construction of SVM-based classifiers that simultaneously relabel observations has many advantages when dealing with label noise datasets, but also when working on problems in which false positives and false negatives have different misclassifying costs. Also, in problems with unbalanced classes (as for instance in datasets on fraud with credit card transactions in which around a 99.9% of the observations are not fraudulent transactions (Maldonado et al., 2017) or in the number of claims in non-life insurances (Boucher et al., 2009)). In Figure 4.1 we illustrate this situation. One can observe in the left picture the projection on the plane of a set of observations labeled by fraudulent (red) and non fraudulent (green) transactions. Linear separators seems to be impossible to construct for this instance, but also non linear classifiers may result in overfitting. However, as shown in the right picture, if one allows a few of the labels to be changed, one can obtain better classifiers. Note that in this case, false positives are more costly than false negatives (since asking for a little more of information via text message on the phone normally solves this true negative cases). It is also important to remark that this separating hyperplane could not have been obtained through standard SVM since all the support vectors belong to the same class (green points). Figure 4.1: Original data (left) and optimal hyperplane separating re-labeled classes with our method (right). In this chapter we propose two different approaches. We present a model in which re-labeling observations depends on the errors of the SVM-based method itself searching for a compromise between the gain obtained in misclassification error and
4.2. Mathematical Programming formulations 107 margin and the penalty paid for each change of labels. On the other hand, we will also introduce two models in which re-labeled observations will come from similarity measures on the data. Our method is distribution-free so that it does not assume any distribution on the dataset and the detection of mislabeled observations and the construction of the classifier is optimal based on solving an add-hoc Mathematical Programming problem. To asses the validity of these methods we have performed a battery of computational experiments on 7 different real datasets. For these datasets we have repeated the experiments for 5 different scenarios, by randomly flipping a 0%, 20%, 30%, 40% or 50% of the labels in the original data. When comparing our method with respect to classical SVM, and with SVM-CL from (Northcutt et al., 2021), we can see that ours gets better results on noisy label datasets. 4.2 Mathematical Programming formulations 4.2.1 Preliminaries In this chapter we consider a given training sample X={(x1, y1),...,(xn, yn)} ⊆ Rp×{−1,+1}and a SVM-based separating hyperplane H={z∈Rp:ω0z+ω0= 0}, ω∈Rp,ω0∈R. We show in Figure 4.2 an example where we can see a set of points belonging to two different, blue and green, classes (left picture) and its SVM optimal solution for a given parameter c(right picture). The black line is the separating hyperplane while the other two parallel lines are delimiting the strip, S, between classes. The points that lie on these parallel lines, the boundary of the strip, are the so called support vectors, and they verify that |ω0xi+ω0|= 1. Finally, we represent in red color the magnitude of the errors induced by margin violations. Figure 4.2: Original set of points (left) and optimal SVM solution on these points (right). If we further analyze the above dataset, we can see that there are four blue observations at the very right of the dataset, and two green observations on the left
108 Chapter 4. SVM-based classification with label noise that have a strong impact when building the classifier. These observations do not allow one to construct a SVM separator of the dataset as the one we can see in Figure 4.3, since that would lead to very big misclassification errors with a very tiny margin. Figure 4.3: Not optimal solution on the SVM problem. Moreover, there are another two green observations, besides the two on the left, that are closer to the blue cloud of points than to the green one. Hence, if we could consider that these four green points and the four blue ones on the right were wrongly labeled (because of their closeness to the rest of points), we might consider a separating hyperplane with a slope like the one presented on the left of Figure 4.4 as a better classifier. However, this separating hyperplane would be impossible to obtain with the SVM model since all the support vectors belong to the same class and to avoid huge misclassification errors the model would forbid such a slope. Motivated by the above kind of configurations, we have studied different models in which a separating hyperplane is obtained not only based on the original labels but also on the possibility of relabeling some of the original observations of the training sample at a given penalty cost. We say that an observation is relabelled if one of the following assumptions occurs: yi=±1 but our model considers that yi=∓1. We will use the notation ˆyito represent the class that the model is considering for observation i. Hence, an observation is said to be relabelled if yi6= ˆyi. Following the example shown in Figures 4.2 and 4.3, we can see on the right of Figure 4.4 the solution of our model, with a separating hyperplane with the desired
4.3. Experiments 115 data. We have performed a 5-fold cross validation scheme. Thus, data have been split into 5 train-test random partitions. In each of these folds we have trained our models and we have used the other four folds for testing. Moreover, we have repeated this 5-fold cross validation 5 times for each dataset, in order to avoid beneficial starting partitions, and we report the average results obtained. For all the instances we have trained our three models and we have compared them with standard SVM and SVMCL (Northcutt et al., 2021). We have considered standard SVM as benchmark since, despite the good results provided by SVM-CL for some experiments, standard SVM provided a better performance on average among all the experiments (see Table 4.1 and Figure 4.5). The measure used to evaluate the performance of the models have been the accuracy, in percentage, on out of sample data. The parameters that appears in the different methods that we compare are validated as usual, that is, for each of the instances we perform a grid search on the cost parameters and the best result obtained in the validation sample among these parameters is the one reported. More specifically, the grids used in the experiments are the following: SVM: c∈10i:i=−5,...,5. RE-SVM: c1, c2∈10i:i=−5,...,5. 2-medians-SVM: c1, c2∈10i:i=−5,...,5,c3∈10i:i=−3,...,0. 2-means-SVM: c1, c2∈10i:i=−5,...,5,c3∈10i:i=−3,...,0. SVM-CL: Default tuning parameters (see (Northcutt et al., 2021)). The Mathematical Programming models were coded in Python 3.6, and solved using Gurobi 7.5.2 on a PC Intel Core i7-7700 processor at 2.81 GHz and 16GB of RAM. Due to the complexity of the 2-means-SVM, we have helped the solver uploading an initial feasible solution that was obtained in the 2-medians-SVM problem. We have not solved to optimality all the instances, especially those with the 2-means-SVM in which the problem becomes nonlinear, and hence we have established a time limit of 30 seconds for all the experiments. This training time has sufficed to obtain rather good classifiers. Indeed, as one can observe from the results obtained, this time limit is adequate to construct robust classifiers under noisy labels. Note that not guarantying the optimality of the solutions of our models does not necessarily imply that the classifiers are not adequate. In Table 4.1 we report the average accuracy results obtained in all the experiments for the different models and the different levels of label-noise. In such a table we have used the yellow-green color to indicate the results in which we are a 3%−5%
116 Chapter 4. SVM-based classification with label noise better than the benchmark, the green color to indicate whether we are a 5% −10% better than the benchmark, and the cyan color to highlight the results in wich we are at least a 10% above the benchmark. Also, we show in Figure 4.5 the accuracy boxplots of the 625 instances per dataset (5 partitions ×5 scenarios ×5 folds ×5 models). Regarding to the results, several conclusions can be pointed out: •Our three models perform consistently better than classical SVM when the training dataset is corrupted. Besides, the stronger the percentage of flipped labels, the bigger the difference between our models’ results and SVM’s results. In Figure 4.5 one can check how SVM model has lower tails and wider boxes than RE-SVM. •2-medians-SVM and 2-means-SVM perform better than RE-SVM for heavy attacks (40% −50% of flipped observations). In contrast, the cluster-based models require more time to be trained than RE-SVM, both because the problems are harder to solve (apart from relabeling, the distances to the centroids and the assignments observations-to-centroids are modeled) and the number of parameters that must be tuned. In Figure 4.5 one can easily check that RE-SVM has wider boxes than 2-medians-SVM and 2-means-SVM, which are explained by the behavior of these models against the attacks. •Our models have a better performance than the rest of approaches even for the original datasets in which no labels are flipped. This is due to the flexibility offered by methodologies, because some of the observations are allowed to be relabeled looking for a better classifier. The original datasets may contain outliers that contaminate the sample and so they deteriorate the classifier. This situation is automatically detected and fixed by our methods, by adequately relabeling observations. •Our methods outperform SVM-CL, which is a specialized method, designed to detect noisy labels. The rationale under these results is that SVM-CL seems to wrongly identify the right distributions of the data. These mistakes propagate to the construction of the classifier since it is built on some incomplete data. This fact also results in worse accuracies than standard SVM that works with the entire dataset without paying attention to the existence of outliers. Overall, as one may expect and it is confirmed in our computational experiments, it is better to construct the classifier without identifying incorrectly the noise labels (as SVM does) than using inadequate flips to build the classifier (as SVM-CL seems to do in the tested datasets). Obviously, the results in the experiments also show
4.3. Experiments 117 that it is rather advantageous the correct identification of the wrong labels since it improves significantly the classification rates.
118 Chapter 4. SVM-based classification with label noise Percentage of Flipped Labels Dataset Method 0% 20% 30% 40% 50% SVM-CL 84.55 82.10 71.12 58.93 49.68 Australian SVM 86.11 85.43 79.23 68.13 59.47 (690,14) RE-SVM 86.42 85.68 83.37 76.97 66.13 2-medians-SVM 86.08 85.84 84.67 78.95 69.54 2-means-SVM 85.97 85.74 82.65 77.14 67.70 SVM-CL 95.73 91.87 87.37 78.49 58.36 BreastCancer SVM 96.49 93.47 89.96 85.94 68.16 (683,9) RE-SVM 96.88 96.20 94.97 90.36 77.00 2-medians-SVM 96.63 95.31 94.46 91.10 87.31 2-means-SVM 96.96 95.93 95.39 93.11 90.01 SVM-CL 78.70 71.03 60.09 56.01 49.66 Heart SVM 82.23 76.86 69.68 63.79 56.90 (270,13) RE-SVM 82.84 78.38 73.16 68.86 61.25 2-medians-SVM 82.01 78.75 77.29 75.38 71.99 2-means-SVM 82.06 78.81 77.40 75.97 72.90 SVM-CL 78.18 65.56 59.47 55.58 49.29 Parkinson SVM 81.66 74.74 70.17 62.28 57.82 (240,40) RE-SVM 82.43 77.64 73.22 67.29 62.97 2-medians-SVM 80.32 78.62 78.12 77.51 76.28 2-means-SVM 80.47 79.22 78.78 78.20 77.03 SVM-CL 81.62 78.86 74.07 56.78 46.78 QSARbiodeg SVM 82.12 78.07 74.09 63.38 48.97 (1055,40) RE-SVM 84.53 79.61 75.00 66.42 54.58 2-medians-SVM 84.08 78.79 74.32 67.87 67.02 2-means-SVM 83.61 78.55 74.42 67.86 66.81 SVM-CL 80.94 72.79 68.54 60.53 50.69 Vertebral SVM 84.51 75.43 71.34 66.78 57.47 (310,6) RE-SVM 85.10 79.61 74.83 72.33 67.92 2-medians-SVM 85.31 82.62 80.80 78.30 76.31 2-means-SVM 86.28 84.32 81.77 79.91 76.76 SVM-CL 88.98 85.40 78.03 57.19 45.42 Wholesale SVM 90.08 85.30 79.74 72.23 57.73 (440,7) RE-SVM 90.39 88.77 85.97 80.12 69.07 2-medians-SVM 90.58 89.54 87.79 82.78 73.54 2-means-SVM 91.23 89.56 87.39 85.88 82.92 Table 4.1: Accuracy results of our computational experiments.
4.4. Conclusions 119 Figure 4.5: Accuracy Boxplots of the obtained accuracies. 4.4 Conclusions This chapter presents a methodology to construct a classification rule that at the same time incorporates the detection of label noise in the datasets. Our methodology combines the power of SVM and the features of clustering analysis to simultaneously identify wrong labels to build a separating hyperplane maximizing the margin, minimizing the misclassification errors and penalizing relabeling. The rationale is simple: observations identified as wrongly labeled will be relabelled only if the gain in margin or the decrease in misclassification error compensate the flipping. In spite of its theoretical simplicity we show the exceptional performance of our methodology in a number of databases taken from the UCI repository. These models are implemented using Mathematical Programming formulations with some integer variables (MIP). In all cases, they give rise to models that are simple and that enjoy the quality of being solvable by nowadays off-the-shelf commercial solvers (Gurobi, CPLEX, XPRESS...) Our findings are not only of theoretical interest. Its practical performance when applied to databases is remarkable. In all tested cases, our methods are superior to the considered benchmark that in our case is standard SVM. Thus, they are directly applicable to datasets in which flipped labels are suspected, resulting in robust classifiers to noisy labels.
Chapter 5 Robust Optimal Classification Trees under Noisy Labels 122
124 Chapter 5. Robust Optimal Classification Trees under Noisy Labels In this chapter we propose a novel methodology to construct Optimal Classification Trees that takes into account that noisy labels may occur in the training sample. The motivation of this new methodology is based on the superaditive effect of combining together margin based classifiers and outlier detection techniques. Our approach rests on two main elements: (1) the splitting rules for the classification trees are designed to maximize the separation margin between classes applying the paradigm of SVM; and (2) some of the labels of the training sample are allowed to be changed during the construction of the tree trying to detect the label noise. Both features are considered and integrated together to design the resulting Optimal Classification Tree. We present a Mixed Integer Non Linear Programming formulation for the problem, suitable to be solved using any of the available off-the-shelf solvers. The model is analyzed and tested on a battery of standard datasets taken from UCI Machine Learning repository, showing the effectiveness of our approach. Our computational results show that in most cases the new methodology outperforms both in accuracy and AUC the results of the benchmarks provided by OCT and OCT-H. 5.1 Introduction Interpretability is a crucial requisite demanded to Machine Learning methods provoked by the tremendous amount of methodologies that have arised in the last decade (Du et al., 2019). It is expected that the model that results when applying a Machine Learning methodology using a training sample, apart from being able to adequately predict the behaviour of out-of-sample observations, can be interpreted. Different tools have been applied to derive interpretable Machine Learning methods. One of the most popular strategies to simplify the obtained models is feature selection, in which a reduced set of attributes is to be selected without loosing quality in the predictions. Reducing the number of parameters to analyze, the models can be easier to understand, yielding higher descriptive accuracy. One could also consider models that can be modulated, in the sense that a great proportion of its prediction-making process can be interpreted independently. This is the case of generalized linear models (Hastie and Tibshirani, 2017). Other methods incorporate interpretability as a synonym of being able to be reproduced by humans in its entire construction(Letham et al., 2015). This is the case of Decision Trees with small depth which can be visualized and interpreted easily by users even not familiar with the tools behind their construction. CART is the most popular Decision Trees method. In CART, one constructs the decision rule based on a hierarchical relation amongst a set of nodes which is used to define paths that lead observations from the root node (highest node in the hierar-
5.3. Optimal Classification Trees with SVM splits and Relabeling (OCTSVM) 131 two methods are obtaining a perfect classification on the training sample, Figure 5.2 shows that OCTSVM provides a more balanced solution than OCT-H since it has wider margins between both classes, what could be translated into a higher accuracy for out-of-sample observations. Figure 5.2: Optimal solutions for OCT-H with D= 2 (left) and OCTSVM with D= 1 (right). In order to formulate the OCTSVM as a MINLP, we will start describing the rationale of its objective function that must account for the margins induced by the splits at all nodes, the classification errors, the penalty paid for relabelling observations and the cost-complexity of the final classification tree. To formulate the above concepts we need different sets of variables. First of all, we consider continuous variables: ωt∈Rp, ωt0∈R, t = 1, . . . , T, which represent the coefficients and the intercept of the hyperplane performing the split at node t. Taking into account that the margin of the hyperplane Ht={z:ω0 tz+ωt0= 0}is given by 2 ||ωt|| , maximizing the minimum margin between classes induced by the splits can be done introducing an auxiliary variable ι∈R(that will be minimized in the objective function) which is enforced by the following constraints: 1 2||ωt||2≤ι∀t= 1, . . . , T. (5.1) Once the maximization of the margin is set, we have to model the minimization of the errors at the nodes, whereas at the same time we minimize the number of relabelled observations. These two tasks are accomplished by the variables dit ∈R, i = 1, . . . , n, t = 1, . . . , T , that account for the misclassification error of observation iat node t, and ξit ∈ {0,1}binary variables modelling whether observation iis relabeled or not at node t. If c1and c2are the unit costs of misclassification and relabelling, respectively, our goal is achieved adding to the objective function
132 Chapter 5. Robust Optimal Classification Trees under Noisy Labels the following two terms: c1 n X i=1 T X t=1 dit +c2 n X i=1 T X t=1 ξit. The correct meaning of these two sets of variables must be enforced by some families of constraints that we describe next. Nevertheless, for the sake of readability before describing those constraints modeling these dit and ξit, we must introduce another family of variables the βit ∈Rp, βi0∈R, i = 1, . . . , n, t = 1, . . . , T, which are continuous variables equal to the coefficients of the separating hyperplane at node t when observation iis relabelled, and equal to zero otherwise. In addition, we consider binary variables zit ∈ {0,1}needed to control whether observation ibelongs to node tof the tree. Now, putting all these elements together, as it is done in RE-SVM, we can properly define the splits and their errors at each node of the tree using the following constraints: yi(ω0 txi+ωt0)−2yi(β0 txi+βit0)≥1−dit −M(1 −zit), ∀i= 1, . . . , n, t= 1, . . . , T, (5.2) βitj =ξitωtj,∀i= 1, . . . , n, t = 1, . . . , T, j = 0, . . . , p. (5.3) Constraints (5.3), which can be easily linearized, are used to define the βit variables: they equal the ωtvariables when the observation iis relabelled (ξi= 1), and are equal to zero otherwise (ξi= 0). On the other hand, constraints (5.2) control the relabelling at each node of the tree. If an observation iis in node t(zit = 1), and ξi= 0, we obtain the standard SVM constraints for the separating hyperplane. Nevertheless, if ξi= 1, then the separating constraints are applied for the observation ias if its class were the opposite to its actual one, i.e., as if observation iis relabeled. Moreover, since Mis a big enough constant, these constraints do not have any kind of impact in the error variables of node tif observation idoes not belong to this node (zit = 0). On the other hand, there are still some details left that must be imposed to make the model to work as required. In the decision tree, observations start at the root node and they advance descending through the levels of the tree until they reach a leaf or a pruned node. Hence, we have to guarantee that observations must belong to one, and only one, node per level. By means of the zit variables, this can be easily done by the usual assignment constraints applied in each level, u∈U, of the tree: X t∈u zit = 1,∀i= 1, . . . , n, u ∈U. (5.4)
5.3. Optimal Classification Trees with SVM splits and Relabeling (OCTSVM) 133 Moreover, for consistency in the relation between a node and its ancestor, it is clear that if observation iis in node t(zit = 1), then, observation imust be also in the parent of node t(zip(t)= 1), with the only exception of the root node. Besides, if observation iis not in node t(zit = 0), then ican not be in its successors, and this is modeled by adding the following constraints to the problem: zit ≤zip(t),∀i= 1, . . . , n, t = 2, . . . , T. (5.5) So far, the OCTSVM model has everything it needs to properly perform the splits by following the RE-SVM rationale described in the previous chapter, taking into consideration the tree complexity, and maintaining the hierarchical relationship amongst nodes. The last element that we need to take care of, to assure the correct performance of the whole model, is to define how observations follow their paths inside the tree. We get from constraints (5.5) that observations move from parent to children (nodes), but every non terminal node has a left and a right child node, and we need to establish how observations take the left or the right branch. Since the splits are made by the separating hyperplane, we force observations that lie on the positive half space of a hyperplane to follow the right branch of the parent node, and observations that lie on the negative one to take the left branch. This behavior is modeled with the binary variables θit ∈ {0,1}, that are used to identify whether observation ilies in the positive half space of the separating hyperplane at node t, θit = 1, or if observation ilies on the negative half space, θit = 0. By considering M a big enough constant, the correct behavior of the path followed by the observations is enforced by the following constraints: ω0 txi+ωt0≥ −M(1 −θit),∀i= 1, . . . , n, t = 1, . . . , T, (5.6) ω0 txi+ωt0≤Mθit,∀i= 1, . . . , n, t = 1, . . . , T. (5.7) Hence, by making use of these θit variables, and distinguishing between nodes that come from left splits, τbl (nodes indexed by even numbers), and right splits, τbr (nodes indexed by odd numbers), we control that the observations follow the paths through the branches in the way we described above throughout the following constraints: zip(t)−zit ≤θip(t),∀i= 1, . . . , n, t ∈τbl,(5.8) zip(t)−zit ≤1−θip(t),∀i= 1, . . . , n, t ∈τbr.(5.9) According to constraints (5.8), if an observation iis on the parent node of an even node t(zip(t)= 1), and ilies on the negative half space of the hyperplane defining the split on p(t) (θip(t)= 0), then zit is forced to be 1. Hence, θip(t)= 0 implies that
134 Chapter 5. Robust Optimal Classification Trees under Noisy Labels observation itakes the left branch to the child node t∈τbl. Moreover, we can see that this constraint is consistent since if zip(t)= 1, but observation iis not in the left child node, zit = 0, t∈τbl, then θip(t)equals 1, what means that observation ilies on the positive half space of the hyperplane of p(t). On the other hand, constraints (5.9) are similar but for the right child nodes, τbr. If an observation iis in the parent node of an odd node t∈τbr, and ilies on the positive half space of the hyperplane of p(t) (θip(t)= 1), then, zit = 1 what means that observation ihas to be on node t. The final term to be included in the objective function of the problem is the complexity of the resulting tree. Following the approach in OCT and OCT-H, we consider binary variables δt∈ {0,1}, t = 1, . . . , T, that control whether a separating split is applied at node t. Thus, to control the tree complexity resulting of the process, we minimize the sum of these variables multiplied by a cost penalty c3. Gathering all the components together, the objective function to be minimized in our problem results in ι+c1 n X i=1 T X t=1 dit +c2 n X i=1 T X t=1 ξit +c3 T X t=1 δt. According to this, it is important to make a distinction between the methods that use, or the ones that do not use, SVM based splits. When a model does not use SVM based splits, taking into account that the binary variable δt= 1 implies that a hyperplane is being used to split the points in node t, complexity can be easily regulated by just imposing kωtk2≤Mδtin all the branch nodes. Nevertheless, when using a SVM based method such as the one we are presenting here, the previous constraint would be in conflict with constraints (5.2). This is due to the fact that δt= 0 would imply kωtk2= 0, and under this scenario, observations in this node would have to pay for the margin violation error dit = 1, even though these errors are not justified since points are not being separated at this node. To overcome this issue, we need to add some other constraints to the model, but before getting into the mathematical formulation, we would like to emphasize the different aspects that explain the difficulty of the problem. The objective function accounts for the complexity cost of node t(δt= 1) when a hyperplane is actually built so as to split the points in node t. In case we do not pay the complexity cost in node t(δt= 0), it does not necessarily mean that a hyperplane is not built at this node, it could simply turn out to be a hyperplane that leaves all the observations at one of the half spaces it creates, i.e., a hyperplane that is not splitting the points. These non splitting hyperplanes do not affect the training set at all, and more importantly, they would not affect out of sample observations since predictions will occur in leaf nodes or in the first branch node in which δt= 0. Going back to the formulation, the first
5.3. Optimal Classification Trees with SVM splits and Relabeling (OCTSVM) 135 step now is to introduce some binary variables, hit ∈ {0,1}, i = 1, . . . , n, t ∈τb, that will be relaxed afterwards, defined by hit =zitθit. These variables tell us which observations belong to node tand at the same time lie on the positive half space of the split created in t. It is important to define these variables because they are used to measure whether the hyperplane is splitting the points in node tor not. This is going to be done by means of some new binary variables vt∈ {0,1}, t ∈τb, that will be equal to zero in case all the points in node tbelong to the positive half space of the hyperplane built in this node, Ht, and one otherwise. For the definition of the δtvariables to make sense, the following constraints must simultaneously be considered: hit =zitθit,∀i= 1, . . . , n, t = 1, . . . , T, (5.10) n X i=1 (zit −hit)≤Mvt,∀t= 1, . . . , T, (5.11) n X i=1 hit −M(1 −vt)≤Mδt,∀t= 1, . . . , T, (5.12) where Mis a big enough constant. The reader should observe that (5.10) is the definition of the hit variables as the product of zit by θit. These constraints can be easily linearized by the usual tricks. On the other hand, taking into account that vtvariables would tend to be zero (by the effects on constraints (5.12)), constraints (5.11) stand for the definition of the vtvariables, since these variables could be equal to zero just in case all the hit are equal to one, what means that all the observations belong to the positive half space of Ht. Finally, we obtain through constraints (5.12) the definition of the δtvariables. Since these variables are being minimized in the objective function, they would try to be equal to zero. Nevertheless, this can just happen if vt= 0, what means that all the observations belong to the positive half space of Ht, or if n X i=1 hit = 0 in case vt= 1, what means that all the observations belong to the negative half space of Ht. Hence, δtis equal to one if and only if the points in node tare being actually separated by Ht. Finally, another point that it is important to remark about the δtvariables is that once a non-leaf node does not split (that is, the corresponding branch is pruned at this node), the successors of node tcan not make splits either to maintain the correct hierarchical structure of the tree. Recalling that p(t) is the parent node of node t, we can guarantee this relationship throughout the following constraints δt≤δp(t),∀t= 2, . . . , T. (5.13) Some of the constraints discussed above require big-Mconstants. The value of
136 Chapter 5. Robust Optimal Classification Trees under Noisy Labels these constants can be adjusted once we know the dimension (n, p) of the dataset involved in the optimization problem. Recall that xi∈[0,1]p, therefore it is well known that the maximum distance between two points in such a domain is √p. Hence, the maximum distance from a point to a hyperplane is also √p. According to this, the big-Mconstant in constraints (5.2), (5.5) and (5.6) is actually √p. On the other hand, the hit variables are upper bounded by one, and therefore n X i=1 hit is upper bounded by n. As a result, the big-Mconstant involved in constraints (5.10) and (5.11) is nas well. Gathering all the constraints together, and substituting the generic big-M constants by those estimated in our discussion above, the OCTSVM is obtained by solving the following MINLP: min ι+c1 n X i=1 T X t=1 dit +c2 n X i=1 T X t=1 ξit +c3 T X t=1 δt(OCTSVM) s.t. 1 2||ωt||2≤ι, ∀t= 1, . . . , T, yi(ω0 txi+ωt0)−2yi(β0 txi+βit0)≥1−dit −√p(1 −zit),∀i= 1, . . . , n, t = 1, . . . , T, βitj =ξitωtj,∀i= 1, . . . , n, t = 1, . . . , T, j = 0, . . . , p, (5.14) X t∈u zit = 1,∀i= 1, . . . , n, u ∈U, zit ≤zip(t),∀i= 1, . . . , n, t = 2, . . . , T, ω0 txi+ωt0≥ −√p(1 −θit),∀i= 1, . . . , n, t = 1, . . . , T, ω0 txi+ωt0≤√pθit,∀i= 1, . . . , n, t = 1, . . . , T, zip(t)−zit ≤θip(t),∀i= 1, . . . , n, t ∈τbl, zip(t)−zit ≤1−θip(t),∀i= 1, . . . , n, t ∈τbr, hit =zitθit,∀i= 1, . . . , n, t = 1, . . . , T, (5.15) n X i=1 (zit −hit)≤nvt,∀t= 1, . . . , T, n X i=1 hit −n(1 −vt)≤nδt,∀t= 1, . . . , T, δt≤δp(t),∀t= 2, . . . , T, dit ∈R+, βit ∈Rp, βit0∈R, ξit, zit, θit, hit ∈ {0,1},∀i= 1, . . . , n, t = 1, . . . , T, ωt∈Rp, ωt0∈R, δt, vt∈ {0,1},∀t= 1, . . . , T, ι∈R. The reader may note that the above formulation has two families of bilinear
5.4. Experiments 137 constraints, namely (5.14) and (5.15), that can be easily linearized giving rise to a mixed integer second order cone formulation. 5.4 Experiments In this section we present the results of our computational experiments. Five different classification tree-based methods are compared, CART, OCT, OCT-H, OCT+SVM and OCTSVM, on nine popular real-life datasets from UCI Machine Learning Repository Lichman et al. (2013). Notice that OCT+SVM is a modification of our OCTSVM in which the ξit variables are fixed to zero, i.e., no relabeling is allowed in the model. This method is included so as to assess the isolated effect of using margin-based splits, without relabeling, within the classification trees. The considered datasets together with their names and dimensions (n: number of observations, p: number of features) are reported in the three tables of this section. Our computational experiments focus on the analysis of the accuracy and AUC (area under the ROC curve) of the different classification tree-based methods. This analysis is based in four different experiments for each one of the nine considered dataset. Our goal is to analyze the usefulness of the different methods for classifying data affected by label noise. Therefore, in our experiments we use, apart from the original datasets, three different modifications where in each one of them a percentage of the labels in the sample are randomly flipped. The percentages of flipped labels range in {0%,20%,30%,40%}, where 0% indicates that the original training data set is used to construct the tree. We perform a 4-fold cross-validation scheme, i.e., datasets are split into four random train-test partitions. One of the folds is used for training the model while the rest are used for testing. When testing the classification rules, we compute the accuracy, in percentage, on out of sample data, as well as the AUC, which provides an aggregate measure of performance across all possible classification thresholds. One way of interpreting AUC is as the probability that the model ranks a random positive example more highly than a random negative example. The CART method was coded in Rusing the rpart library. OCT, OCT- H, OCT+SVM and OCTSVM Mathematical Programming models were coded in Python and solved using Gurobi 8.1.1 on a PC Intel Core i7-7700 processor at 2.81 GHz and 16GB of RAM. A time limit of 30 seconds was set for training. Not all the problems were solved to optimality, nevertheless, this time limit was enough in order to obtain good classifiers. The average gap of the experiments is reported in table 5.3. The calibration of the parameters of the different optimization-based models compared in these experiments was set as follows:
138 Chapter 5. Robust Optimal Classification Trees under Noisy Labels •For OCTSVM we used a grid on 10i:i=−5,...,5for the constants c1and c2, and a smaller one 10i:i=−2,...,2for c3. For OCT+SVM c2is not used since all the ξit variables are set to zero. •For OCT, in order to check every possible optimal subtree of a maximal tree for a given depth D, we did the predictions in two steps. We first set up a grid on the parameter c1=1,...,2D−1and solved a slightly different OCT model in which the objective function did not take into account the complexity term, since complexity was already considered in the model by adding the following constraint X t∈τB δt≤c1. The resulting solutions that were optimal for a certain c1of the original OCT problem were afterwards computed. This methodology was not used in OCT-H since the grid of the modified problem should have been extended to c1=1, . . . , p(2D−1), therefore we used the same grid that we used for c1in OCTSVM directly into OCT-H formulation. On the other hand, the minimum number of observations per leaf was set to a 5% of the training sample size. •CART trees were post-processed to satisfy any depth constraint as it was done with OCT. The minimum number of observations per leaf was the same 5% of the training sample size that we used for OCT and OCT-H. Last to mention, the depth, D, considered in these experiments was set equal to three for CART, OCT and OCT-H, whereas for OCTSVM and OCT+SVM we fixed depth equal to two, creating consequently trees with 3 levels, to set a fair comparison amongst the different methods. For each dataset, we replicate each experiment four times. In Table 5.1 we show the average (±standard deviation) accuracies for each one of the methods and datasets. In addition, Table 5.2 shows the average (±standard deviation) AUC results. The best results are highlighted (boldfaced). The first column stands for the percentage of flipped labels (FL) of the training sets. On the other hand, the last column shows the mean difference (Diff) between OCTSVM and the best result amongst OCT and OCT-H. Finally, in Table 5.3 we report the average MINLP gaps obtained by the models at the end of the training time limit. 5.5 Conclusions One of the classification methods that has experienced a more in depth transformation in the last years is classification trees. Since the pioneer contribution by Breiman et al. Breiman et al. (1984), where CART was proposed, this technique has included tools from mathematical optimization giving rise to the so called OCT
5.5. Conclusions 139 % FL CART OCT OCT-H OCT+SVM OCTSVM Diff 0 84.91 ±1.31 85.48 ±0.99 85.16 ±1.28 85.62 ±1.08 86.16 ±0.68 0.68 ±0.76 Australian 20 83.35 ±5.30 85.31 ±1.12 80.05 ±4.57 81.93 ±7.83 85.25 ±1.34 -0.06 ±1.03 (690,14) 30 70.50 ±11.89 79.27 ±7.55 72.01 ±4.25 75.24 ±10.46 82.57 ±7.46 3.29 ±11.96 40 54.34 ±7.94 70.36 ±9.56 64.63 ±4.25 64.38 ±7.93 77.92 ±8.19 7.56 ±11.71 0 91.64 ±2.09 87.94 ±1.54 98.60 ±0.47 92.68 ±4.86 96.06 ±3.19 -2.54 ±3.32 Banknote 20 88.00 ±2.63 85.99 ±1.91 70.01 ±9.18 66.46 ±10.01 89.97 ±4.71 3.97 ±4.63 (1372,5) 30 85.79 ±3.38 85.96 ±2.27 63.18 ±9.91 67.36 ±12.85 84.03 ±12.46 -1.93 ±12.72 40 69.13 ±10.33 77.73 ±9.25 59.15 ±8.53 58.45 ±12.66 87.16 ±12.46 9.42 ±8.67 0 91.81 ±1.51 93.71 ±0.94 94.21 ±1.06 96.15 ±0.94 95.09 ±5.27 2.44 ±1.07 BreastCancer 20 90.38 ±2.63 92.88 ±1.72 89.45 ±4.24 82.96 ±5.29 91.85 ±3.07 -1.03 ±3.08 (683,9) 30 87.52 ±3.38 92.15 ±1.84 84.41 ±5.41 79.16 ±4.94 89.90 ±2.75 -2.24 ±3.85 40 73.95 ±12.53 90.44 ±2.68 75.36 ±6.81 74.58 ±5.15 87.28 ±4.71 -3.15 ±5.26 0 72.34 ±3.57 75.02 ±1.82 78.95 ±2.29 82.31 ±2.05 83.36 ±1.41 4.41 ±1.88 Heart 20 70.53 ±6.16 74.10 ±3.57 74.66 ±4.24 73.95 ±5.98 80.13 ±4.39 5.46 ±4.92 (270,13) 30 71.87 ±3.72 71.87 ±6.23 71.19 ±5.61 73.7 ±4.69 77.25 ±6.31 5.37 ±7.24 40 71.72 ±3.63 64.18 ±5.32 65.18 ±5.04 66.82 ±5.15 74.26 ±3.39 9.07 ±4.94 0 81.09 ±4.73 85.34 ±2.36 84.77 ±1.97 83.15 ±4.46 85.55 ±2.65 0.20 ±2.70 Ionosphere 20 74.46 ±6.77 78.74 ±3.14 74.88 ±5.78 77.06 ±3.63 78.05 ±4.65 -0.68 ±5.43 (351,34) 30 65.80 ±6.76 74.30 ±5.70 73.61 ±4.47 74.84 ±4.46 76.83 ±5.21 2.53 ±7.27 40 60.48 ±7.39 70.89 ±5.30 66.06 ±6.69 71.73 ±4.64 75.68 ±2.62 4.79 ±5.77 0 59.19 ±3.04 60.30 ±2.90 61.24 ±2.40 60.48 ±2.59 62.85 ±1.72 1.61 ±2.90 MONK’s 20 55.40 ±7.01 59.59 ±3.42 58.96 ±4.93 59.49 ±1.77 62.40 ±2.15 2.80 ±3.79 (415,17) 30 55.10 ±4.88 57.90 ±5.00 59.21 ±2.91 58.73 ±0.97 60.57 ±4.67 1.36 ±4.93 40 50.30 ±5.99 56.86 ±3.04 56.86 ±3.61 58.84 ±2.11 61.91 ±2.30 5.05 ±4.41 0 69.19 ±8.95 76.09 ±4.82 74.82 ±3.18 80.66 ±2.26 81.67 ±4.26 5.57 ±4.26 Parkinson 20 61.95 ±12.39 73.95 ±4.34 66.40 ±5.20 74.96 ±4.69 77.29 ±3.88 3.33 ±5.27 (240,40) 30 57.82 ±10.61 71.05 ±7.03 62.63 ±6.46 70.27 ±5.42 75.57 ±3.70 4.52 ±6.56 40 51.76 ±8.18 67.95 ±7.44 59.39 ±3.98 62.38 ±5.51 69.59 ±5.10 1.64 ±9.46 0 62.11 ±10.86 69.17 ±4.82 71.53 ±4.69 74.15 ±3.70 76.86 ±3.86 5.32 ±6.16 Sonar 20 58.06 ±9.03 64.41 ±5.65 67.41 ±4.07 75.21 ±5.59 71.62 ±4.29 4.20 ±4.84 (208,60) 30 53.35 ±6.40 60.62 ±4.78 65.11 ±5.76 67.36 ±7.27 70.07 ±4.36 4.96 ±5.38 40 50.04 ±6.76 58.79 ±5.18 60.24 ±4.15 62.44 ±4.13 64.95 ±11.88 4.71 ±10.68 0 89.61 ±2.22 90.35 ±1.28 90.52 ±1.27 87.54 ±3.80 88.59 ±1.82 -1.93 ±2.63 Wholesale 20 81.22 ±10.12 87.80 ±3.56 85.14 ±2.95 75.21 ±5.59 81.63 ±4.67 -6.17 ±5.67 (440,7) 30 77.40 ±12.03 83.99 ±5.15 78.01 ±7.01 74.89 ±6.35 78.95 ±5.06 -5.03 ±8.30 40 61.46 ±14.41 73.14 ±16.11 72.31 ±5.70 72.01 ±10.52 76.06 ±4.94 3.76 ±16.90 Table 5.1: Average accuracies (±standard deviations) obtained in our computational experiments.
140 Chapter 5. Robust Optimal Classification Trees under Noisy Labels % FL CART OCT OCT-H OCT+SVM OCTSVM Diff 0 85.31 ±1.72 86.10 ±0.93 85.18 ±1.52 85.92 ±1.16 86.30 ±0.78 0.20 ±0.70 Australian 20 83.78 ±5.81 85.88 ±1.23 80.33 ±5.19 81.87 ±9.42 85.52 ±8.77 -0.35 ±1.15 (690,14) 30 70.26 ±12.49 79.44 ±7.70 71.87 ±4.32 74.89 ±11.96 82.70 ±8.77 3.26 ±13.08 40 53.98 ±6.12 70.04 ±10.26 64.91 ±5.58 63.07 ±9.80 77.39 ±9.39 7.35 ±12.53 0 91.37 ±2.15 87.52 ±1.47 98.56 ±0.43 90.84 ±7.11 96.01 ±3.23 -2.55 ±3.34 Banknote 20 87.80 ±2.69 85.87 ±1.76 66.91 ±10.95 63.38 ±12.15 89.47 ±5.26 3.60 ±5.17 (1372,5) 30 85.33 ±3.90 85.54 ±2.40 60.13 ±11.04 64.97 ±13.83 84.17 ±11.63 -1.37 ±11.84 40 67.54 ±12.06 76.64 ±10.95 56.05 ±8.14 55.68 ±12.33 85.85 ±10.02 9.20 ±9.98 0 91.36 ±1.92 93.10 ±1.33 93.60 ±1.54 95.73 ±1.27 94.20 ±7.98 2.62 ±1.41 BreastCancer 20 89.62 ±2.86 91.95 ±2.46 88.52 ±5.76 79.58 ±8.69 90.32 ±3.40 -1.09 ±4.14 (683,9) 30 85.28 ±8.18 90.82 ±2.46 82.96 ±8.50 74.00 ±7.94 87.06 ±4.37 -3.75 ±6.23 40 67.61 ±15.04 88.90 ±3.96 69.71 ±10.28 67.46 ±8.33 85.48 ±7.57 -3.42 ±9.10 0 71.91 ±2.89 74.49 ±1.97 78.62 ±2.19 82.17 ±1.53 82.87 ±1.21 4.24 ±1.96 Heart 20 69.88 ±6.09 73.31 ±3.42 74.24 ±3.91 73.04 ±6.08 79.40 ±4.29 5.15 ±4.43 (270,13) 30 71.47 ±3.04 71.16 ±6.40 70.14 ±5.82 72.75 ±5.15 73.97 ±6.87 4.73 ±7.89 40 71.45 ±3.01 63.91 ±5.34 64.07 ±6.14 66.29 ±5.85 73.97 ±3.08 9.89 ±6.29 0 78.29 ±6.21 81.90 ±3.84 81.46 ±2.04 77.49 ±5.66 82.38 ±3.71 0.48 ±4.26 Ionosphere 20 69.22 ±9.11 74.39 ±3.87 70.69 ±6.64 73.18 ±4.46 72.64 ±7.39 -1.70 ±7.48 (351,34) 30 60.11 ±8.73 70.20 ±6.56 69.09 ±4.62 70.62 ±6.42 71.37 ±5.07 1.16 ±7.39 40 56.09 ±7.09 68.23 ±5.56 60.00 ±7.84 64.32 ±6.35 69.44 ±3.51 1.21 ±6.80 0 58.49 ±4.06 59.83 ±3.21 59.84 ±3.85 59.90 ±2.78 60.92 ±2.46 1.08 ±3.16 MONK’s 20 54.93 ±6.45 59.06 ±3.96 58.29 ±3.85 57.51±2.83 60.68 ±2.36 1.62 ±4.50 (415,17) 30 53.56 ±5.26 57.25 ±5.98 57.85 ±2.95 55.86 ±2.91 59.58 ±3.09 1.73 ±3.39 40 49.51 ±4.56 56.08 ±3.41 55.67 ±4.11 57.33 ±2.57 59.79 ±2.96 4.12 ±3.81 0 69.58 ±8.13 76.15 ±1.77 74.95 ±3.05 80.77 ±2.19 81.81 ±4.07 5.66 ±4.14 Parkinson 20 62.69 ±11.63 74.08 ±4.34 66.48 ±5.27 75.13 ±4.72 77.37 ±3.47 3.29 ±5.34 (240,40) 30 58.15 ±10.43 70.89 ±7.57 62.57 ±6.71 70.18 ±5.54 75.87 ±3.47 4.97 ±7.28 40 51.94 ±8.09 67.89 ±7.61 59.36 ±3.90 62.27 ±5.98 69.32 ±5.30 1.42 ±9.52 0 62.11 ±10.86 68.68 ±4.65 71.19 ±4.77 73.64 ±4.25 76.55 ±3.68 5.35 ±6.34 Sonar 20 58.06 ±9.03 63.75 ±6.20 67.46 ±5.67 68.81 ±4.86 71.61 ±4.41 4.15 ±4.85 (208,60) 30 53.35 ±6.40 59.88 ±5.54 65.20 ±5.67 66.65 ±7.78 69.64 ±4.42 4.43 ±5.10 40 50.04 ±6.76 58.00 ±5.25 59.82 ±4.98 62.05 ±4.40 64.99 ±11.18 5.17 ±10.17 0 89.61 ±2.22 88.82 ±2.47 88.94 ±3.02 86.25 ±3.58 87.36 ±2.98 -1.57 ±4.25 Wholesale 20 81.22 ±10.12 85.80 ±5.37 82.41 ±5.53 67.10 ±12.64 77.12 ±7.94 -8.67 ±7.50 (440,7) 30 77.40 ±12.03 81.80 ±8.12 74.57 ±10.13 67.99 ±12.29 72.40 ±9.12 -9.39 ±13.03 40 61.46 ±14.41 72.30 ±15.57 64.72 ±12.85 59.88 ±10.52 69.66 ±10.35 -2.64 ±9.89 Table 5.2: Average AUC (±standard deviations) obtained in our computational experiments.