Bi-Parameter Space Partition for Cost-Sensitive SVM

1 Introduction

Ever since Vapnik's influential work in statistical learning theory [Vapnik and Vapnik, 1998], Support Vector Machines (SVMs) have been successfully applied to a lot of classification problems due to its good generalization performance. However, in many real-world classification problems such as medical diagnosis [Park et al., 2011], object recognition [Zhang and Zhou, 2010], business decision making [Cui et al., 2012], and so on, the costs of different types of mistakes are naturally unequal. Cost sensitive learning [Sheng and Ling, 2006] takes the unequal misclassification costs into consideration, which has also been deemed as a good solution to class-imbalance learning where the class distribution is highly imbalanced [Elkan, 2001]. There have been several cost-sensitive SVMs, such as the boundary movement [Shawe-Taylor, 1999], biased penalty (2C-SVM [Sch¨olkopf and Smola, 2002] and 2ν-SVM [Davenport et al., 2010]), cost sensitive hinge loss [Masnadi-Shirazi and Vasconcelos,

2010], and so on. In this paper, we focus on the most popular one (2C-SVM 1 ) of them.

Given a training set S = {(x1, y1), · · · , (xl, yl)} where xi ∈ Rd and yi ∈ {+1,-1}, 2C-SVM introduces two cost parameters C+ and C- to denote the costs of false negative and false positive respectively, and considers the following primal formulation:

where φ(xi) denotes a fixed feature-space transformation, S + = {(xi, yi) : yi = +1}, and S - = {(xi, yi) : yi = -1}. The dual problem of (1) is

where Q is a positive semidefinite matrix with Qij = yiyjK(xi, xj) and K(xi, xj) = 〈φ(xi), φ(xj)〉. It is obviously noted that how one tunes the cost parameter pair (C+, C-) to achieve optimal generalization performance (it is also called the problem of model selection) is a central problem of CS-SVM.

A general approach to tackle this problem is to specify some candidate parameter values, and then apply cross validation (CV) to select the best choices. A typical implementation for this approach is grid search [Mao et al., 2014]. However, extensive explorating the optimal parameter values is seldom pursued, because there exist double-sided difficulties. 1) It requires to train the classifier many times under different parameter settings. 2) And testing it on the validation dataset for each parameter setting.

To overcome the first difficulty, solution path algorithms were proposed for many learning models, such as C-SVM [Hastie et al., 2004], ϵ-SVR [Gunter and Zhu, 2007], quantile regression [Rosset, 2009] and so on, to fit the entire solutions for every value of the parameter. It should be noted that there are also several works involving the bi-parametric solution path. Wang et al. [Wang et al., 2008] works with respect to only one parameter of ϵ-SVR while the other parameter is fixed. Bach et al. [Bach et al., 2006] search the space (C+, C-) of 2C-SVM by using a large number of parallel oneparametric solution paths. Rosset's model [Rosset, 2009] follows a large number of one-parametric solution paths simultaneously. Essentially, they all follow one-parametric solution paths in bi-parameter space in different ways, and none of them can explore all solutions for every parameter pair. To address the second difficulty, a global search method [Yang and Ong, 2011] was recently proposed based on the solution path. It can find the global optimal parameters in K-fold CV for C-SVM [Yang and Ong, 2011] and ν-SVM [Gu et al., 2012]. The power of the method is proved by theoretical and empirical analysis for model selection. Therefore, it is highly desirable to design an extension version for CV on the biparametric problem (e.g. CS-SVM) based on fitting all solutions for each parameter pair.

1 Actually, 2ν-SVM is equivalent to 2C-SVM as proved in [Davenport et al., 2010]. For the sake of convenience, we do not distinguish the names of 2C-SVM and CS-SVM hereafter unless explicitly mentioned.

The contributions of this paper can be summarized as follows. (i) We propose a bi-parameter space partition (BPSP) algorithm, which can fit all solutions for every parameter pair (C+,C-). To the best of our knowledge, it is the first such contribution. (ii) Based on the bi-parameter space partition, we propose a K-fold cross validation algorithm for computing the global optimum parameter pairs of CS-SVM. Experimental results demonstrate that the method has better generalization ability than various kinds of grid search methods, however, with less running time.

2 CS-SVMand KKT conditions

We reformulate the primal formulation of 2C-SVM with λ = 1 C++C- and η = C+ C++C- , as presented in (3), and name it (λ, η)-SVM. Fig. 1 shows the relation between the (C+, C-) and (λ, η) coordinate systems. Specifically, the region of C+ > 0, C- > 0, and C+ + C- ≥ 1 corresponds the region of 0 < λ ≤ 1 and 0 ≤ η ≤ 1. Thus, the whole region of (C+, C-) can be explored in 1.5 square units by searching the region of (0, 1] × [0, 1] in the (λ, η) coordinate system, and the lower triangle region of [0, 1]×[0, 1] in the (C+, C-) coordinate system.

s.t. yi (〈w, φ(xi)〉 + b) ≥ 1- ξi, ξi ≥ 0, i = 1, · · · , l The corresponding dual of (3) is

Letting gi = 1λ (∑j∈S αjQij + yib′′) -1, from the KKT theorem [Boyd and Vandenberghe, 2009], we obtain the following KKT conditions for (4):

where b′ = λb′′ , b′′ is the Lagrangian multiplier corresponding to the equality constraint in (4). According to the value of gi, a training sample set S is partitioned as π(λ, η) = (M(λ, η), E(λ, η),R(λ, η)), where M(λ, η) = {i : gi = 0, 0 ≤ αi ≤ 1-yi+2yiη 2 }, E(λ, η) = {i : gi < 0, αi = 1-yi+2yiη 2 }; R(λ, η) = {i : gi > 0, αi = 0}.

Similar to (5)-(6), we can give the KKT conditions for (2). Accordingly, the set S has the partition π(C+,C-) = (M(C+, C-), E (C+, C-), R(C+, C-)).

3 BPSP using Critical Regions

3.1 Detecting the Critical Convex Polygon Region

Given a partition π(λ0, η0), we have the critical region CR(λ0, η0) = {(λ, η) ∈ (0, 1]×[0, 1] : π(λ, η) = π(λ0, η0)} induced by the bi-parametric piecewise linear solution. Theorem 1 shows that CR(λ0, η0) is a convex polygon region.

Theorem 1. The set CR(λ0, η0) is a convex set and its closure is a convex polygon region.

When adjusting λ and η, the weights of the samples inM and the variable b′ should also be adjusted accordingly. From gi = 0, ∀i ∈ M, and the equality constraint (6), and let ˜ gi = λ(gi + 1), we have the following linear system:

If 1M defined as the |M|-dimensional column vector with all ones, and let yM = [y1, · · · , y|M|]T , the linear system

(7)-(8) can be rewritten as:

Let R = ˜Q-1 , the linear relationship between [∆b′ ∆αTM]T and [∆λ ∆η]T can be obtained as follows:

Substituting (10) into (7), we can get the linear relationship between ∆˜ gi (∀i ∈ S) and [∆λ ∆η]T as follows:

When adjusting λ and η, meanwhile keeping all the samples satisfying the KKT conditions, the following constraints should be kept:

Obviously, the set of above inequalities is a convex polygon region. The compact representation can be obtained after removing redundant inequalities from (12)-(14) by the vertex enumeration algorithm [Avis and Fukuda, 1992]. Here, we use CR(λ0, η0) to denote the compact representation of (12)- (14).

Similarly, given a partition π(C0+, C0-), we can obtain the critical region CR(C0+, C0-) as (12)-(14). Obviously, it is also a convex polygon region.

3.2 Critical-Regions-Based BPSP Algorithm

This section tries to find all critical convex polygon regions in (0, 1] × [0, 1] for the (λ, η) parameter space, and the lower triangle region of [0, 1] × [0, 1] for the (C+,C-) parameter space. It means all solutions of CS-SVMwould be obtained.

An intuitive idea to find all the regions is using a progressive construction method. Before designing this progressive construction algorithm, there are three problems which should be answered. (i) How do we give an initial solution of the first critical convex polygon region for (λ, η)-SVM and 2C-SVM? (ii) How do we handle the problem of overlapped critical convex polygon regions? (iii) How do we find the next critical convex polygon regions based on the current one? Our answers to the three problems are as follows, which derive a recursive bi-parameter space partition algorithm (i.e., Algorithm 1).

Initialization A simple strategy for initialization is directly using the SMO technology [Cai and Cherkassky, 2012] or other quadratic programming solvers to find solution for a parameter pair in (0, 1]×[0, 1] for the (λ, η) parameter space, or the lower triangle region of [0, 1] × [0, 1] for the (C+, C-) parameter space, respectively. Here, a method without requiring any numerical solver is presented in Lemma 1 and 2, which can directly give the solutions for the parameter pairs of (λ, η)-SVM and 2C-SVM, respectively, under some conditions.

Lemma 1. All the αi = 1-yi+2yiη 2 , which is the optimal solution of the minimization problem (4) with η = |S-| |S| and λ ≥ 12 (maxi∈S+ ∑j∈S αjQij + maxi∈S- ∑j∈S αjQij ).

Lemma 2. All the αi = C++C-+yi(C+-C-) 2 , which is the optimal solution of the minimization problem (2) with C+ C- =

Partitioning the Parameter Space The minimization problem (4) or (2) can not be guaranteed to be strict convex in many real-world problems. There exists the phenomenon of overlapped critical convex polygon regions (see Fig. 2 (a)). This makes it difficult to find all critical convex polygon regions by a progressive construction method. A parameter space partition method is presented by Theorem 2 [Borrelli, 2003], where A and b are issued from the compact representation of inequalities (12)-(14), which can be computed by the vertex enumeration algorithm [Avis and Fukuda, 1992]. m is the number of inequalities in the compact representation. (ρ, ϱ) is the shorthand implying (λ, η) and (C+,C-) hereafter.

Theorem 2. Let X ⊆ R2 be a convex polygon region, and R0 = {(ρ, ϱ) ∈ X : A[ρ ϱ]T ≤ b} be a convex polygon subregion of X , where A ∈ Rm×2 , b ∈ Rm×1 , R0 ̸= ∅. Also let

then {R0,R1, · · · ,Rm} is a partition of X, i.e., ⋃mi=0 Ri = X, and Ri ∩Rj = ∅, ∀i ̸= j, i, j ∈ {0, 1, · · · ,m}.

Theorem 2 defines a partition procedure which consists of considering one by one the inequalities of R0. See Fig. 2 (b), the four inequalities of R0 induce four disjoint subregions of X (i.e., R1, R2, R3, and R4), respectively, and ⋃4i=0 Ri = X . Obviously, this partition method can be used to handle the problem of overlapped critical convex polygon regions.

Computing Solution for a Parameter Pair in Ri For each convex subregion Ri, similar to the above Initialization, we need to find the solution for a parameter pair (ρi, ϱi) in Ri, and compute the corresponding CR(ρi, ϱi), then partition Ri based on CR(ρi, ϱi) ∩ Ri. Repeat the above steps until the full parameter space are partitioned with critical convex polygon regions (see Fig. 2c and 2d). Obviously, how to find the solution for a parameter pair in Ri is the key to compute the next critical convex polygon region. A simple strategy is using the SMO technology [Cai and Cherkassky, 2012] or other quadratic programming solvers similar to the Initialization. Instead, Theorem 3 allows us to compute α and ˜ g (g) for a parameter pair in the subregion of Ri adjacent to R0 directly.

Theorem 3. Supposing X ⊆ R2 is a convex polygon region, CR(λ0, η0)∩X def = R0 or CR(C 0+,C0-)∩X def = R0, R0 has the partition π, and {R0, R1, · · · , Rm} is a partition of X as Theorem 2. ∀i ∈ {1, · · · ,m}, if Ri ̸= ∅, the i-th inequality Ai [ρ ϱ]T ≤ bi of CR only corresponds to the t-th sample of S,

1. from the left part of (12), there will exist a subregion of Ri adjacent to R0 with the partition π = (M \ {t}, E , R ∪ {t});
2. from the right part of (12), there will exist a subregion of Ri adjacent to R0 with the partition π = (M \ {t}, E ∪ {t}, R);
3. from (13), there will exist a subregion of Ri adjacent to R0 with the partition π = (M∪ {t}, E \ {t}, R);
4. and from (14), there will exist a subregion of Ri adjacent to R0 with the partition π = (M ∪ {t}, E , R \ {t}).

The partition π for the subregion of Ri adjacent to R0 is given by Theorem 3. Thus, we can update the inverse matrix R corresponding to the extended kernel matrix ˜Q in time O(|M|2 ) as the method described in [Laskov et al., 2006], and compute the linear relationships between ∆b′ (∆b′′ ), ∆αTM, ∆˜ g (∆g) and [∆ρ ∆ϱ] as (10)-(11). Further, α(ρi, ϱi) and ˜ g(λi, ηi) (g(C+, C-)), where (ρi, ϱi) is a parameter pair in the subregion of Ri adjacent to R0 with the partition π, can be computed by (12)-(14) directly.

Algorithm 1 CR-BPSP (CRs-based BPSP algorithm)

Input: α(ρ0, ϱ0), ˜ g(λ0, η0) (g(C0+,C0-)), π(ρ0, ϱ0), a convex polygon region X with (ρ0, ϱ0) ∈ X .

• 1: Detect CR(ρ0, ϱ0) according to (12)-(14).

Output: P (a partition of X in a nested set structure).

• 2: Let R0 := CR(ρ0, ϱ0) ∩X, and P := {R0}.
• 3: Partition the parameter spaceX with {R0,R1, · · · ,Rm} (c.f. Theorem 2).
• 4: while i ≤m&Ri ̸= ∅ do
• 6: Compute α and ˜ g (g) for a parameter pair (ρi, ϱi) in the subregion of Ri adjacent to R0.
• 5: Update π, R, for the subregion of Ri adjacent to R0.
• 7: Pi:=CR-BPSP(α(ρi, ϱi), ˜ g(λ0, η0) (g(C 0+,C0-)), π(ρi, ϱi), Ri). {Pi is the partition of Ri.}
• 9: i := i +1.
• 8: Update P := P ∪ {Pi}.
• 10: end while

4 BPSP using Invariant Regions

4.1 Invariant-Regions-Based BPSP Algorithm

Based on the above bi-parameter space partition, the decision function of CS-SVM can be obtained as f(λ, η)(x) = 1λ (∑j∈S αj (λ, η)yjK(xj , x) + b′ (λ, η)) for all (λ, η) in (0, 1] × [0, 1], and f(C+,C-)(x) = ∑j∈S αj(C+,C-)yjK(xj, x) + b′′ (C+,C-) for all (C+,C-) in the lower triangle region of [0, 1] × [0, 1]. Given a validation set V = {(˜ x1, ˜ y1), · · · , (˜ xn, ˜ yn)}, and assuming C(-,+) and C(+,-) are the misclassification costs of false negative and false positive respectively (no costs for the true positive and the true negative), the empirical error on the validation set can be computed as C(ρ, ϱ) = 1n ∑ni=1 C (sign (f (ρ, ϱ)(˜ xi)) , ˜ yi). To select the parameter pairs with the lowest empirical error, we need to investigate empirical errors for all parameter pairs.

Detecting the Invariant Convex Polygon Region According to the sign of f(˜ xi), the validation set V can be partitioned as:

Based on the partition, we have the invariant region IR(ρ0, ϱ0) = {(ρ, ϱ) ∈ CR(ρ0, ϱ0) : ˜ π(ρ, ϱ) = ˜ π(ρ0, ϱ0)}, in which the empirical error obviously remains unchanged. Theorem 4 shows that IR(ρ0, ϱ0) is also a convex polygon region. Thus, we can compute all empirical errors though finding invariant convex polygon regions in the two parameter spaces.

Theorem 4. The sets IR(ρ0, ϱ0) is a convex set and its closure is a convex polygon region.

∀(λ, η) ∈ CR(λ0, η0), according to (10), we can get the linear relationship between ∆f(˜ xi) and [∆λ ∆η] as follows:

Combining (16) with the constraint of ˜ π(λ0, η0), we can get the following constraints:

∀i ∈ I+(λ0, η0) :

Obviously, the closure of inequalities (17)-(18) is a convex polygon region, and the compact representation is denoted IR(λ0, η0). The same analysis can be extended to IR(C0+, C0-).

Partitioning Each CR with IRs To find all invariant convex polygon regions in the whole parameter space, we use the strategy of divide and conquer (i.e., find all invariant convex polygon regions for each critical convex polygon region).

Thus, similar to Algorithm 1, a recursive procedure (called IR-BPSP) can be designed to find all invariant convex polygon regions and compute the corresponding empirical errors for each critical convex polygon region. A nested set structure for the output of IR-BPSP can be retained based on Theorem 2. The nested set structure will speed up finding the global optimal solution for K-fold CV in Section 4.2. Combining all results of the critical convex polygon regions based on the framework of Algorithm 1, we can obtain the empirical errors for all parameter pairs of (λ, η) in the region of (0, 1] × [0, 1] as shown in Fig. 3a, and the empirical errors for all parameter pairs of (C+, C-) in the lower triangle region of [0, 1]×[0, 1] as shown in Fig. 3c.

4.2 Computing the Superposition of K BPSPs

The validation set V is randomly partitioned into K equal size subsets. For each k = 1, · · · ,K, we fit the CS-SVM model with a parameter pair (ρ, ϱ) to the other K - 1 parts, which produces the decision function f(ρ, ϱ)(x) and compute its empirical error in predicting the k part C k (ρ, ϱ) = 1 |Vk | ∑i∈Vk C (sign (f k (ρ, ϱ)(˜ xi)) , ˜ yi ). This gives the CV error CVC(ρ, ϱ) = 1K ∑Kk=1 C k (ρ, ϱ). The superposition of K invariant region partitions can be easily computed for selecting the best parameter pairs of (C+, C-) in R+ ×R+ (see Fig. 3b and 3d), which is omitted here.

5 Experiments

Design of Experiments We compare the generalization ability and runtime of BPSP with other three typical model selection methods of CS-SVM: (1) grid search (GS): a two-step grid search strategy is used for 2C-SVM. The initial search is done on a 20 × 20 coarse grid linearly spaced in the region {(log2 C+, log2 C-)| - 9 ≤ log2 C+ ≤ 10,-9 ≤ log2 C- ≤ 10}, followed by a fine search on a 20 × 20 uniform grid linearly spaced by 0.1; (2) a hybrid method of one-parametric solution path searching on η and grid searching on λ (SPη+GSλ): λ is selected by a two-step grid search in the region {log2 λ| - 9 ≤ log2 λ ≤ 10} with the granularity 1 and followed by 0.1; (3) a hybrid method of oneparametric solution path searching on λ and grid searching on η (SPλ+GSη): η is selected by a two-step grid search in the region {η|0 ≤ η ≤ 1} with the granularity 0.1 and followed by 0.01.

Implementation We implemented SPη+GSλ, SPλ+GSη, and our BPSP in MATLAB. To compare the run-time in the same platform, we did not directly modify the LIBSVM software package [Chang and Lin, 2011] as stated in [Davenport et al., 2010], but implemented the SMO-type algorithm of 2C-SVM in MATLAB. All experiments were performed on a 2.5-GHz Intel Core i5 machine with 8GB RAM and MATLAB 7.10 platform. C(-, +) and C(+, -) are the misclassification costs of false negative and false positive respectively. To investigate how the performance of an approach changes with different settings in misclassification cost, C(-, +) was set to 2, 5, 10 for normal datasets of binary classification, and the class imbalance ratio for imbalanced datasets, respectively, while C(+,-) was fixed at 1. Gaussian kernel K(x1, x2) = exp(-κ‖x1 - x2‖2 ) was used in all the experiments with κ ∈ {10-3 , 10-2 , 10-1 , 1, 10, 102 , 103}, where the value of κ having the lowest CV error was adopted. For BPSP (or SPη+GSλ, SPλ+GSη), our implementation returns a center point from the region (or the line segment) with the minimum error.

Datasets The sonar (Son), ionosphere (Ion), diabetes (Dia), breast cancer (BC), heart (Hea), and hill-valley (HV) datasets were obtained from the UCI benchmark repository [Bache and Lichman, 2013]. The spine image (SI) dataset collected by us is to diagnose degenerative disc disease depending on five image texture features quantified from magnetic resonance imaging. The dataset contains 350 records, where 157 are normal and 193 are abnormal. They are normal datasets for binary classification. Ecoli1, Ecoli3, Vowel0, and Vehicle0 are the imbalanced datasets from the KEEL-dataset repository 2 . Their class imbalance ratios are varying from 3.25 to 9.98.

We selected 30% from a dataset once as a validation set. The validation set was used with a 5-fold CV procedure to determine the optimal parameters. We then randomly partitioned each dataset into 75% for training and 25% for testing for many times. Each time, we removed the instances appearing in the validation set from the testing set to guarantee that the test set of each run is disjoint from the validation set.

Experimental Results The CV errors are presented in Table 1 for 5-fold CV of each method. It is easily observed that BPSP obtains lowest CV error for all datasets and settings of C(-,+). This is reasonable because GS and SPη+GSλ, SPλ+GSη are points-based and lines-based grid search method, respectively, however, BPSP is a regionsbased method which covers all candidate values in the biparameter space, and give the best choices from them. Noted that SPη+GSλ, SPλ+GSη and BPSP can have the same CV error for some datasets, because both of them find the optimal on these datasets. BPSP always find an optimal parameter pair, and SPη+GSλ, SPλ+GSη can also find an optimal sometimes. Based on the optimal parameters in Table 1, the empirical errors on each dataset in different methods over 50 trials are presented in Figure 4 as C(-,+) = 2, 5, 10, and ratio, respectively. The results show that BPSP has better generalization ability than GS, SPη+GSλ, and SPλ+GSη generally. Especially, BPSP has the best stability, because it returns a center point from the optimal region with the minimum error.

2 http://sci2s.ugr.es/keel/imbalanced.php.

The empirical running time (in minutes) for different algorithms on each dataset is also presented in Table 1, which is the average result on the seven different values of κ. It is easy to find that GS method has the longest running time. Because SPη+GSλ and SPλ+GSη searche a large number of parallel one-parametric solution paths, we find that BPSP has the less running time than SPη+GSλ and SPλ+GSη.

6 Conclusion

We proposed a bi-parameter space partition algorithm for CS-SVM which can fit all solutions for each parameter pair (C+,C-). Based on the space partition, a K-fold cross validation algorithm was proposed which can find the global optimum parameter pair. The experiments indicate that our method has better generalization ability than various kinds of grid search methods, however, with less running time. In future work, we plan to extend this framework to a more general formulation which can cover bi-parametric learning models, and even multi-parametric learning models [Mukhopadhyay et al., 2014].

Acknowledgments

This work was supported by the NSF of China (61232016 and 61202137) and the U.S. NSF (IIS-1115417).