Nearest labelset using double distances for multi-label classification

Hypercube view of a multi-label problem

In multi-label classification, we are given a set of possible output labels $\mathcal{L}=\left\{1,2,\dots ,L\right\}$. Each instance with a feature vector $\overrightarrow{x}\in {\mathbb{R}}^d$ is associated with a subset of these labels. Equivalently, the subset can be described as $\overrightarrow{y}=\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(L\right)}\right)$, where y⁽ⁱ⁾ = 1 if label i is associated with the instance and y⁽ⁱ⁾ = 0 otherwise. A multi-label training data set is described as $T=\left\{\left({\overrightarrow{x}}_i,{\overrightarrow{y}}_i\right)\right\},$i=1 , {2, …, N}.

Any labelset $\overrightarrow{y}$ can be described as a vertex in the L-dimensional unit hypercube (Tai & Lin, 2012). Each component y⁽ⁱ⁾ of $\overrightarrow{y}$ represents an axis of the hypercube. As an example, Fig. 1 illustrates the label space of a multi-label problem with three labels (y⁽¹⁾, y⁽²⁾, y⁽³⁾).

Assume that the presence or absence of each label is modeled independently with a probabilistic classifier. For a new instance, the classifiers provide the probabilities, p⁽¹⁾, …, p^(L), that the corresponding labels are associated with the instance. Using the probability outputs, we may obtain a L-dimensional vector $\hat{\overrightarrow{p}}=\left(p^{\left(1\right)},p^{\left(2\right)},\dots ,p^{\left(L\right)}\right)$. Every element of $\hat{\overrightarrow{p}}$ has a value from 0 to 1 and the vector $\hat{\overrightarrow{p}}$ is an inner point in the hypercube (see Fig. 1). Given $\hat{\overrightarrow{p}}$ the prediction task is completed by assigning the inner point to a vertex of the cube.

For the new instance, we may calculate the Euclidean distance, $D_{{\overrightarrow{y}}_i}$, between $\hat{\overrightarrow{p}}$ and each ${\overrightarrow{y}}_i$ (i.e., the labelset of the ith training instance). In Fig. 1, three training instances ${\overrightarrow{y}}_1$, ${\overrightarrow{y}}_2$ and ${\overrightarrow{y}}_3$ and the corresponding distances are shown. A small distance $D_{{\overrightarrow{y}}_i}$ indicates that ${\overrightarrow{y}}_i$ is likely to be the labelset for the new instance.

Nearest labelset using double distances (NLDD)

In addition to computing the distance in the label space, $D_{{\overrightarrow{y}}_i}$, we may also obtain the (Euclidean) distance in the feature space, denoted by $D_{{\overrightarrow{x}}_i}$. The proposed method, NLDD, uses both $D_{\overrightarrow{x}}$ and $D_{\overrightarrow{y}}$ as predictors to find a training labelset that minimizes the expected loss. For each test instance, we define loss as the number of misclassified labels out of L labels. The expected loss is then Lθ where $\theta =g\left(D_{\overrightarrow{x}},D_{\overrightarrow{y}}\right)$ represents the probability of misclassifying each label. The predicted labelset, ${\hat{\overrightarrow{y}}}^{\ast }$, is the labelset observed in the training data that minimizes the expected loss: (1)${\displaystyle {\hat{\overrightarrow{y}}}^{\ast }=\underset{\overrightarrow{y}\in T}{argmin}Lg\left(D_{\overrightarrow{x}},D_{\overrightarrow{y}}\right)}$ The loss follows a binomial distribution with L trials and a parameter θ. We model $\theta =g\left(D_{\overrightarrow{x}},D_{\overrightarrow{y}}\right)$ using binomial regression. Specifically, (2)${\displaystyle log\left(\frac{\theta }{1-\theta }\right)={\beta }_0+{\beta }_1{D}_{\overrightarrow{x}}+{\beta }_2{D}_{\overrightarrow{y}}}$ where β₀, β₁ and β₂ are the model parameters. Greater values for β₁ and β₂ imply that θ becomes more sensitive to the distances in the feature and label spaces, respectively. The misclassification probability decreases as $D_{\overrightarrow{x}}$ and $D_{\overrightarrow{y}}$ approach zero.

A test instance with $D_{\overrightarrow{x}}=D_{\overrightarrow{y}}=0$ has a duplicate instance in the training data (i.e., with identical features). The predicted probabilities for the test instance are either 0 or 1 and the match the labels of the duplicate training observation. For such a “double”-duplicate instance (i.e., $D_{\overrightarrow{x}}=D_{\overrightarrow{y}}=0$), the probability of misclassification is 1∕(1 + e^−β₀) > 0. As expected, the uncertainty of a test observation with a “double-duplicate” training observation is greater than zero. This is not surprising: duplicate training observations do not necessarily have the same response, and neither do double-duplicate observations.

The model in Eq. (2) implies $g\left(D_{\overrightarrow{x}},D_{\overrightarrow{y}}\right)=1∕\left(1+e^{-\left({\beta }_0+{\beta }_1{D}_{\overrightarrow{x}}+{\beta }_2{D}_{\overrightarrow{y}}\right)}\right)$. Because $log\left(\frac{\theta }{1-\theta }\right)$ is a monotone transformation of θ and L is a constant, the minimization problem in (1) is equivalent to (3)${\displaystyle {\hat{\overrightarrow{y}}}^{\ast }=\underset{\left(\overrightarrow{x},\overrightarrow{y}\right)\in T}{argmin}{\beta }_1{D}_{\overrightarrow{x}}+{\beta }_2{D}_{\overrightarrow{y}}}$ That is, NLDD predicts by choosing the labelset of the training instance that minimizes the weighted sum of the distances. For prediction, the only remaining issue is how to estimate the weights.

 
________________________________________________________________________________________________________________________________________________________________ 
Algorithm 1 The training process of NLDD 
________________________________________________________________________________________________________________________________________________________________ 
  Input: training data T, number of labels L 
  Output: probabilistic classifiers h(i), binomial regression g 
  Split T into T1 and T2 
  for i = 1 to L do 
      train probabilistic classifier h(i) based on T 
      train probabilistic classifier h(i)∗ based on T1 
  end for 
  S,W ←∅ 
  for each instance in T2 do 
      obtain ˆ p = (h(1)∗ (x),...,h(L)∗ (x)) 
      for each instance in T1 do 
          compute Dx and Dy 
          W ← W ∪ (Dx,Dy) 
      end for 
      find m1,m2 ∈ W 
      update S ← S ∪{m1,m2} 
  end for 
  Fit log (θ __1−θ) = β0 + β1Dx + β2Dy to S 
  Obtain g : S → ˆ θ =   e ˆf1+e ˆf where  ˆ f= ˆ β0 + ˆ β1Dx + ˆ β2Dy 
________________________________________________________________________________________________________________________________________________________________    

 
_______________________________________________________________________________________________________ 
Algorithm 2 The classification process of NLDD 
_______________________________________________________________________________________________________ 
Input: new instance x, binomial model g, probabilistic classifiers h(i), training data T of 
size N 
Output: multi-label classification vector ˆ y 
for j = 1 to N do 
    compute ˆ p = (h(1)(x),...,h(L)(x)) 
    compute Dxj  and Dyj 
    obtain ˆ θj ← g(Dxj,Dyj) 
end for 
return ˆ y ← argmin  yj∈T   ˆ θj 
________________________________________________________________________________________________________ 

Estimating the relative weights of the two distances

The weights β₀, β₁ and β₂ can be estimated using binomial regression. Binomial regression can be fit by running separate logistic regressions, one for each of the L labels. To run the regressions $D_{\overrightarrow{x}}$ and $D_{\overrightarrow{y}}$ need to be computed on the training data. For this purpose we split the training data (T) equally into an internal training data set, T₁, and an internal validation data set, T₂.¹ We next fit a binary classifier to each of the L labels separately and obtain the labelset predictions (i.e., probability outcomes) for the instances in T₂. In principle, each observation in T₂ can be paired with each observation in T₁, creating a $\left(D_{\overrightarrow{x}},D_{\overrightarrow{y}}\right)$ pair, and the regressions can be run on all possible pairs. Note that matching any single instance in T₂ to those in T₁ results in N∕2 distance pairs. However, most of the pairs are uninformative because the distance in either the feature space or the label space is very large. Since candidate labelsets for the final prediction will have a small $D_{\overrightarrow{y}}$ and a small $D_{\overrightarrow{x}}$, it is reasonable to focus more on the behavior of the loss especially at small values of $D_{\overrightarrow{x}}$ and $D_{\overrightarrow{y}}$ than considering the loss at the entire range of the distances. Moreover, since T₂ contains N∕2 instances, the number of possible pairs is potentially large (N²/4). Therefore, to reduce computational complexity, for each instance we only identify two pairs: the pair with the smallest distance in $\overrightarrow{x}$ and the pair with the smallest distance in $\overrightarrow{y}$. In case of ties in one distance, the pair with the smallest value in the other distance is chosen. More formally, we identify the first pair m_i1 by ${\displaystyle m_{i_1}=\underset{\left(D_x,D_y\right)\in W_{ix}}{argmin}{D}_y}$ where W_ix is the set of pairs that are tied; i.e., that each corresponds to the minimum distance in D_x. Similarly, the second pair m_i2 is found by ${\displaystyle m_{i_2}=\underset{\left(D_x,D_y\right)\in W_{iy}}{argmin}{D}_x.}$ where W_iy is the set of labels that are tied with the minimal distance in D_y. Figure 2 illustrates an example of how to identify m_i1 and m_i2 for N = 20. Our goal was to identify the instance with the smallest distance in $\overrightarrow{x}$ and the instance with the smallest distance in $\overrightarrow{y}$. Note that m_i1 and m_i2 may be the same instance. If we find a single instance that minimizes both distances, we use just that instance. (A possible duplication of that instance is unlikely to make any difference in practice).

The two pairs corresponding to the ith instance in T₂ are denoted as the set $S_i=\left\{m_{i_1},m_{i_2}\right\}$, and their union for all instances is denoted as $S={\bigcup }_{i=1}^{N∕2}{S}_i$. The binomial regression specified in Eq. (2) is performed on the $2\frac{N}{2}=N$ instances in S. Algorithm 1 outlines the training procedure.

For the classification of a new instance, we first obtain $\hat{\overrightarrow{p}}$ using the probabilistic classifiers fitted to the training data T. $D_{{\overrightarrow{x}}_j}$ and $D_{{\overrightarrow{y}}_j}$ are obtained by matching the instance with the jth training instance. Using the MLEs ${\hat{\beta }}_0$, ${\hat{\beta }}_1$ and ${\hat{\beta }}_2$, we calculate ${\hat{\theta }}_j=\frac{e^{{\hat{f}}_j}}{1+e^{{\hat{f}}_j}}$ where ${\hat{f}}_j={\hat{\beta }}_0+{\hat{\beta }}_1{D}_{{\overrightarrow{x}}_j}+{\hat{\beta }}_2{D}_{{\overrightarrow{y}}_j}$. The final prediction of the new instance is obtained by ${\displaystyle \hat{\overrightarrow{y}}=\underset{{\overrightarrow{y}}_j\in T}{argmin}\hat{E}\left(loss\right)=\underset{{\overrightarrow{y}}_j\in T}{argmin}{\hat{\theta }}_j.}$ The second equality holds because $\hat{E}\left(loss\right)=L\hat{\theta }$ and L is a constant. As in LP, NLDD chooses a training labelset as the predicted vector. Algorithm 2 outlines the classification procedure.

The training time of NLDD is O(L(f(d, N) + f(d, N∕2) + g(d, N∕2)) + N²(d + L) + Nlog(k)) where O(f(d, N)) is the complexity of each binary classifier with d features and N training instances, O(g(d, N∕2)) is the complexity for predicting each label for T₂, N²(d + L) is the complexity for obtaining the distance pairs for the regression and O(Nlog(k)) is the complexity for fitting a binomial regression. T₁ and T₂ have N∕2 instances respectively. O(Lf(d, N∕2)) is the complexity for fitting binary classifiers using T₁ and obtaining the probability results for T₂ takes O(Lg(d, N∕2)). For each instance of T₂, we obtain N∕2 numbers of distance pairs. This has complexity O((N∕2)(d + L)). Since there are N∕2 instances, overall it takes O((N∕2)(N∕2)(d + L)) or O(N²(d + L)) when omitting the constant. Among the N∕2 pairs for each instance of T₂, we only identify at most 2 pairs. This implies N∕2 ≤ s ≤ N where s is the number of elements in S. Each iteration of the Newton–Raphson method has a complexity of O(N). For k-digit precision complexity O(logk) is required (Ypma, 1995). Combined, the complexity for estimating the parameters with k-digit precision is O(Nlog(k)). In practice, however, this term is dominated by N²(d + L) as we can set k <  < N.

Data sets

We evaluated the proposed approach using nine commonly used multi-label data sets from different domains. Table 1 shows basic statistics for each data set including its domain, numbers of labels and features. In the text data sets, all features are categorical (i.e., binary). The last column “lcard”, short for label cardinality, represents the average number of labels associated with an instance. The data sets are ordered by (|L|⋅|X|⋅|E|).

The emotions data set (Trohidis et al., 2008) consists of pieces of music with rhythmic and timbre features. Each instance is associated with up to 6 emotion labels such as “sad-lonely”, “amazed-surprised” and “happy-pleased”. The scene data set (Boutell et al., 2004) consists of images with 294 visual features. Each image is associated with up to 6 labels including “mountain”, “urban” and “beach”. The yeast data set (Elisseeff & Weston, 2001) contains 2,417 yeast genes in the Yeast Saccharomyces Cerevisiae. Each gene is represented by 103 features and is associated with a subset of 14 functional labels. The medical data set consists of documents that describe patient symptom histories. The data were made available in the Medical Natural language Processing Challenge in 2007. Each document is associated with a set of 45 disease codes. The slashdot data set consists of 3,782 text instances with 22 labels obtained from Slashdot.org. The enron data set (Klimt & Yang, 2004) contains 1,702 email messages from the Enron corporation employees. The emails were categorized into 53 labels. The ohsumed data set (Hersh et al., 1994) is a collection of medical research articles from MEDLINE database. We used the same data set as in (Read et al., 2011) that contains 13,929 instances and 23 labels. The tmc2007 data set (Srivastava & Zane-Ulman, 2005) contains 28,596 aviation safety reports associated with up to 22 labels. Following Tsoumakas, Katakis & Vlahavas (2011), we used a reduced version of the data set with 500 features. The bibtex data set (Katakis, Tsoumakas & Vlahavas, 2008) consists of 7,395 bibtex entries for automated tag suggestion. The entries were classified into 159 labels. All data sets are available online at: MULAN (http://mulan.sourceforge.net/datasets-mlc.html) and MEKA (http://meka.sourceforge.net/#datasets).

Evaluation metrics

Multi-label classifiers can be evaluated with various loss functions. Here, four of the most popular criteria are used: Hammingloss, 0∕1loss, multi-label accuracy and F-measure. These criteria are defined in the following paragraphs.

Let L be the number of labels in a multi-label problem. For a particular test instance, let $\overrightarrow{y}=\left(y^{\left(1\right)},\dots ,y^{\left(L\right)}\right)$ be the labelset where y^(j) = 1 if the jth label is associated with the instance and 0 otherwise. Let $\hat{\overrightarrow{y}}=\left({\hat{y}}^{\left(1\right)},\dots ,{\hat{y}}^{\left(L\right)}\right)$ be the predicted values obtained by any machine learning method. Hammingloss refers to the percentage of incorrect labels. The Hammingloss for the instance is ${\displaystyle Hammingloss=1-\frac{1}{L}\sum _{j=1}^L\mathbb{1}\left\{y^{\left(j\right)}={\hat{y}}^{\left(j\right)}\right\}}$ where 1 is the indicator function. Despite its simplicity, the Hammingloss may be less discriminative than other metrics. In practice, an instance is usually associated with a small subset of labels. As the elements of the L-dimensional label vector are mostly zero, even the empty set (i.e., zero vector) prediction may lead to a decent Hammingloss.

The 0∕1loss is 0 if all predicted labels match the true labels and 1 otherwise. Hence, ${\displaystyle 0∕1loss=1-\mathbb{1}\left\{\overrightarrow{y}=\hat{\overrightarrow{y}}\right\}.}$ Compared to other evaluation metrics, 0∕1 loss is strict as all the L labels must match to the true ones simultaneously. The multi-label accuracy (Godbole & Sarawagi, 2004) (also known as the Jaccard index) is defined as the number of labels counted in the intersection of the predicted and true labelsets divided by the number of labels counted in the union of the labelsets. That is, ${\displaystyle Multi\text{-}labelaccuracy=\frac{|\overrightarrow{y}\cap \hat{\overrightarrow{y}}|}{|\overrightarrow{y}\cup \hat{\overrightarrow{y}}|}.}$ The multi-label accuracy measures the similarity between the true and predicted labelsets.

The F-measure is the harmonic mean of precision and recall. The F-measure is defined as ${\displaystyle F\text{-}measure=\frac{2|\overrightarrow{y}\cap \hat{\overrightarrow{y}}|}{\left|\overrightarrow{y}\right|+\left|\hat{\overrightarrow{y}}\right|}.}$

The metrics above were defined for a single instance. On each metric, the overall value for an entire test data set is obtained by averaging out the individual values.

Experimental setup

We compared our proposed method against BR, SMBR, ECC, RAKEL, HOMER, RF-PCT, MLKNN and CBM. To train multi-label classifiers, the parameters recommended by the authors were used, since they appear to give the best (or comparable to the best) performance in general. In the case of MLKNN, we set the number of neighbors and the smoothing parameter to 10 and 1 respectively. For RAKEL, we set the number of separate models to 2L and the size of each sub-labelset to 3. For ECC, the number of CC models for each ensemble was set to 10. For HOMER, the number of clusters was set to 3 as used in Liu et al. (2015). On the larger data sets (ohsumed, tmc2007 and bibtex), we fit ECC using reduced training data sets (75% of the instances and 50% of the features) as suggested in Read et al. (2011). On the same data sets, we ran NLDD using 70% of the training data to reduce redundancy in learning.

For NLDD, we used support vector machines (SVM) (Vapnik, 2000) as the base classifier on unscaled variables with a linear kernel and tuning parameter C = 1. The SVM scores were converted into probabilities using Platt’s method (Platt, 2000). The analysis was conducted in R (R Core Team, 2014) using the e1071 package (Meyer et al., 2014) for SVM. For the data sets with less than 5,000 instances 10-fold cross validations (CV) were performed. On the larger data sets, we used 75/25 train/test splits. For fitting binomial regression models, we divided the training data sets at random into two parts of equal sizes.

For RF-PCT, we used the Clus (http://clus.sourceforge.net) system. In the pre-pruning strategy of PCT, the significance level for the F-test was automatically chosen from {0.001, 0.005, 0.01, 0.05, 0.1, 0.125} using a reserved prune-set.

For CBM, we used the authors’ Java program (https://github.com/cheng-li/pyramid). The default settings (e.g., logistic regression and 10 iterations for the EM algorithm) were used on non-large data sets. For the large data sets tmc2007 and bibtex, the number of iterations was set to 5 and random feature reduction was applied as suggested by the developers. On each data set we used the train/test split recommended at their website (https://github.com/cheng-li/pyramid).

To test the hypothesis that all classifiers perform equally, we used the Friedman test as recommended by Demšar (2006). We then compared NLDD with each of the other methods using Wilcoxon signed-rank tests (Wilcoxon, 1945). We adjusted p-values for multiple testing using Hochberg’s method (Hochberg, 1988).

In NLDD, when calculating distances in the feature spaces we used the standardized features so that no particular features dominated distances. For a numerical feature variable x, the standardized variable z is obtained by $z=\left(x-\bar{x}\right)∕\text{sd}\left(x\right)$ where $\bar{x}$ and sd(x) are the mean and standard deviation of x in the training data.

Results

Tables 2 to 5 summarize the results in terms of Hammingloss, 0∕1loss, multi-label accuracy and F-measure, respectively. We also ranked the algorithms for each metric.

According to the Friedman tests, the classifiers are not all equal (p < 0.05). The post-hoc analysis - adjusted for multiple testing - showed that NLDD performed significantly better than BR and SMBR on all metrics, significantly better than RAKEL and MLKNN on all but Hammingloss, significantly better than HOMER on Hammingloss and 0∕1loss, and significantly better than ECC and RF-PCT on 0∕1loss. No method performed statistically significantly better than NLDD on any evaluation metric.

NLDD achieved lowest (i.e., best) average ranks on 0∕1loss and multi-label accuracy, while ECC and RF-PCT achieved the lowest average ranks on the F-measure and Hammingloss, respectively. On both F-measure and Hammingloss, NLDD achieved the second lowest (i.e., best) average ranks. CBM achieved the second lowest average rank on 0∕1loss and multi-label accuracy. The performance of CBM on the 0∕1loss was very variable achieving the lowest rank on five out of nine data sets and the second worst on two data sets.

We next look at the performance of NLDD by whether or not the true labelsets were observed in the training data. A labelset has been observed if the exact labelset can be found in the training data and unobserved otherwise. Since NLDD makes a prediction by choosing a training labelset, a predicted labelset can only be partially correct on an unobserved labelset. Table 6 compares the evaluation results of BR and NLDD on two separate subsets of the test set of the bibtex data. The bibtex data were chosen because the data set contains by far the largest percentage of unobserved labelsets (33%) among the data sets investigated. The test data set was split into subsets A and B; if the labelset of a test instance was an observed labelset, the instance was assigned to A; otherwise the instance was assigned to B. For all of the four metrics, NLDD outperformed BR even though 33% of the labelsets in the test data were unobserved labelsets.

We next look at the three regression parameters the proposed method (NLDD) estimated (Eq. (2)) for each data set in more detail. Table 7 displays the MLE of the parameters of the binomial model in each data set. In all data sets, the estimates of β₁ and β₂ were all positive. The positive slopes imply that the expected loss (or, equivalently the probability of misclassification for each label) decreases as $D_{\overrightarrow{x}}$ or $D_{\overrightarrow{y}}$ decreases.

From the values of ${\hat{\beta }}_0$ we may infer how low the expected loss is when either $D_{\overrightarrow{x}}$ or $D_{\overrightarrow{y}}$ is 0. For example, ${\hat{\beta }}_0=-3.5023$ in the scene data set. If $D_{\overrightarrow{x}}=0$ and $D_{\overrightarrow{y}}=0$, $\hat{p}=0.0292$ because $log\frac{\hat{p}}{1-\hat{p}}=-3.5023$. Hence $\hat{E}\left(loss\right)=L\hat{p}=6\cdot 0.0292=0.1752$. This is the expected number of mismatched labels for choosing a training labelset whose distances to the new instance are zero in both feature and label spaces. The results suggest the expected loss would be very small when classifying a new instance that had a duplicate in the training data ($D_{\overrightarrow{x}}=0$) and whose labels are predicted with probability 1 and the predicted labelset was observed in the training data ($D_{\overrightarrow{y}}=0$).