An Ensemble of Bayesian Networks for Multilabel Classification

1 Introduction

In traditional classification each instance belongs only to a single class. Multilabel classification generalizes this idea by allowing each instance to belong to different relevant classes at the same time. A simple approach to deal with multilabel classification is binary relevance (BR): a binary classifier is independently trained for each class and then predicts whether such class is relevant or not for each instance. BR ignores dependencies among the different classes; yet it can be quite competitive [Dembczynski et al., 2012] especially under loss functions which decompose label-wise, such as the Hamming loss.

Instead the global accuracy (also called exact match) does not decompose label-wise; a classification is regarded as accurate only if the relevance of every class has been correctly identified. To obtain good performance under global accuracy, it is necessary identifying the maximum a posteriori (MAP) configuration, that is, the mode of the joint probability distribution of the classes given the features, which in turn requires to properly model the stochastic dependencies among the different classes. A state-of-the-art approach to this end is the classifier chain (CC) [Read et al., 2011]. Although CC has not been designed to minimize a specific loss function, it has been recently pointed out [Dembczynski et al., 2010; 2012] that it can be interpreted as a greedy approach for identifying the mode of the joint distribution of the classes given the value of the features; this might explain its good performance under global accuracy.

∗Work supported by the Swiss NSF grant no. 200020 134759 / 1 and by the Hasler foundation grant n. 10030.

Bayesian networks (BNs) are a powerful tool to learn and perform inference on the joint distribution of several variables; yet, they are not commonly used in multilabel classification. They pose two main problems when dealing with multilabel classification: learning from data the structure of the BN model and performing inference on the learned model in order to issue a classification. In [Bielza et al., 2011], the structural learning problem is addressed by partitioning the arcs of the structure into three sets: links among the class variables (class graph), links among features (features graph) and links between class and features variables (bridge graph). Each subgraph is separately learned and can have different topology (tree, polytree, etc). Inference is performed either by an optimized enumerative scheme or by a greedy search, depending on the number of classes.

Thoroughly searching for the best structure introduces the issue of model uncertainty: several structures, among the many analyzed, might achieve similar scores and thus the model selection can become too sensitive on the specificities of the training data, eventually increasing the variance component of the error. In traditional classification, this problem has been successfully addressed by instantiating all the BN classifiers whose structure satisfy some constraints and then combining their inferences, avoiding thus an exhaustive search for the best structure [Webb et al., 2005]. The resulting algorithm, called AODE, is indeed known to be a highperformance classifier.

In this paper, we extend to the multilabel case the idea of averaging over a constrained family of classifiers. We assume the graph connecting the classes to be a tree. Rather than searching for the optimal tree, we instantiate a different classifier for each class; each classifier adopts a different class as root. We moreover assume the independence of the features given the classes, thus generalizing to the multilabel case the naive Bayes assumption. This implies that the subgraph connecting the features is empty. Thanks to the naivemultilabel assumption, we can identify the optimal subgraph linking classes and the features in polynomial time [de Campos and Ji, 2011]. This optimal subgraph is the same for all the models of the ensemble, which speeds up inference and learning. Summing up: we perform no search for the class graph, as we take the ensemble approach; we perform no search for the feature graph; we perform an optimal search for the bridge graph.

Eventually, we combine the inferences produced by the members of the ensemble via geometric average (logarithmic opinion pool), which minimizes the average KL-divergence between the true distribution and the elements of an ensemble of probabilistic models [Heskes, 1997]. It is moreover well known that simple combination schemes (such as arithmetic average or geometric average) often dominates more refined combination schemes [Yang et al., 2007; Timmermann, 2006].

The inferences are tailored to the loss function being used. To minimize Hamming loss, we query each classifier of the ensemble about the posterior marginal distribution of each class variable, and then pool such marginals to obtain an ensemble probability of each class being relevant. To instead maximize the global accuracy, we compute the MAP inference in the pooled model. Although both the marginal and joint inferences are NP-hard in general, they can often be performed quite efficiently by state-of-the-art algorithms.

2 Probabilistic Multilabel Classification

We denote the array of class events as C = (C1, . . . ,Cn); this is an array of Boolean variables which expresses whether each class label is relevant or not for the given instance. Thus, C takes its values in {0, 1}n . A instantiation of class events is denoted as c = (c1, . . . , cn). We denote the set of features as F = (F1, . . . , Fm). An instantiation of the features is denoted as f = (f1, . . . , fm). A set of complete training instances D = {(c, f )} is supposed to be available.

For a generic instance, we denote as ˆc and ˜c respectively the labels assigned by the classifier and the actual labels; similarly, with reference to class ci, we denote as ˆ ci and ˜ c∗i its predicted and actual labels. The most common metrics to evaluate multilabel classifiers are the global accuracy (a.k.a. exact match), denoted as acc, and the Hamming loss (HL). Rather than HL, we use in the following the Hamming accuracy, defined as Hacc=1-HL. This simplifies the interpretation of the results as both acc and Hacc are better when they are higher. The metrics are defined as follows:

A probabilistic multilabel classifier classifies data according to a conditional joint probability distribution P(C|F). Our approach consists in learning an ensemble of probabilistic classifiers that jointly approximate the (unknown) true distribution P (C|F). When classifying an instance with features f , two different kinds of inference are performed depending on whether the goal is maximizing global accuracy or Hamming accuracy. The most probable joint class ˆc maximizes global accuracy (under the assumption that P(C|F) is the true distribution), and is given by

The maximizer of Hamming accuracy instead consists in selecting the most probable marginal class label separately for each Ci ∈ C:

with P(ci|f) = ∑C\{Ci} P(c|f). The choice of the appropriate inference to be performed given a specific loss function is further discussed in [Dembczynski et al., 2012].

3 An Ensemble of Bayesian Classifiers

As with traditional classification problems, building a multilabel classifier involves a learning phase, where the training data are used to tune parameters and perform model selection, and a test phase, where the learned model (or in our case, the ensemble of models) is used to classify instances.

Our approach consists in learning an ensemble of Bayesian networks from data, and then merging these models via geometric average (logarithmic pooling) to produce classifications. Bayesian networks [Koller and Friedman, 2009] provide a compact specification for the joint probability mass function of class and feature variables. Our ensemble contains n different Bayesian networks, one for each class event Ci (i = 1, . . . , n), each providing an alternative model of the joint distribution P (C, F).

3.1 Learning

The networks in the ensemble are indexed by the class events. With no lack of generality, we describe here the network associated to class Ci. First, a directed acyclic graph Gi over C such that Ci is the unique parent of all the other class events, and no more arcs are present, is considered (see Fig. 1.a). Following the Markov condition, 1 this corresponds to assume that, given Ci, the other classes are independent.

The features are assumed to be independent given the joint class C. This extends the naive Bayes assumption to the multilabel case. At the graph level, the assumption can be accomodated by augmenting the graph Gi, defined over C, to a graph Gi defined over the whole set of variables (C, F) as follows. The features are leaf nodes (i.e., they have no children) and their parents are, for each feature, all the class events (see Fig. 1.b).

A Bayesian network associated to Gi has maximum indegree (i.e., maximum number of parents) equal to the number of classes n: there are therefore conditional probability tables whose number of elements is exponential in n. This is not practical for problems with many classes. A structural learning algorithm can be therefore adopted to reduce the maximum in-degree. In practice, to remove some classto-feature arcs in Gi, we evaluate the output of the following optimization:

1 The Markov condition for directed graphs specifies the graph semantics in Bayesian networks: any variable is conditionally independent of its non-descendant non-parents given its parents.

where set inclusions among graphs should be intended in the arcs space, and

where Pa(Fj) denotes the parents of Fj according to G, and similarly Pa(Ci) are the parents of Ci. Moreover, ψα is BDEu score [Buntine, 1991] with equivalent sample size α. For instance, the score ψα[Fj,Pa(Fj)] is

where the second sum is over the possible states of Fj, and the first over the joint states of its parents. Moreover, Njik is the number of records such that Fj is in its k-th state and its parents in their i-th configuration, while Nji = ∑k Njik. Finally, αji is equal to α divided by the number of states of Fj and by the number of (joint) states of the parents of Fj, while αj = ∑i αji.

If we consider a network associated to a specific class event, the links connecting the class events are fixed; thus, the class variables have the same parents for all the graphs in the search space. This implies that the first sum in the right-hand side of (4) is constant. The optimization in (3) can be therefore achieved by considering only the features. Any subset of C is a candidate for the parents set of a feature, this reducing the problem to m independent local optimizations. In practice, the parents of Fj according to G∗ are

for each j = 1, . . . ,m. This is possible because of the bipartite separation between class events and features, while in general cases the graph maximizing all the local scores can include directed cycles. Moreover, optimization in (6) becomes more efficient by considering the pruning techniques proposed in [de Campos and Ji, 2011] (see Section 3.3).

According to Equation (6), the subgraph linking the classes does not impact on the outcome of optimization; in other words, all models of the ensemble have the same class-tofeature arcs.

Graph G∗ i is the qualitative basis for the specification of the Bayesian network associated to Ci. Given the data, the network probabilities are indeed specified following a standard Bayesian learning procedure. We denote by Pi(C,F) the corresponding joint probability mass function indexed by Ci. The graph induces the following factorization:

with the values of fj, ci, cl, cFj consistent with those of c, f . Note that according to the above model assumptions, the data likelihood term

is the same for all models in the ensemble.

(c) Removing class-to-feature arcs from

G1 to obtain G∗

1 .

Figure 1: The construction of the graphs of the Bayesian networks in the ensemble. The step from (b) to (c) is based on the structural learning algorithm in (6).

3.2 Classification

The models in the ensemble are combined to issue a classification according to the evaluation metric being used, in agreement with Equations (1) and (2). For global accuracy we compute the mode of the geometric average of joint distribution of each Bayesian network in the ensemble, given by

The above inference is equivalent to the computation of the MAP configuration in the Markov random field over variables C with potentials φi,l(ci, cl) = Pi(cl|ci), φcFj (cFj ) = P(fj|cFj ) and φi(ci) = Pi(ci), i, l = 1, . . . , n (l ̸= i) and j = 1, . . . ,m. Such an inference can be solved exactly by junction tree algorithms when the treewidth of the corresponding Bayesian network is sufficiently small (e.g., up to 20 in standard computers) or transformed into a mixedinteger linear program and subsequently solved by standard solvers or message-passing methods [Koller and Friedman,

2009; Sontag et al., 2008]. For very large models (n > 1000), the problem can be solved approximately by message-passing algorithms such as TRW-S and MPLP [Kolmogorov, 2006; Globerson and Jaakkola, 2007], or yet by network flow algorithms such as extended roof duality [Rother et al., 2007].

For Hamming accuracy, we instead classify an instance f as ˆc = (ˆ c1, . . . , ˆ cn) such that for each l = 1, . . . , n

where

is the marginal probability of class event cl according to model Pi(c, f ).

The computation of the marginal probability distributions in (11) can also be performed by junction tree methods, when the treewidth of the corresponding Bayesian network is small, or approximated by message-passing algorithms such as belief propagation and its many variants [Koller and Friedman, 2009].

Overall we have defined two types of inference for multilabel classification, based on the models in (9) and (10), which will be called, respectively, EBN-J and EBN-M (these acronyms standing for ensemble of Bayesian networks with the joint and the marginal).

A linear, instead of log-linear, pooling could be considered to average the models. The complexity of marginal inferences is equivalent whether we use linear or log-linear pooling, and the same algorithms can be used. The theoretical complexity of MAP inference is also the same for either pooling method (their decision versions are both NP-complete, see [Park and Darwiche, 2004]), but linear pooling requires different algorithms to cope with the presence of a latent variable [Mau´a and de Campos, 2012]. In our experiments, log-linear pooling produced slightly more accurate classifications and we report only these results due to the limited space.

3.3 Computational Complexity

Let us evaluate the computational complexity of the algorithms EBN-J and EBN-M presented in the previous section. We distinguish between learning and classification complexity, where the latter refers to the classification of a single instantiation of the features. Both space and time required for computations are evaluated. The orders of magnitude of these descriptors are reported as a function of the (training) dataset size k, the number of classes n, the number of features m, the average number of states for the features f = m-1 ∑mj=1 |Fj|, and the maximum in-degree g = maxj=1,...,m ‖CFj ‖ (where ‖ · ‖ returns the number elements in a joint variable). We first consider a single network in the ensemble. Regarding space, the tables P(Cl), P(Cl|Ci), P(Fj|CFj ) need, respectively 2, 4 and |Fi| · 2‖CFi ‖ . We already noticed that the tables associated to the features are the same for each network. Thus, when considering the whole ensemble only the (constant) terms associated to the classes are multiplied by n and the overall space complexity is O(4n+f2g ). These tables should be available during both learning and classification for both classifiers.

Regarding time, let us start from the learning. Structural learning requires O(2gm), while network quantification require to scan the dataset and learn the probabilities, i.e., for the whole ensemble, O(nk). For classification, both the computation of a marginal in Bayesian network and the MAP inference in the Markov random field obtained by aggregating the models can be solved exactly by junction tree algorithms in time exponential in the network treewidth, or approximated reasonably well by message-passing algorithms in time O(ng2g ). Since, as proved by [de Campos and Ji, 2011], g = O(log k), the overall complexity is polynomial (approximate inference algorithms are only used when the treewidth is too high).

3.4 Implementation

We have implemented the high-level part of the algorithm in Python. For structural learning, we modified the GOBNILP package 2 in order to consider the constraints in (3) during the search in the graphs space. The inferences necessary for classification have been performed using the junction tree and belief propagation algorithms implemented in libDAI 3 , a library for inference in probabilistic graphical models.

4 Experiments

We compare our models against the binary relevance method, using naive Bayes as base classifiers (BR-NB); the ensemble of chain classifier using Naive Bayes (ECC-NB) and J48 (ECC-J48) as base classifier. ECC stands for 'ensemble of chain classifiers', where the ensemble contains a number of models which equals the number of classes. We set to 5 the number of members of the ensemble, namely the number of chains which are implemented using different random label orders. Therefore, ECC runs 5n base models, where n is the number of classes. We use the implementation of these methods provided by MEKA package. 4 It has not been possible comparing also with the Bayesian networks models of [Bielza et al., 2011] or with the conditional random fields in [Ghamrawi and McCallum, 2005], because of the lack of public domain implementations.

Regarding the parameters of our model, we set the equivalent sample size for the structural learning to α = 5. Moreover, we applied a bagging procedure to our ensemble, generating 5 different bootstrap replicates of the original training set. In this way, also our approach run 5n base models.

On each bootstrap replicate, the ensemble is learned from scratch. To combine the output of the ensemble learned on the different replicates, for EBN-J, we assign to each class the value appearing in the majority of the outputs (note that the number of bootstraps is odd). For EBN-M, instead, we take the arithmetic average of the marginals of the different bootstrapped models.

On each data set we perform feature selection using the correlation-based feature selection (CFS) [Witten et al., 2011, Chap. 7.1], which is a method developed for traditional classification. We perform CFS n times, once for each different class variables. Eventually, we retain the union of the features selected by the different runs of CFS.

2 http://www.cs.york.ac.uk/aig/sw/gobnilp

3 http://www.libdai.org

4 http://meka.sourceforge.net

The characteristics of the data sets are in Table 1. The density is the average number of classes in the true state over the total number of classes. We perform 10-fold cross validation, stratifying the folds according to the most rare class. For Slashdot, given the high number instances, we validate the models by a 70-30% split.

Since Bayesian networks need discrete features, we discretize each numerical feature into four bins. The bins are delimited by the 25-th, the 50-th and 75-th percentile of the value of the feature.

The results for the different data sets are provided in Tables 2-7. We report both global accuracy (acc) and Hamming accuracy (Hacc). Regarding our ensembles, we evaluate EBN-J in terms of global accuracy and EBN-M in terms of Hamming accuracy. Thus, we perform a different inference depending on the accuracy function which has to be maximized. Empirical tests confirmed that, as expected (see the discussion in Section 3.2), EBN-J has higher global accuracy and lower Hamming accuracy than EBN-M. In the special case of the Cal500 data set, the very high number of classes makes unnecessary the evaluation of the global accuracy which is almost zero for all the algorithms. Only Hamming accuracies are therefore reported.

Regarding computation times, EBN-J and EBN-M are one or two orders of magnitude slower than the MEKA implementations of ECC-J48/NB (which, in turn, are slower than BR-NB). This gap is partially due to our implementation, which at this stage is only prototypal.

We summarize the results in Table 8, which provides the mean rank of the various classifiers across the various data sets. For each indicator of performance, the best-performing approach is boldfaced. Although we cannot claim statistical significance of the difference among ranks (more data sets need to be analyzed to this purpose), the proposed ensemble achieves the best rank on both metrics; it is clearly competitive with state-of-the art approaches.

5 Conclusions

A novel approach to multilabel classification with Bayesian networks is proposed, based on an ensemble of models averaged by logarithmic opinion pool. Two different inference algorithm based respectively on the MAP solution on the joint model (EBN-J) and on the average of the marginals (EBN-M)

are proposed, respectively to maximize the exact match and the Hamming accuracy. Empirical validation shows competing results with the state of the art. As future work we intend to support missing data both in the learning and in the classification step. Also more sophisticated mixtures could be considered, for instance by learning the maximum a posteriori weights for the mixture of the models.

Acknowledgments

We thank Cassio P. de Campos and David Huber for valuable discussions.