Multi-Objective Evolutionary Rule-Based Classification with Categorical Data

The ease of interpretation of a classification model is essential for the task of validating it. Sometimes it is required to clearly explain the classification process of a model’s predictions. Models which are inherently easier to interpret can be effortlessly related to the context of the problem, and their predictions can be, if necessary, ethically and legally evaluated. In this paper, we propose a novel method to generate rule-based classifiers from categorical data that can be readily interpreted. Classifiers are generated using a multi-objective optimization approach focusing on two main objectives: maximizing the performance of the learned classifier and minimizing its number of rules. The multi-objective evolutionary algorithms ENORA and NSGA-II have been adapted to optimize the performance of the classifier based on three different machine learning metrics: accuracy, area under the ROC curve, and root mean square error. We have extensively compared the generated classifiers using our proposed method with classifiers generated using classical methods such as PART, JRip, OneR and ZeroR. The experiments have been conducted in full training mode, in 10-fold cross-validation mode, and in train/test splitting mode. To make results reproducible, we have used the well-known and publicly available datasets Breast Cancer, Monk’s Problem 2, Tic-Tac-Toe-Endgame, Car, kr-vs-kp and Nursery. After performing an exhaustive statistical test on our results, we conclude that the proposed method is able to generate highly accurate and easy to interpret classification models.


Introduction
Supervised Learning is the branch of Machine Learning (ML) [1] focused on modeling the behavior of systems that can be found in the environment. Supervised models are created from a set of past records, each one of which, usually, consists of an input vector labeled with an output. A supervised model is an algorithm that simulates the function that maps inputs with outputs [2]. The best models are those that predict the output of new inputs in the most accurate way. Thanks to modern computing capabilities, and to the digitization of ever-increasing quantities of data, nowadays, supervised learning techniques play a leading role in many applications. The first classification systems date back to the 1990s; in those days, researchers were focused on both precision and interpretability, and the systems to be modeled were relatively simple. Years later, when it became necessary to model more difficult behaviors, the researchers focused on developing more and more precise models, leaving aside the interpretability. Artificial Neural Networks (ANN) [3], and, more recently, Deep Learning Neural Networks (DLNN) [4], as well as Support Vector Machines (SVM) [5], and Instance-based Learning (IBL) [6] are archetypical examples of this approach. A DLNN, for example, is a large mesh of ordered nodes arranged in a hierarchical manner and composed of a huge number of variables. DLNNs are capable of modeling very complex behaviors, but it is extremely difficult to understand the logic behind their predictions, and similar considerations can be drawn for SVNs and IBLs, although the underlying principles are different. These models are known as black-box methods. While there are applications in which knowing the ratio behind a prediction is not necessarily relevant, (e.g., predicting a currency's future value, whether or not a user clicks on an advert or the amount of rain in a certain area), there are other situations where the interpretability of a model plays a key role.
The interpretability of classification systems refers to the ability they have to explain their behavior in a way that is easily understandable by a user [7]. In other words, a model is considered interpretable when a human is able to understand the logic behind its prediction. In this way, Interpretable classification models allow external validation by an expert. Additionally, there are certain disciplines such as medicine, where it is essential to provide information about decision making for ethical and human reasons. Likewise, when a public institution asks an authority for permission to investigate an alleged offender, or when the CEO of a certain company wants to take a difficult decision which can seriously change the direction of the company, some kind of explanations to justify these decisions may be required. In these situations, using transparent (also called grey-box) models is recommended. While there is a general consensus on how the performance of a classification system is measured (popular metrics include accuracy, area under the ROC curve, and root mean square error), there is no universally accepted metric to measure the interpretability of the models. Nor is there an ideal balance between the interpretability and performance of classification systems but this depends on the specific application domain. However, the rule of thumb says that the simpler a classification system is, the easier it is to interpret. Rule-based Classifiers (RBC) [8,9] are among the most popular interpretable models, and some authors define the degree of interpretability of an RBC as the number of its rules or as the number of axioms that the rules have. These metrics tend to reward models with fewer rules as simple as possible [10,11]. In general, RBCs are classification learning systems that achieve a high level of interpretability because they are based on a human-like logic. Rules follow a very simple schema: IF (Condition 1) and (Condition 2) and . . . (Condition N) THEN (Statement) and the fewer rules the models have and the fewer conditions and attributes the rules have, the easier it will be for a human to understand the logic behind each classification. In fact, RBCs are so natural in some applications that they are used to interpret other classification models such as Decision Trees (DT) [12]. RBCs constitute the basis of more complex classification systems based on fuzzy logic [13] such as LogitBoost or AdaBoost [14].
Our approach investigates the conflict between accuracy and interpretability as a multi-objective optimization problem. We define a solution as a set of rules (that is, a classifier), and establish two objectives to be maximized: interpretability and accuracy. We decided to solve this problem by applying multi-objective evolutionary algorithms (MOEA) [15,16] as meta-heuristics, and, in particular, two known algorithms: NSGA-II [15] and ENORA [17]. They are both state-of-the-art evolutionary algorithms which have been applied, and compared, on several occasions [18][19][20]. NSGA-II is very well-known and has the advantage of being available in many implementations, while ENORA generally has a higher performance. In the current literature, MOEAs are mainly used for learning RBCs based on fuzzy logic [18,[21][22][23][24][25][26]. However, Fuzzy RBCs are designed for numerical data, from which fuzzy sets are constructed and represented by linguistic labels. In this paper, on the contrary, we are interested in RBCs for categorical data, for which a novel approach is necessary. This paper is organized as follows. In Section 2, we introduce multi-objective constrained optimization, the evolutionary algorithms ENORA and NSGA-II, and the well-known rule-based classifier learning systems PART, JRip, OneR and ZeroR. In Section 3, we describe the structure of an RBC for categorical data, and we propose the use of multi-objective optimization for the task of learning a classifier. In Section 4, we show the result of our experiments, performed on the well-known publicly accessible datasets Breast Cancer, Monk's Problem 2, Tic-Tac-Toe-Endgame, Car, kr-vs-kp and Nursery. The experiments allow a comparison among the performance of the classifiers learned by our technique against those of classifiers learned by PART, JRip, OneR and ZeroR, as well as a comparison between ENORA and NSGA-II for the purposes of this task. In Section 5, the results are analyzed and discussed, before concluding in Section 6. Appendices A and B show the tables of the statistical tests results. Appendix C shows the symbols and the nomenclature used in the paper.

Multi-Objective Constrained Optimization
The term optimization [27] refers to the selection of the best element, with regard to some criteria, from a set of alternative elements. Mathematical programming [28] deals with the theory, algorithms, methods and techniques to represent and solve optimization problems. In this paper, we are interested in a class of mathematical programming problems called multi-objective constrained optimization problems [29], which can be formally defined, for l objectives and m constraints, as follows: where f i (x) (usually called objectives) and g j (x) are arbitrary functions. Optimization problems can be naturally separated into two categories: those with discrete variables, which we call combinatorial, and those with continuous variables. In combinatorial problems, we are looking for objects from a finite, or countably infinite, set X , where objects are typically integers, sets, permutations, or graphs. In problems with continuous variables, instead, we look for real parameters belonging to some continuous domain. In Equation (1), x = {x 1 , x 2 , . . . , x w } ∈ X w represents the set of decision variables, where X is the domain for each variable x k , k = 1, . . . , w. Now, let F = {x ∈ X w | g j (x) ≤ 0, j = 1, . . . , m} be the set of all feasible solutions to Equation (1). We want to find a subset of solutions S ⊆ F called non-dominated set (or Pareto optimal set). A solution x ∈ F is non-dominated if there is no other solution x ∈ F that dominates x, and a solution x dominates x if and only if there exists i (1 ≤ i ≤ l) such that f i (x ) improves f i (x), and for every i (1 ≤ i ≤ l), f i (x) does not improve f i (x ). In other words, x dominates x if and only if x is better than x for at least one objective, and not worse than x for any other objective. The set S of non dominated solutions of Equation (1) can be formally defined as: Once the set of optimal solutions is available, the most satisfactory one can be chosen by applying a preference criterion. When all the functions f i are linear, then the problem is a linear programming problem [30], which is the classical mathematical programming problem and for which extremely efficient algorithms to obtain the optimal solution exist (e.g., the simplex method [31]). When any of the functions f i is non-linear then we have a non-linear programming problem [32]. A non-linear programming problem in which the objectives are arbitrary functions is, in general, intractable. In principle, any search algorithm can be used to solve combinatorial optimization problems, although it is not guaranteed that they will find an optimal solution. Metaheuristics methods such as evolutionary algorithms [33] are typically used to find approximate solutions for complex multi-objective optimization problems, including feature selection and fuzzy classification.

The Multi-Objective Evolutionary Algorithms ENORA and NSGA-II
The MOEA ENORA [17] and NSGA-II [15] use a (µ + λ) strategy (Algorithm 1) with µ = λ = popsize, where µ corresponds to the number of parents and λ refers to the number of children (popsize is the population size), with binary tournament selection (Algorithm 2) and a rank function based on Pareto fronts and crowding (Algorithms 3 and 4). The difference between NSGA-II and ENORA is how the calculation of the ranking of the individuals in the population is performed. In ENORA, each individual belongs to a slot (as established in [34]) of the objective search space, and the rank of an individual in a population is the non-domination level of the individual in its slot. On the other hand, in NSGA-II, the rank of an individual in a population is the non-domination level of the individual in the whole population. Both ENORA and NSGA-II MOEAs use the same non-dominated sorting algorithm, the fast non-dominated sorting [35]. It compares each solution with the rest of the solutions and stores the results so as to avoid duplicate comparisons between every pair of solutions. For a problem with l objectives and a population with N solutions, this method needs to conduct l · N · (N − 1) objective comparisons, which means that it has a time complexity of O(l · N 2 ) [36]. However, ENORA distributes the population in N slots (in the best case), therefore, the time complexity of ENORA is O(l · N 2 ) in the worst case and O(l · N) in the best case. Algorithm 1 (µ + λ) strategy for multi-objective optimization.  1: if rank (P, I) < rank (P, J) then 2: return True 3: end if 4: if rank(P, J) < rank(P, I) then 5: return False 6: end if 7: return Crowding_distance (P, I) > Crowding_distance (P, J)

Require
The main reason ENORA and NSGA-II behave differently is as follows. NSGA-II never selects the individual dominated by the other in the binary tournament, while, in ENORA, the individual dominated by the other may be the winner of the tournament. Figure 1 shows this behavior graphically. For example, if individuals B and C are selected for a binary tournament with NSGA-II, individual B beats C because B dominates C. Conversely, individual C beats B with ENORA because individual C has a better rank in his slot than individual B. In this way, ENORA allows the individuals in each slot to evolve towards the Pareto front encouraging diversity. Even though in ENORA the individuals of each slot may not be the best of the total individuals, this approach generates a better hypervolume than that of NSGA-II throughout the evolution process.
ENORA is our MOEA, on which we are intensively working over the last decade. We have applied ENORA to constrained real-parameter optimization [17], fuzzy optimization [37], fuzzy classification [18], feature selection for classification [19] and feature selection for regression [34]. In this paper, we apply it to rule-based classification. NSGA-II algorithm was designed by Deb et al. and has been proved to be a very powerful and fast algorithm in multi-objective optimization contexts of all kinds. Most researchers in multi-objective evolutionary computation use NSGA-II as a baseline to compare the performance of their own algorithms. Although NSGA-II was developed in 2002 and remains a state-of-the-art algorithm, it is still a challenge to improve on it. There is a recently updated improved version for many-objective optimization problems called NSGA-III [38].

PART
PART (Partial DT Method [39]) is a widely used rule learning algorithm that was developed at the University of Waikato in New Zealand [40]. Experiments show that it is a very efficient algorithm in terms of both computational performance and results. PART combines the divide-and-conquer strategy typical of decision tree learning with the separate-and-conquer strategy [41] typical of rule learning, as follows. A decision tree is first constructed (using C4.5 algorithm [42]), and the leaf with the highest coverage is converted into a rule. Then, the set of instances that are covered by that rule are discarded, and the process starts over. The result is an ordered set of rules, completed by a default rule that applies to instances that do not meet any previous rule.

JRip
JRip is a fast and optimized implementation in Weka of the famous RIPPER (Repeated Incremental Pruning to Produce Error Reduction) algorithm [43]. RIPPER was proposed in [44] as a more efficient version of the incrementally reduced error pruning (IREP) rule learner developed in [45]. IREP and RIPPER work in a similar manner. They begin with a default rule and, using a training dataset, attempt to learn rules that predict exceptions to the default. Each rule learned is a conjunction of propositional literals. Each literal corresponds to a split of the data based on the value of a single feature. This family of algorithms, similar to decision trees, has the advantage of being easy to interpret, and experiments show that JRip is particularly efficient in large datasets. RIPPER and IREP use a strategy based on the separate-and-conquer method to generate an ordered set of rules that are extracted directly from the dataset. The classes are examined one by one, prioritizing those that have more elements. These algorithms are based on four basic steps (growing, pruning, optimizing and selecting) applied repetitively to each class until a stopping condition is met [44]. These steps can be summarized as follows. In the growing phase, rules are created taking into account an increasing number of predictors until the stopping criterion is satisfied (in the Weka implementation, the procedure selects the condition with the highest information gain). In the pruning phase redundancy is eliminated and long rules are reduced. In the optimization phase, the rules generated in the previous steps are improved (if possible) by adding new attributes or by adding new rules. Finally, in the selection phase, the best rules are selected and the others discarded.

OneR
OneR (One Rule) is a very simple, while reasonably accurate, classifier based on a frequency table. First, OneR generates a set of rules for each attribute of the dataset, and, then, it selects only one rule from that set-the one with the lowest error rate [46]. The set of rules is created using a frequency table constructed for each predictor of the class, and numerical classes are converted into categorical values.

ZeroR
Finally, ZeroR (Zero Rules [40]) is a classifier learner that does not create any rules and uses no attributes. ZeroR simply creates the class classification table by selecting the most frequent value. Such a classifier is obviously the simplest possible one, and its capabilities are limited to the prediction of the majority class. In the literature, it is not used for practical classifications tasks, but as a generic reference to measure the performance of other classifiers.

Multi-Objective Optimization for Categorical Rule-Based Classification
In this section, we propose a general schema for an RBC specifically designed for categorical data. Then, we propose and describe a multi-objective optimization solution to obtain optimal categorical RBCs.

Rule-Based Classification for Categorical Data
Let Γ be a classifier composed by M rules, where each rule R Γ i , i = 1, . . . , M, has the following structure: where for j = 1, . . . , p the attribute b Γ ij (called antecedent) takes values in a set {1, . . . , v j } (v j > 1), and c Γ i (called consequent) takes values in {1, . . . , w} (w > 1). Now, let x = {x 1 , . . . , x p } be an observed example, with x j ∈ {1, . . . , v j }, for each j = 1, . . . , p. We propose maximum matching as reasoning method, where the compatibility degree of the rule R Γ i for the example x (denoted by ϕ Γ i (x)) is calculated as the number of attributes whose value coincides with that of the corresponding antecedent in R Γ i , that is where: The association degree for the example x with a class c ∈ {1, . . . , w} is computed by adding the compatibility degrees for the example x of each rule R Γ i whose consequent c Γ i is equal to class c, that is: where: Therefore, the classification (or output) of the classifier Γ for the example x corresponds to the class whose association degree is maximum, that is:

A Multi-Objective Optimization Solution
Let D be a dataset of K instances with p categorical input attributes, p > 0, and a categorical output attribute. Each input attribute j can take a category x j ∈ 1, . . . , v j , v j > 1, j = 1, . . . , p, and the output attribute can take a class c ∈ {1, . . . , w}, w > 1. The problem of finding an optimal classifier Γ, as described in the previous section, can be formulated as an instance of the multi-objective constrained problem in Equation (1) with two objectives and two constraints: In the problem (Equation (3)), the function F D (Γ) is a performance measure of the classifier Γ over the dataset D, the function N R(Γ) is the number of rules of the classifier Γ, and the constraints N R(Γ) ≥ w and N R(Γ) ≤ M max limit the number of rules of the classifier Γ to the interval [w, M max ], where w is the number of classes of the output attribute and M max is given by a user. Objectives F D (Γ) and N R(Γ) are in conflict. The fewer rules the classifier has, the fewer instances it can cover, that is, if the classifier is simpler it will have less capacity for prediction. There is, therefore, an intrinsic conflict between problem objectives (e.g., maximize accuracy and minimize model complexity) which cannot be easily aggregated to a single objective. Both objectives are typically optimized simultaneously in many other classification systems, such as neural networks or decision trees [47,48]. Figure 2 shows the Pareto front of a dummy binary classification problem described as in Equation (3), with M max = 6 rules, where F D (Γ) is maximized. This front is composed of three non-dominated solutions (three possible classifiers) with two, three and four rules, respectively. The solutions with five and six rules are dominated (both by the solution with four rules). Both ENORA and NSGA-II have been adapted to solve the problem described in Equation (3) with variable-length representation based on a Pittsburgh approach, uniform random initialization, binary tournament selection, handling constraints, ranking based on non-domination level with crowding distance, and self-adaptive variation operators. Self-adaptive variation operators work on different levels of the classifier: rule crossover, rule incremental crossover, rule incremental mutation, and integer mutation.

Representation
We use a variable-length representation based on a Pittsburgh approach [49], where each individual I of a population contains a variable number of rules M I , and each rule R I i , i = 1, . . . , is codified in the following components: • Integer values associated to the antecedents b I ij ∈ {1, . . . , v j }, for i = 1, . . . , M I and j = 1, . . . , p.

•
Integer values associated to the consequent c I i ∈ {1, . . . , w}, for i = 1, . . . , M I . Additionally, to carry out self-adaptive crossing and mutation, each individual has two discrete parameters d I ∈ {0, . . . , δ} and e I ∈ {0, . . . , } associated with crossing and mutation, where δ ≥ 0 is the number of crossing operators and ε ≥ 0 is the number of mutation operators. Values d I and e I for self-adaptive variation are randomly generated from {0, δ} and {0, }, respectively. Table 1 summarizes the representation of an individual.

Constraint Handling
The constraints N R(Γ) ≥ w and N R(Γ) ≤ M max are satisfied by means of specialized initialization and variation operators, which always generate individuals with a number of rules between w and M max .

Initial Population
The initial population (Algorithm 5) is randomly generated with the following conditions:

•
Individuals are uniformly distributed with respect to the number of rules with values between w and M max , and with an additional constraint that specifies that there must be at least one individual for each number of rules (Steps 4-8). This ensures an adequate initial diversity in the search space in terms of the second objective of the optimization model.

•
All individuals contain at least one rule for any output class between 1 and w (Steps 16-20).
Algorithm 5 Initialize population. if k ≤ M max − w + 1 then 5: end if 9: {Random rule R I i } 10: for i = 1 to M I do 11: {Random integer values associated with the antecedents} 12: for j = 1 to p do 13: b I ij ←Random(1,v j ) 14: end for 15: {Random integer value associated with the consequent} 16: if i < w then 17:

Fitness Functions
Since the optimization model encompasses two objectives, each individual must be evaluated with two fitness functions, which correspond to the objective functions F D (Γ) and N R(Γ) of the problem (Equation (3)). The selection of the best individuals is done using the Pareto concept in a binary tournament.

Variation Operators
We use self-adaptive crossover and mutation, which means that the selection of the operators is made by means of an adaptive technique. As we have explained (cf. Section 3.2.1), each individual I has two integer parameters d I ∈ {0, . . . , δ} and e I ∈ {0, . . . , } to indicate which crossover or mutation is carried out. In our case, δ = 2 and = 2 are two crossover operators and two mutation operators, so that d I , e I ∈ {0, 1, 2}. Note that value 0 indicates that no crossover or no mutation is performed. Self-adaptive variation (Algorithm 6) generates two children from two parents by self-adaptive crossover (Algorithm 7) and self-adaptive mutation (Algorithm 8). Self-adaptive crossover of individuals I, J and self-adaptive mutation of individual I are similar to each other. First, with a probability p v , the values d I and e I are replaced by a random value. Additionally, in the case of crossover, the value d J is replaced by d I . Then, the crossover indicated by d I or the mutation indicated by e I is performed. In summary, if an individual comes from a given crossover or a given mutation, that specific crossover and mutation are preserved to their offspring with probability p v , so the value of p v must be small enough to ensure a controlled evolution (in our case, we use p v = 0.1). Although the probability of the crossover and mutation is not explicitly represented, it can be computed as the ratio of the individuals for which crossover and mutation values are set to 1. As the population evolves, individuals with more successful types of crossover and mutation will be more common, so that the probability of selecting the more successful crossover and mutation types will increase. Using self-adaptive crossover and mutation operators helps to realize the goals of maintaining diversity in the population and sustaining the convergence capacity of the evolutionary algorithm, also eliminating the need of setting an a priori operator probability to each operator. In other approaches (e.g., [50]), the probabilities of crossover and mutation vary depending on the fitness value of the solutions.
Both ENORA and NSGA-II have been implemented with two crossover operators, rule crossover (Algorithm 9) and rule incremental crossover (Algorithm 10), and two mutation operators: rule incremental mutation (Algorithm 11) and integer mutation (Algorithm 12). Rule crossover randomly exchanges two rules selected from the parents, and rule incremental crossover adds to each parent a rule randomly selected from the other parent if its number of rules is less than the maximum number of rules. On the other hand, rule incremental mutation adds a new rule to the individual if the number of rules of the individual is less than the maximum number of rules, while integer mutation carries out a uniform mutation of a random antecedent belonging to a randomly selected rule.

Experiment and Results
To ensure the reproducibility of the experiments, we have used publicly available datasets. In particular, we have designed two sets of experiments, one using the Breast Cancer [51] dataset, and the other using the Monk's Problem 2 [52] dataset.

The Breast Cancer Dataset
Breast Cancer encompasses 286 instances. Each instance corresponds to a patient who suffered from breast cancer and uses nine attributes to describe each patient. The class to be predicted is binary and represents whether the patient has suffered a recurring cancer event. In this dataset, 85 instances are positive and 201 are negative. Table 2 summarizes the attributes of the dataset. Among all instances, nine present some missing values; in the pre-processing phase, these have been replaced by the mode of the corresponding attribute. In July 1991, the monks of Corsendonk Priory attended a summer course that was being held in their priory, namely the 2nd European Summer School on Machine Learning. After a week, the monks could not yet clearly identify the best ML algorithms, or which algorithms to avoid in which cases. For this reason, they decided to create the three so-called Monk's problems, and used them to determine which ML algorithms were the best. These problems, rather simple and completely artificial, became later famous (because of their peculiar origin), and have been used as a comparison for many algorithms on several occasions. In particular, in [53], they have been used to test the performance of state-of-the-art (at that time) learning algorithms such as AQ17-DCI, AQ17-HCI, AQ17-FCLS, AQ14-NT,  AQ15-GA, Assistant Professional, mFOIL, ID5R, IDL, ID5R-hat, TDIDT, ID3, AQR, CN2, WEB CLASS, ECOBWEB, PRISM, Backpropagation, and Cascade Correlation. For our research, we have used the Monk's Problem 2, which contains six categorical input attributes and a binary output attribute, summarized in Table 3. The target concept associated with the Monk's Problem 2 is the binary outcome of the logical formula: Exactly two of: {heap_shape= round, body_shape=round, is_smiling=yes, holding=sword, jacket_color=red, has_tie=yes} In this dataset, the original training and testing sets were merged to allow other sampling procedures. The set contains a total of 601 instances, and no missing values.

Optimization Models
We have conducted different experiments with different optimization models to calculate the overall performance of our proposed technique and to see the effect of optimizing different objectives for the same problem. First, we have designed a multi-objective constrained optimization model based on the accuracy: Max.
where ACC D (Γ) is the proportion of correctly classified instances (both true positives and true negatives) among the total number of instances [54] obtained with the classifier Γ for the dataset D. ACC D (Γ) is defined as: where K is the number of instances of the dataset D, and T D (Γ, i) is the result of the classification of the instance i in D with the classifier Γ, that is: whereĉ Γ i is the predicted value of the ith instance in Γ, and c i D is the corresponding true value in D. Our second optimization model is based on the area under the ROC curve: where AU C D (Γ) is the area under the ROC curve obtained with the classifier Γ with the dataset D.
The ROC (Receiver Operating Characteristic) curve [55] is a graphical representation of the sensitivity versus the specificity index for a classifier varying the discrimination threshold value. Such a curve can be used to generate statistics that summarize the performance of a classifier, and it has been shown in [54] to be a simple, yet complete, empirical description of the decision threshold effect, indicating all possible combinations of the relative frequencies of the various kinds of correct and incorrect decisions.
The area under the ROC curve can be computed as follows [56]: where S D (Γ, t) (sensitivity) is the proportion of positive instances classified as positive by the classifier Γ in D, 1 − E D (Γ, t) (specificity) is the proportion of negative instances classified as negative by Γ in D, and t is the discrimination threshold. Finally, our third constrained optimization model is based on the root mean square error (RMSE): where RMSE D (Γ) is defined as the square root of the mean square error obtained with a classifier Γ in the dataset D: whereĉ Γ i is the predicted value of the ith instance for the classifier Γ, and c i D is the corresponding output value in the database D. Accuracy, area under the ROC curve, and root mean square error are all well-accepted measures used to evaluate the performance of a classifier. Therefore, it is natural to use such measures as fitting functions. In this way, we can establish which one behaves better in the optimization phase, and we can compare the results with those in the literature.

Choosing the Best Pareto Front
To compare the performance of ENORA and NSGA-II as metaheuristics in this particular optimization task, we use the hypervolume metric [57,58]. The hypervolume measures, simultaneously, the diversity and the optimality of the non-dominated solutions. The main advantage of using hypervolume against other standard measures, such as the error ratio, the generational distance, the maximum Pareto-optimal front error, the spread, the maximum spread, or the chi-square-like deviation, is that it can be computed without an optimal population, which is not always known [15]. The hypervolume is defined as the volume of the search space dominated by a population P, and is formulated as: where Q ⊆ P is the set of non-dominated individuals of P, and v i is the volume of the individual i. Subsequently, the hypervolume ratio (HVR) is defined as the ratio of the volume of the non-dominated search space over the volume of the entire search space, and is formulated as follows: where VS is the volume of the search space. Computing HVR requires reference points that identify the maximum and minimum values for each objective. For RBC optimization, as proposed in this work, the following minimum (F D lower , N R lower ) and maximum (F D upper , N R upper ) points, for each objective, are set in the multi-objective optimization models in Equations (4)- (6): A first single execution of all six models (three driven by ENORA, and three driven by NSGA-II), over both datasets, has been designed for the purpose of showing the aspect of the final Pareto front, and compare the hypervolume ratio of the models. The results of this single execution, with population size equal to 50 and 20,000 generations (1,000,000 evaluations in total), are shown in Figures 3 and 4 (by default, M max is set to 10, to which we add 2, because both datasets have a binary class). Regarding the configuration of the number of generations and the size of the population, our criterion has been established as follows: once the number of evaluations is set to 1,000,000, we can decide to use a population size of 100 individuals and 10,000 generations, or to use a population size of 50 individuals and 20,000 generations. The first configuration (100 × 10,000) allows a greater diversity with respect to the number of rules of the classifiers, while the second one (50 × 20,000) allows a better adjustment of the classifier parameters and therefore, a greater precision. Given the fact that the maximum number of rules of the classifiers is not greater than 12, we think that 50 individuals are sufficient to represent four classifiers on average for each number of rules (4 × 12 = 48∼50). Thus, we prefer the second configuration (50 × 20,000) because having more generations increases the chances of building classifiers with a higher precision.  Experiments were executed in a computer x64-based PC with one processor Intel64 Family 6 Model 60 Stepping 3 GenuineIntel 3201 Mhz, RAM 8131 MB. Table 4 shows the run time for each method over both datasets. Note that, although ENORA has less algorithmic complexity than NSGA-II, it has taken longer in experiments than NSGA-II. This is because the evaluation time of individuals in ENORA is higher than that of NSGA-II since ENORA has more diversity than NSGA-II, and therefore ENORA evaluates classifiers with more rules than NSGA-II. From these results, we can deduce that, first, ENORA maintains a higher diversity of the population, and achieves a better hypervolume ratio with respect to NSGA-II, and, second, using accuracy as the first objective generates better fronts than using the area under the ROC curve, which, in turn, performs better than using the root mean square error.

Comparing Our Method with Other Classifier Learning Systems (Full Training Mode)
To perform an initial comparison between the performance of the classifiers obtained with the proposed method and the ones obtained with classical methods (PART, JRip, OneR and ZeroR), we have executed again the six models in full training mode.
The parameters have been configured as in the previous experiment (population size equal to 50 and 20,000 generations), excepting the M max parameter that was set to 2 for the Breast Cancer dataset (this case), while, for the Monk's Problem 2, it was set to 9. Observe that, since M min = 2 in both cases, executing the optimization models using M max = 2 leads to a single objective search for the Breast Cancer dataset. In fact, after the preliminary experiments were run, it turned out that the classical classifier learning systems tend to return very small, although not very precise, set of rules on Breast Cancer, and that justifies our choice. On the other hand, executing the classical rule learners on Monk's Problem 2 returns more diverse sets of rules, which justifies choosing a higher M max in that case. To decide, a posteriori, which individual is chosen from the final front, we have used the default algorithm: the individual with the best value on the first objective is returned. In the case of Monk's Problem 2, that individual has seven rules. The comparison is shown in Tables 5 and 6, which show, for each classifier, the following information: number of rules, percent correct, true positive rate, false positive rate, precision, recall, F-measure, Matthews correlation coefficient, area under the ROC curve, area under precision-recall curve, and root mean square error. As for the Breast Cancer dataset (observe that the best result emerged from the proposed method), in the optimization model driven by NSGA-II, with root mean square error as the first objective (see Table 7), only PART was able to achieve similar results, although slightly worse, but at the price of having 15 rules, making the system clearly not interpretable. In the case of the Monk's Problem 2 dataset, PART returned a model with 47 rules, which is not interpretable by any standard, although it is very accurate. The best interpretable result is the one with seven rules returned by ENORA, driven by the root mean square error (see Table 8). The experiments for classical learners have been conducted using the default parameters.  To test the capabilities of our methodology in a more significant way, we proceeded as follows. First, we designed a cross-validated experiment for the Breast Cancer dataset, in which we iterated three times a 10-fold cross-validation learning process [59] and considered the average value of the performance metrics percent correct, area under the ROC curve, and serialized model size of all results. Second, we designed a train/test percentage split experiment for the Monk's Problem 2 dataset, in which we iterated ten times a 66% (training) versus 33% (testing) split and considered, again, the average result of the same metrics. Finally, we performed a statistical test over on results, to understand if they show any statistically significant difference. An execution of our methodology, and of standard classical learners, has been performed to obtain the models to be tested precisely under the same conditions of the experiment Section 4.5. It is worth observing that using two different types of evaluations allows us to make sure that our results are not influenced by the type of experiment. The results of the experiments are shown in Tables 9 and 10. The statistical tests aim to verify if there are significant differences among the means of each metric: percent correct, area under the ROC curve and serialized model size. We proceeded as follows. First, we checked normality and sphericity of each sample by means of the Shapiro-Wilk normality test. Then, if normality and sphericity conditions were met, we applied one way repeated measures ANOVA; otherwise, we applied the Friedman test. In the latter case, when statistically significant differences were detected, we applied the Nemenyi post-hoc test to locate where these differences were. Tables A1-A12 in Appendix A show the results of the performed tests for the Breast Cancer dataset for each of the three metrics, and Tables A13-A24 in Appendix B show the results for the Monk's Problem 2 dataset.

Additional Experiments
Finally, we show the results of the evaluation with 10-fold cross-validation for Monk's problem 2 dataset and for the following four other datasets:

4.
Nursery dataset, with 8 input attributes, 12,960 instances, and 5 output classes (Table 14).    We have used the ENORA algorithm together with the ACC D and RMSE D objective functions in this case because these combinations have produced the best results for the Breast Cancer and Monk's problem 2 datasets evaluated in 10-fold cross-validation (population size equal to 50, 20,000 generations and M max = 10 + number of classes). Table 15 shows the results of the best combination ENORA-ACC or ENORA-RMSE together with the results of the classical rule-based classifiers.

Analysis of Results and Discussion
The results of our tests allow for several considerations. The first interesting observation is that NSGA-II identifies fewer solutions than ENORA on the Pareto front, which implies less diversity and therefore a worse hypervolume ratio, as shown in Figures 3 and 4. This is not surprising: in several other occasions [19,34,60], it has been shown that ENORA maintains a higher diversity in the population than other well-known evolutionary algorithms, with generally positive influence on the final results. Comparing the results in full training mode against the results in cross-validation or in splitting mode makes it evident that our solution produces classification models that are more resilient to over-fitting. For example, the classifier learned by PART with Monk's Problem 2 presents a 94.01% accuracy in full training mode that drops to 73.51% in splitting mode. A similar, although with a more contained drop in accuracy, is shown by the classifier learned with Breast Cancer dataset; at the same time, the classifier learned by ENORA driven by accuracy shows only a 5.57% drop in one case, and even an improvement in the other case (see Tables 5, 6, 9, and 10). This phenomenon is easily explained by looking at the number of rules: the more rules in a classifier, the higher the risk of over-fitting; PART produces very accurate classifiers, but at the price of adding many rules, which not only affects the interpretability of the model but also its resilience to over-fitting. Full training results seem to indicate that when the optimization model is driven by RMSE the classifiers are more accurate; nevertheless, they are also more prone to over-fitting, indicating that, on average, the optimization models driven by the accuracy are preferable.
From the statistical tests (whose results are shown in the Appendixes A and B) we conclude that among the six variants of the proposed optimization model there are no statistical significative differences, which suggests that the advantages of our method do not depend directly on a specific evolutionary algorithm or on the specific performance measure that is used to drive the evolutions. Significant statistical differences between our method and very simple classical methods such as OneR were expectable. Significant statistical differences between our method and a well-consolidated one such as PART have not been found, but the price to be paid for using PART in order to have similar results to ours is a very high number of rules (15 vs. 2 in one case and 47 vs. 7 in the other case).
We would like to highlight that both the Breast Cancer dataset and the Monk's problem 2 dataset are difficult to approximate with interpretable classifiers and that none of the analyzed classifiers obtains high accuracy rates using the cross-validation technique. Even powerful black-box classifiers, such as Random Forest and Logistic, obtain success rates below 70% in 10-fold cross-validation for these datasets. However, ENORA obtains a better balance (trade-off) between precision and interpretability than the rest of the classifiers. For the rest of the analyzed datasets, the accuracy obtained using ENORA is substantially higher. For example, for the Tic-Tac-Toe-Endgame dataset, ENORA obtains a 98.3299% success percentage with only two rules in cross-validation, while PART obtains 94.2589% with 49 rules, and JRip obtains 97.8079% with nine rules. With respect to the results obtained in the datasets Car, kr-vs-kp and Nursery, we want to comment that better success percentage can be obtained if the maximum number of evaluations is increased. However, better success percentages imply a greater number of rules, which is to the detriment of the interpretability of the models.

Conclusions and Future Works
In this paper, we have proposed a novel technique for categorical classifier learning. Our proposal is based on defining the problem of learning a classifier as a multi-objective optimization problem, and solving it by suitably adapting an evolutionary algorithm to this task; our two objectives are minimizing the number of rules (for a better interpretability of the classifier) and maximizing a metric of performance. Depending on the particular metric that is chosen, (slightly) different optimization models arise. We have tested our proposal, in a first instance, on two different publicly available datasets, Breast Cancer (in which each instance represents a patient that has suffered from breast cancer and is described by nine attributes, and the class to be predicted represents the fact that the patient has suffered a recurring event) and Monk's Problem 2 (which is an artificial, well-known dataset in which the class to be predicted represents a logical function), using two different evolutionary algorithms, namely ENORA and NSGA-II, and three different choices as a performance metric, i.e., accuracy, the area under the ROC curve, and the root mean square error. Additionally, we have shown the results of the evaluation in 10-fold cross-validation of the publicly available Tic-Tac-Toe-Endgame, Car, kr-vs-kp and Nursery datasets.
Our initial motivation was to design a classifier learning system that produces interpretable, yet accurate, classifiers: since interpretability is a direct function of the number of rules, we conclude that such an objective has been achieved. As an aside, observe that our approach allows the user to decide, beforehand, a maximum number of rules; this can also be done in PART and JRip, but only indirectly. Finally, the idea underlying our approach is that multiple classifiers are explored at the same time in the same execution, and this allows us to choose the best compromise between the performance and the interpretability of a classifier a posteriori.
As a future work, we envisage that our methodology can benefit from an embedded future selection mechanism. In fact, all attributes are (ideally) used in every rule of a classifier learned by our optimization model. By simply relaxing such a constraint, and by suitably re-defining the first objective in the optimization model (e.g., by minimizing the sum of the lengths of all rules, or similar measures), the resulting classifiers will naturally present rules that use more features as well as rules that use less (clearly, the implementation must be adapted to obtain an initial population in which the classifiers have rules of different lengths as well as mutation operators that allow a rule to grow or to shrink). Although this approach does not follow the classical definition of feature selection mechanisms (in which a subset of features is selected that reduces the dataset over which a classifier is learned), it is natural to imagine that it may produce even more accurate classifiers, and more interpretable at the same time.
Currently, we are implementing our own version of multi-objective differential evolution (MODE) for rule-based classification for inclusion in the Weka Open Source Software issued under the GNU General Public License. The implementation of other algorithms, such as MOEA/D, their adaptation in the Weka development platform and subsequent analysis and comparison are planned for future work. Acknowledgments: This study was partially supported by computing facilities of Extremadura Research Centre for Advanced Technologies (CETA-CIEMAT), funded by the European Regional Development Fund (ERDF). CETA-CIEMAT belongs to CIEMAT and the Government of Spain.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:  Table A3. Nemenyi post-hoc procedure for percent correct metric-Breast Cancer dataset.  Table A4. Summary of statistically significant differences for percent correct metric-Breast Cancer dataset.