Alternating local search based VNS for linear classification

Alternating local search based VNS for linear classification Frank Plastria, Steven De Bruyne MOSI, Vrije Universiteit Brussel, Belgium. {Frank.Plastria,Steven.De.Bruyne}@vub.ac.be Emilio Carrizosa Universidad de Sevilla, Spain. [email protected] June 29, 2007 Abstract We consider the linear classification method consisting of separating two sets of points in d-space by a hyperplane. We wish to determine the hyperplane which minimises the sum of distances from all misclassified points to the hyperplane. To this end two local descent methods are developed, one grid-based and one optimisation-theory based, and are embedded in several ways into a VNS metaheuristic scheme. Computational results show these approaches to be complementary, leading to a single hybrid VNS strategy which combines both approaches to exploit the strong points of each. Extensive computational tests show that the resulting method performs well. Keywords. Data Mining, Classification, Linear Classification, Heuristic Minimisation, Normdistance, Variable Neighbourhood Search, Variable Neighborhood Search, VNS, Local Search, Grid Search, Cell Search. 1 Introduction The goal of classification is to find a simple rule to classify objects into one of given classes. The simplest rules are linear rules, and in this paper we study such rules for the separation of two classes of numerical data A, B ⊂Rd. Several criteria may be applied, e.g. the popular minimisation of the number of misclassifieds. Here we rather explore the criterion proposed by Mangasarian [9]: minimise the (possibly weighted) sum of misclassification distances, i.e. the distance to the separating hyperplane for each misclassified point of Aand B. This is a nonconvex criterion, which for separable Aand B is evidently optimized by any separating hyperplane with objective value 0, but for non separable classes turns out to admit very many local optima. The aim of this work is to develop an efficient algorithm to solve the resulting global optimization problem, so in all what follows we assume implicitly that Aand Bcannot be linearly separated. Mangasarian [9] has described an exact solution algorithm when distances are measured by the L1-norm, based on solving two LP-subproblems per dimension, but this approach cannot be generalized to other norms, in particular to the Euclidean norm. For this latter case, exact solution approaches of branch-and-cut type were developed by Audet et al. [1]; see also Karam [7] for the L∞-norm. To the best of our knowledge the only further algorithmic work on this problem is the heuristic approach of Karam et al [8], who use a Variable Neighbourhood Search (VNS, see [10, 5]) to solve problems with any Lp-norm. 1 Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 2 In this paper we report on an independent study of several heuristic approaches of VNS-type to solve the problem under any norm. First we discuss a relatively simple adaptive grid-based VNS approach, which turns out to be very quick, but does not behave too well in all cases. Then we develop a more complex local optimisation method based on two theoretical necessary optimality conditions, which yields local optima from any starting point, but, when built into a VNS framework, is much more time-consuming with somewhat disappointing results. Finally it is by combining these two search strategies into a single framework that we obtain a quite stable method which yields high quality results in acceptable times. 2 Problem formulation Let there be given two finite datasets A, B ⊂Rd. Together these data form the training set of the rule, and their number will be denoted by pdef =|A∪B|. We work with halfspaces and hyperplanes in Rd. These are defined by some pair σdef = (u, β)∈ Rd\{0}×Rby H#(σ) = H#(u, β)def ={x∈Rd| h u;xi#b} where # ∈ {≤, <, =,≥, >}, and hu;xidenotes the scalar product. Note that these sets all remain the same when σis multiplied by any strictly positive constant, but not when its sign is inversed. Any σ= (u, β) defines the following linear classification rule: For a new object x∈Rd: If hu;xi> b then we classify xin A, If hu;xi< b then we classify xin B, If hu;xi=bthen we consider xas non classified. This means in fact that we use the halfspaces H>(σ) and H<(σ) to discriminate elements from A, considered to have to be Above the hyperplane H=(σ), i.e. in H>(σ), as opposed to the elements of B, that should be Below, i.e. in H<(σ). Therefore we say that we separate by a halfspace (or oriented hyperplane), and the sign of σis thus part of the information. We will also speak of the ‘rule’ or ‘halfspace σ’ by an abuse of terminology. In order to derive a good rule σwe must be able to differentiate between well and wrongly classified points by σ. Therefore we also need to define similarly following sets for any subset C⊂Rd C#(σ) = C∩H#(σ) The set A>(σ) (resp. B<(σ)) is thus the set of correctly classified points from A(resp. B), the set A<(σ) (resp. B>(σ)) is the set of misclassified points from A(resp. B), and the set A=(σ) (resp. B=(σ)) contains the non-classified points from A(resp. B). According to the criterion proposed in [9], we say that a halfspace σ∗is an optimal separating halfspace, if it minimizes the sum of distances from all misclassified points A<(σ∗)∪B>(σ∗) to the boundary hyperplane H=(σ∗) . Thus, determining an optimal rule means solving the following optimization problem σ∗∈arg min{f(σ)|σ∈Rd\{0}×R} where f(σ)def =X a∈A<(σ) d(a, H=(σ)) + X b∈B>(σ) d(b, H=(σ)) (1) We denote the set of all halfspaces in Rdby H. As mentioned above, for any λ > 0 we have H#(λσ) = H#(σ), so the function H<:Rd\{0}×R→ H :σ= (u, β)7→ H<(σ) Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 3 is surjective, with inverse image of each halfspace a ray R+ 0(u, β) for some u6= 0 in Rd. In particular each such ray contains a single (u, β) with kuk= 1. We will assume throughout that distance is measured by a given norm γ. It is known that the distance of a point a∈Rdto a hyperplane H(u, β) (u6= 0) is then calculated as (see e.g. [13]) dγ(a, H(u, β)) = |h u;ai−β| γ◦(u)(2) where γ◦is the dual norm of γ. This expression allows to rewrite the objective function at σ= (u, β) as f(u, β) = 1 γ◦(u)  X a∈A<(σ) (β−h u;ai) + X b∈B>(σ) (hu;bi−β) (3) and it can be seen there will always be an optimal solution (u, β) for which f(u, β) = X a∈A<(σ) (β−h u;ai) + X b∈B>(σ) (hu;bi−β) by choosing γ◦(u) = 1. 3 Local descent methods In this section we describe three local descent approaches which attempt to improve upon some starting solution. All three can be considered as simple descent strategies through a finite space of candidate solutions, a subset of all possible solutions but of quite different type for each method. The first method may be seen as a quick but ‘blind’ grid method, while the others are more involved, exploiting structural properties an optimal solution is known to satisfy. Observe that each objective value evaluation f(u, β) by (3) involves checking the sign of hu;ci − βfor each datapoint c∈A∪B. Therefore all of these values have to be calculated, regardless of whether they will be used or not. However, the values not used are exactly those that are to be used in the evaluation of the complementary half-space f(−u, −β). This means that each pair of complementary halfspaces may be evaluated together at virtually the same cost as a single evaluation of each of them. Therefore in all our codes both complementary halfspaces are always evaluated together and compared, so as to always choose the better one. 3.1 Grid-descent We consider regular grids in Rd\ {0} × Rdefined by some ‘maximum coefficient’ parameter m∈N. This search space H(m) consists of all (u, β) with integer coefficients in the interval [−m, m] (excluding those with u= 0). It should be noted, that H(m) does not necessarily contain an optimal separating halfspace. But thanks to the surjectivity of the map H<, an arbitrarily good approximation may be found in H(m) by choosing a sufficiently large m. On such a space we search as follows: Grid(m,∆) •Choose a random solution σ= (u, β)∈ H(m) •Sequentially for each steplength δstarting from ∆ and halving downto 1 –Cyclically for every coefficient, until a full cycle without change occurs ∗try adding and subtracting δfrom this coefficient ∗check when this solution is in H(m) and update σif it is better A one-parameter Grid(m) corresponds to Grid(m,m). Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 4 3.2 Cell-descent By extension of the results on median hyperplanes obtained in [13], it was shown in [14] that for any norm distance measure, there always exists an optimal halfspace determined by a hyperplane satisfying the following two properties: 1. it is blocked, i.e. passes through daffinely independent points of the training set A∪B. 2. it balances the misclassified datapoints, i.e. both for Aand Bthe number of their misclassified points cannot exceed the number of non-wellclassified points of the other class Since in Rdany daffinely independent points determine a unique hyperplane that passes through them, which generates two halfspaces, the set of blocked halfspaces Hbis finite. According to the first property we may restrict search to Hbwithout loss of optimality. The second property may be used for any fixed u6= 0 to find the best solution of type (u, β) by translation: choose βamong the set of values hu;A∪Bidef ={ h u;ci | c∈A∪B}in such a way that it balances the number of values of hu;Ailower than βagainst the number of values of hu;Bihigher than β. This can be easily done in O(plog p) by a single sweep after sorting hu;A∪Bi. Another method of linear complexity is described in [2]. Such a translation, even when started from a blocked hyperplane, does not necessarily result in a hyperplane in Hb. Therefore we also need a blocking step that allows to construct a blocked hyperplane starting from any σ∈ H, preferably one that does not deteriorate the objective value. This may be obtained by a cell-move as explained next. For any σ0∈ H we define the cell C(σ0)⊂ H as all halfspaces that classify points similarly to σ, or, more precisely, that do not classify correctly any points of A∪Bwhich were also misclassified by σ0, and do not misclassify any other points: C(σ0)def ={σ∈ H | A<(σ)⊂A≤(σ0), A≥(σ)⊂A≥(σ0), B>(σ)⊂B≥(σ0), B≤(σ)⊂B≤(σ0)} For a halfspace (u, β)∈ H to belong to this cell is expressed by the following linear inequalities: hu;ai ≤ β∀a∈A<(σ0) hu;ai ≥ β∀a6∈ A<(σ0) hu;bi ≥ β∀b∈B>(σ0) hu;bi ≤ β∀b6∈ B>(σ0) Note that σ0∈C(σ0). By (3) the constraint f(u, β) = f(σ0) (4) is a linear equality constraint not satisfied by (0,0), and each halfspace in C(σ0) has exactly one representative satisfying this constraint. Therefore this constraint may be added in the definition of C(σ0) without loss of generality, and it was proven in [14] that the polyhedral subset of Rd 0×R we then obtain is nonempty and bounded. We will also call it C(σ0). By (4) and (3) we see that minimising fon C(σ0) is equivalent to maximising γ◦(u) on C(σ0), and, by convexity of γ◦the optimum will be reached at some extreme point σ∗of C(σ0), and such aσ∗is always a blocked solution (for details see [14]). Furthermore, since σ0∈C(σ0), we will have f(σ∗)≤f(σ0), as sought. However, finding a maximum of γ◦(u) on C(σ0) is a hard global optimisation problem (see e.g. [6], chapter I.2), and therefore we propose to solve the following linear approximation instead. At σ0= (u0, β0) let p0∈∂γ◦(u0) be any subgradient of the dual norm γ◦at u0. Then hp0;ui= hp0;u−u0i+γ◦(u0)≤γ◦(u) for all u, because γ◦is a norm for which it is well-known that hp0;u0i=γ◦(u0). Therefore, maximising the linear function hp0;uion C(σ0) will also yield an extreme point of C(σ0), i.e. a blocked halfspace, with objective value no higher than f(σ0). In case finding a subgradient p0∈∂γ◦(u0) is not easy, one may use the direction of increase p0=u0of γ◦at u0instead. Note that for the euclidean norm this choice is a positive multiple Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 5 of a subgradient, so will perform as expected. For general norms, however, this does not fully guarantee that the new extreme point that will be obtained by solving the LP does not deteriorate the objective as compared to σ0. We can now describe our cell-descent as follows: Cell-Descent •Choose drandom points of A∪B, and find the σ= (u, β)∈ Hbpassing through these points. •Repeat until no new solution is found –Construct by translation the optimal solution σ0= (u0, β0) with fixed u0=u. –Choose p0∈∂γ◦(u0) or, if not available, take p0=u0 –Construct a new solution (u, β)∈ Hbby solving the following LP: max hp0;ui hu;ai ≤ β∀a∈A<(σ0) hu;ai ≥ β∀a6∈ A<(σ0) hu;bi ≥ β∀b∈B>(σ0) hu;bi ≤ β∀b6∈ B>(σ0) f(u, β) = f(σ0) u∈Rd, β ∈R 3.3 Translation-descent One may also consider the following much simpler local search method, solely based on translation. Translation-Descent •Choose a random solution σ= (u, β)∈ H(m) •Construct by translation the optimal solution σ0= (u0, β0) with fixed u0=u. This descent-search will be inoperative if started on a balanced solution, and always ends with a balanced solution. Therefore it is useless to try to repeat it. Grid-descent includes several trials of changes of the β-coefficient in a solution, which may be seen as an approximate form of Translation-descent. Cell-descent fully includes a Translation descent step in each of its loops. One may conclude that Cell-descent will certainly be more powerful than Translation-descent, while Grid descent will probably be so. The algorithm described by Karam et al [8] is based on such a translation-descent. However, it does not use the balancing property to find the best β-value for fixed u, but rather uses an updating method of the objective while sweeping the sorted sets hu;Aiand hu;Bi. 3.4 Comparison of descent-methods All three search methods were first coded in Matlab 7. For solving the LP’s in Cell-descent the built-in function ‘linprog’ could not be used, because it sometimes gave unexpected errors, and systematically malfunctioned when over 500 constraints were present. Therefore our code called CPLEX 9 by way of the CPLEXINT library [4]. Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 6 Tests were performed on the six datasets derived from the UCI Machine Learning Repository [12] for which Audet et al [1] published exact optimal solutions, unfortunately with only 4 significant digits. The details how these data sets were produced are given in that reference. For a correct interpretation of the results it is useful to know that all data sets, including the artificially produced ones discussed later, are always linearly standardised to [0,1]. All our tests are based on the euclidean distance L2. 100 runs were done on each of the six datasets. The average results of these 100 runs are listed in table 1. Here ddenotes the dimension of the data-space, pgives the size of the dataset, Trans descent, Grid(1000) and Cell-descent shows the resulting objective value when using Translation descent, Grid-descent with m= ∆ = 1000 or Cell-descent respectively, and the final column gives the global optimum value as published in [1]. Table 1: Descent methods: average results over 100 trials (Matlab) Data set d p Trans descent Grid(1000) Cell descent Global optimum Cancer 9 683 15.207 2.176 3.249 2.067 Diabetes 8 768 35.51 12.65 28.11 12.24 Echocardiogram 7 74 4.289 1.768 1.699 1.207 Glass 9 214 4.5548 0.60043 0.20693 0.03114 Housing 13 506 25.687 2.7770 4.2598 0.8971 Hepatitis 16 150 11.253 2.1717 1.6327 0.8711 As expected, one observes that Translation-descent always gives worse results than both other descent methods. Grid(1000) and Cell-descent work much better, and give comparable results. They even seem to be somewhat complementary. Both methods remain however far from obtaining consistently a good approximation to the global optimum. It is therefore clear that both descent mechanisms should be extended by a global search framework. 3.5 Comparison of random search methods Table 2 shows the best results obtained within the same 100 runs. This may be interpreted as results of repeated (100 trials) local search methods based on each of the three descent methods. Table 2: Random search methods (100 trials): results (Matlab) Data set d p Trans search Grid(1000) search Cell search Global optimum Cancer 9 683 3.277 2.072 2.079 2.067 Diabetes 8 768 16.65 12.26 15.42 12.24 Echocardiogram 7 74 2.182 1.218 1.318 1.207 Glass 9 214 1.2604 0.11016 0.03275 0.03114 Housing 13 506 9.9316 1.0377 0.9685 0.8971 Hepatitis 16 150 6.5964 1.2494 0.8788 0.8711 The Translation descent-based search remains quite poor. Grid-search and Cell-search now produce quite good results, but still not systematically. Observe, however, that in most cases one of the two search methods finds a solution which is quite close to optimal, but it is not always the same one. The computation times for translation search were in the two-hundredth to one-tenth of second range, for Grid search in the ten to two-hundred seconds range, and for Cell search in the two-tenth to three seconds range. It must be observed that our first implementations of Cell-search using the Matlab function linprog took considerably more time even than Grid-search. Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 7 4 Variable Neighbourhood Searches 4.1 General framework We propose to use the metaheuristic framework of Variable Neighbourhood Search (VNS) [10, 5], described in general below. However, in all of our computational testing we use the simpler Variable Neighbourhood Descent (VND), which consists of a single main loop of VNS (obtained by choosing as stopping condition simply ‘True’). VNS Initialization: •select the set of neighbourhood structures Nk(k= 1, . . . , kmax) that will be used in the search •find an initial solution x •choose a stopping condition Main loop: Repeat the following sequence until the stopping condition is met •Set kto 1. •Repeat the following steps until k > kmax Shaking: generate a point xat random from the kth neighbourhood of x(x∈Nk(x)); Local search: apply some local descent method with xas initial solution, ending at a local optimum x. Move or not if xis better than the incumbent x, then move there (i.e. set xto x), set kto 1 otherwise, set kto k+ 1 By specifying the yet undefined details indicated in italics, we obtain several different heuristics. 4.2 GridVNS A first method GridVNS(m), operating in the searchspace H(m), uses following specifications: Initial solution a random halfspace with integer coefficients in the range [−m, m] Neighbourhoods Nk(σ) contains all halfspaces having d+ 1 −kcoefficients in common with σ and knew ones. We also take kmax =d+ 1 Local descent Grid(m,d√me), and using the current solution as initial solution The particular choice of ∆ = d√mewas dictated by our concern to give a more local search character to the inner Grid runs. A first series of tests with this scheme did not give satisfactory results. Furthermore it was observed that by far most improved solutions were found for very low k-values. We therefore implemented four variants of this search strategy: GridNVNS(m)or Narrow Grid VNS: this is the standard Grid VND as described above (thus with a single main loop). GridBVNS(m)or Broad Grid VNS: now in the innermost loop each value of kis repeated several times, in such a way that the total number of coefficients modified while working within Nk is a constant (set to kmax). This means that k= 1 is used kmax times, k= 2 is used kmax/2 times, etc. Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 8 2S-NVNS(m2)or two-stage Grid NVNS: first GridNVNS(m) is executed , and followed by a Grid(m2,m2) run 2S-BVNS(m2)as the previous, but using the broad version GridBVNS(m) We did a comparative test of the following methods: GridNVNS(100), GridBVNS(100), 2SNVNS(1000), 2S-BVNS(1000), in other words, the two-stage methods work first on a coarser grid and secondly on a finer grid than in their one-step variants. Table 3 shows the average results obtained on the same data sets as before during 5 independent runs of each method, and includes again for comparison the ‘exact’ global optimal values given by [1]. The best heuristically found average value is shown in boldface. The last column indicates the relative deviation of this best value with respect to the ‘true’ optimum. Table 3: GridVNS: average results over 5 runs (Matlab) Data set d p NVNS BVNS 2S-NVNS 2S-BVNS Global opt. error Cancer 9 683 2.086 2.078 2.081 2.081 2.067 0.5% Diabetes 8 768 12.41 12.30 12.45 12.39 12.24 0.5% Echocardiogram 7 74 1.267 1.232 1.258 1.286 1.207 2.1% Glass 9 214 0.2880 0.1370 0.1931 0.1688 0.03114 340.0% Housing 13 506 2.198 2.175 2.450 2.201 0.8971 142.4% Hepatitis 16 150 1.122 1.053 1.415 1.157 0.8711 20.9% We can observe that the best results were systematically obtained with GridBVNS. The improvement of BVNS upon NVNS indicates that the more intense search within small neighbourhoods is effective, whereas the second stage using a final finer grid-descent seems almost useless. Observe also that the quality of the solutions found on the first three data sets is relatively good, but quite bad for the three last data sets. These conclusions drawn from the objective values reached should be somewhat revised in view of the computational times taken by the different methods, as shown in table 4. Table 4: GridVNS: average computation times (s) over 5 runs (Matlab) Data set d p NVNS BVNS 2S-NVNS 2S-BVNS Cancer 9 683 205 608 193 442 Diabetes 8 768 349 700 271 406 Echocardiogram 7 74 27 56 17 33 Glass 9 214 69 197 58 126 Housing 13 506 875 1537 394 983 Hepatitis 16 150 175 864 125 388 4.3 PointVNS The second method PointVNS operates in the searchspace of blocked halfspaces Hb. Recall that this means that every halfspace is determined by daffinely independent datapoints, and the choice of an orientation. Initial solution a random hyperplane going through daffinely independent data points Neighbourhoods Nk(σ) contains all halfspaces determined by d−kdata points common with σand knew ones. Evidently kmax =d. Local descent Cell descent, with the current solution as initial solution Plastria/De Bruyne/Carrizosa / Alternating local search based VNS for linear classification 9 Note that choosing a new solution in Nk(σ) is not so simple. We have encountered many difficulties due to degenerate situations. Indeed, simply replacing kdatapoints that determine σby kother datapoints often yielded a (nearly) affinely dependent set of points, which either led to numerical difficulties or did not define a single hyperplane. Therefore we had to include affine dependency tests in the code which was quite detrimental to its computational efficiency. We have implemented two main variants of this search strategy: PointNVNS or narrow Point VNS: this is the standard VND as described above. PointBVNS(r)or broad Point VNS: now in the innermost loop each value of kis repeated as long as the total number of datapoints modified while working within Nkdoes not exceed rd. Each method was coded in Matlab and run 5 times on each data set. Table 5 shows the results obtained when applying the first and three instances of the second variant (taking k= 1,3,10) on the same datasets as before, using the same presentation as in table 3. Table 5: PointVNS: average results over 5 runs (Matlab) Data set d p NVNS BVNS(1) BVNS(3) BVNS(10) Global opt. error Cancer 9 683 2.279 2.131 2.096 2.101 2.067 1.4% Diabetes 8 768 16.74 14.93 14.28 13.49 12.24 10.2% Echocardiogram 7 74 1.455 1.321 1.225 1.233 1.207 1.5% Glass 9 214 0.03262 0.03191 0.03154 0.03147 0.03114 1.1% Housing 13 506 1.045 0.9956 0.9575 0.9194 0.8971 2.5% Hepatitis 16 150 0.9918 0.9085 0.8917 0.8854 0.8711 1.6% Here we observe that the broad VNS methods are always better than the narrow version. Usually (but not systematically) the most intensive search strategy in each neighbourhood gives the best results, as was to be expected. Concerning the quality, one may see that it is usually quite good, with a notable exception on the second data set. Table 6: PointVNS: average computation times (s) over 5 runs (Matlab) Data set d p NVNS BVNS(1) BVNS(3) BVNS(10) Cancer 9 683 86 301 675 3724 Diabetes 8 768 47 149 349 1097 Echocardiogram 7 74 3 8 14 104 Glass 9 214 17 55 204 486 Housing 13 506 108 147 936 2728 Hepatitis 16 150 9 47 83 295 The corresponding average computational times are given in table 6. Computation times for Point-BVNS(k) increase almost linearly with k. 5 Hybrid method: Grid-Cell As compared to the results obtained with GridBVNS, shown in table 3, the results found by all PointVNS, shown in table 5, are worse for the two first data sets, but very much better for the three last ones. This suggests that the two approaches GridVNS and PointVNS are complementary: when one performs badly, the other one performs well and vice-versa. This prompted us to combine their features into a single hybrid method.