Accelerated Asynchronous Greedy Coordinate Descent Algorithm for SVMs

1 Introduction

Support vector machines (SVMs) [Vapnik and Vapnik, 1998] is an elegant and widely applied method with in-depth theoretical analysis in machine learning. However, the traditional SVM solvers (e.g. LIBSVM [Chang and Lin, 2011]) are not scalable in the current big data era because of the ordinary handle of large kernel matrix. To speed up the computation of SVMs, the techniques of kernel approximation and parallel computations have been exploited.

Kernel approximation approaches are trying to approximate the kernel matrix, then solving a linear SVM [Williams and Seeger, 2001; Rahimi and Recht, 2008; Gu et al., 2018b]. However, the kernel approximation approaches cannot obtain the exact solution of SVMs, since the approximate kernel is used. In this paper, we focus to obtain the exact solution of SVMs.

Most parallel algorithms for training SVMs are based on the synchronous model [Chang and Lin, 2011; Zhao and Magoules, 2011; You et al., 2015]. However, without waiting sync point the asynchronous parallel algorithm [Gu et al., 2018a] is much more efficient. To the best of our knowledge, the only asynchronous parallel kernel learning method is the asynchronous parallel greedy coordinate descent (AsyGCD) algorithm for C-SVC [You et al., 2016]. Specifically, at each iteration, workers of AsyGCD asynchronously conduct greedy coordinate descent updates on a block of variables. Comparing with the asynchronously stochastic coordinate descent algorithms [Liu et al., 2015; Duchi et al., 2015], the AsyGCD algorithm can achieve much faster convergence speed due to the greedy selection of updated coordinates.

As mentioned above, AsyGCD is a state of the art solver for solving C-SVC [You et al., 2016]. Thus, it is natural to have the following two questions: 1) Can AsyGCD be further speeded up? 2) Can AsyGCD solve other SVMs, besides C-SVC? To address these issues, in this paper, we propose an asynchronous accelerated greedy coordinate descent algorithm (AsyAGCD) to further accelerate AsyGCD by the technique of active set and extend it to handle generalized SVMs. Active set (or called as shrinking) is a technique to improve the efficiency of dual coordinate descent methods [Hsieh et al., 2008] which has been successfully used in many works [Chang and Lin, 2011; Chiang et al., 2016; De Santis et al., 2016; Birgin and Mart´ınez, 2002] including LIBSVM. Experiment results on a variety of datasets and learning applications confirm that our AsyAGCD is much faster than the existing SVM solvers (including AsyGCD).

2 Generalized SVMs and AsyGCDAlgorithm

In this section, we first give a generalized SVM formulation, and then give a brief review of the AsyGCD algorithm.

2.1 Generalized SVMs

In this paper, we consider a generalized SVM formulation as follows.

where Q is an ℓ × ℓ symmetric and positive semi-definite matrix, p and α are the vectors with the size of ℓ. In the following, we show that C-SVC and ϵ-SVR can be formulated as the form of (1).

Given a training set S = {(xi, yi)}li=1 with xi ∈ Rd , and yi ∈ {+1,-1} for binary classification or yi ∈ R for regression. To make the dual formulations of C-SVC and ϵ-SVR free of equality constraints, we append x with an additional dimension ˆxTi ← [xTi , 1] as mentioned in [Hsieh et al., 2008]. Thus, we have a new training sample set ˆ S = {(ˆxi, yi)}li=1 . For C-SVC, we let Qij = yiyjK(ˆxi, ˆxj), K(ˆxi, ˆxj) = 〈φ(ˆxi), φ(ˆxj)〉 and p = -1. Thus, we can verify that C-SVC is a special case of (1).

For ϵ-SVR, we introduce another vector z with zi ∈ R to denote the targets associated to xi. To present the ϵ-SVR in a compact form, we double the set {(ˆxi, zi)}li=1 , and define an extended training set as ˆ S = ˆ S1 ∪ ˆ S2 , where ˆ S1 = {(ˆxi, zi, yi = 1)}li=1 and ˆ S2 = {(ˆxi, zi, yi = -1)}li=1 . Thus, the dual ϵ-SVR can be presented as follows.

where ϵ defines the ϵ insensitive zone loss in the primal formulation of ϵ-SVR. Thus, it is easy to verify that the dual ϵ-SVR (2) is a special case of (1), if setting p = [ϵ1T - z; ϵ1T + z].

2.2 AsyGCDAlgorithm

AsyGCD Algorithm was originally proposed by [You et al., 2016] to solve the optimization problem (1) with p = -1. In this section, we modify the AsyGCD algorithm such that it can solve the generalized SVM formulation (1). Firstly, we present the i-th coordinated subproblem of (1) as follows.

where ei = [0, · · · , 0 ︸ ︷︷ ︸ i-1 , 1, 0, · · · , 0]T is the i-th indicator vec-

tor. For each iteration of AsyGCD, we update the solution α with αi = αi + δ∗ , where δ∗ is obtained by solving the corresponding i-th coordinated subproblem (3). Fortunately, δ∗ has the closed form solution as follows.

where PΩ(α) is the Euclidean projection of α onto Ω, as mentioned in [You et al., 2016] the step size is set as 1/Qii.

We summarize our AsyGCD algorithm in Algorithm 1. For each iteration of AsyGCD, the i-th coordinate is selected greedily according to |∇+i F(α)|, where ∇+i F(α) is defined as ∇+i F(α) = αi - P[0,C](αi -∇iF(α)). To make AsyGCD run efficiently, we maintain the gradient ∇F (α) as G = ∇F (α) = Qα + p.

Note that, the asynchronously parallel algorithm would overwrite the variables in the shared memory due to without any lock. To ease the phenomenon of overwriting in AsyGCD, [You et al., 2016] partitioned the dual variables into p groups (i.e., S1 ∪ S2 ∪ · · · ∪ Sp = {1, · · · , l}, and Si ∩ Sj = ∅, ∀i, j, i ̸= j), where p is the number of cores in a multi-core shared memory machine.

• Input : The training set ˆ S, an initial α = 0 and G = -1. Output: The vector α. 1 Each thread repeatedly performs the following updates in parallel: 2 for t = 1, 2, · · · do 3 Pick i = argmaxi |∇+i F(α)|; 4 Compute δ∗ i according to (4); 5 for j = 1, 2, · · · , l do 6 Gj = Gj + δ∗ i Qj,i using atomic update; 7 end 8 αi = αi + δ∗ i . 9 end 10 return α.

Algorithm 1: AsyGCD Algorithm

3 Accelerated AsyGCDAlgorithm

In this section, we first introduce the active set technique. Based on the active set technique, we propose our AsyAGCD. Finally, we compare the time complexity between AsyAGCD and AsyGCD.

3.1 Active Set Technique

To present the active set technique smoothly, we first give the Karush-Kuhn-Tucker (KKT) conditions [Boyd and Vandenberghe, 2004] to the generalized SVM formulation (1), then give our active set technique to (1).

KKTConditions: According to theKKT theory [Bertsekas, 1999], the KKT conditions of (1) are presented as follows.

where ∇F (α) = Qα+p is the gradient to the function F (α) with respect to α.

According to the gradient ∇F(α) used in (5), we can define the projected gradient ∇P f(α) as follows.

If α satisfies the KKT conditions, we have that

Thus, Eq. (7) is another form of the KKT conditions (5).

If α does not satisfy the KKT conditions, it is easy to conclude that maxi ∇Pi F(α) > 0 or mini ∇Pi F(α) < 0. These two values (i.e., maxi ∇Pi F(α) and mini ∇Pi F(α), later we will define them as M and m respectively in (9)) measure how the solution α violates the KKT conditions.

Active Set Technique: As mentioned before, AsyGCD algorithm is directly solving the optimization problem of CSVC. In this section, we will use the active set technique to accelerate the AsyGCD algorithm by solving a smaller optimization problem.

Given an active set A (a subset of {1, · · · , ℓ}), correspondingly, we can define an inactive set ¯ A as ¯ A = {1, · · · , l}-A. If we have the prior knowledge that the variables α ¯ A are fixed, and only the variables αA is active, the original optimization (1) can be reduced as a smaller optimization problem (8) as follows.

where QAA and QA ¯ A are the sub-matrices of Q. Thus, we can save the computational time because of solving a smaller optimization problem (8). Now the only question is that how to set the active set A and the inactive set ¯ A. To answer this question, we first give Theorem 1 which can be proved similarly according to the proof in [Hsieh et al., 2008].

Theorem 1. Assume {αk }∞k=1 is an infinite sequence of (1) generated by Algorithm 1. Specifically, αk+1 = αk + δ∗eik , where the index ik is selected by Algorithm 1 for the k-th iteration, and δ∗ is the optimal solution of

Assume every limit point of {αk }∞k=1 is a stationary point, and let α∗ be the convergent point of {αk }∞k=1 . We have that

1. If α∗i = 0 and ∇iF(α∗ ) > 0, then ∃ki such that ∀k ≥ ki,αki = 0 .
2. If α∗i = C and ∇iF(α∗ ) < 0, then ∃ki such that ∀k ≥ ki,αki = C .
3. And limk→∞∇P F(αk ) = 0.

Theorem 1 shows that the variables αi at upper and lower bounds would be fixed from a certain iteration. Thus, we can set the active set A and the inactive set ¯ A according to Theorem 1. Specifically, we define M and m as follows.

If one of the following two conditions is held, we can define the inactive set ¯ A as the collection of the elements satisfying the following conditions.

where

Notice that ¯M and ¯m are the auxiliary variables toM andm, such that ¯M > 0 and ¯m< 0. Correspondingly, the active set A is defined as A = {1, · · · , l}- ¯ A.

During the optimization process, the gradients for i ∈ A can be obtained directly during solving problem (8) due to the following relation:

However, the gradients for i ∈ ¯ A need to be recalculated and this step may be very time-consuming if we shrink many elements. In order to reduce the computational cost of this reconstruction, we maintain a vector ¯G ∈ Rℓ throughout iterations [Chang and Lin, 2011].

We use the fact that αj = 0 or αj = C if j ∈ ¯ A. Thus, for i ∈ ¯ A, we have

Eq. (13) can be used to reduce the computational cost of Gi because normally there exists a large proportion of samples are bounded (i.e., αj = 0 or αj = C). Thus, we can just focus on the computation on the elements in the active set A.

3.2 AsyAGCD

In this section, we introduce the accelerated asynchronous greedy coordinate descent algorithm (AsyAGCD) based on the active set technique. The AsyAGCD algorithm is presented in Algorithm 2. Algorithm 2 includes two parts, i.e., the inner loop and the outer loop.

Inner loop: The inner loop (see lines 5-21 of Algorithm 2) is mainly to optimize the subproblem (8) with the active set technique. In order to simplify the complexity of reconstructing gradients (see line 22 of Algorithm 2), we maintain the vector ¯G (see line 13 of Algorithm 2) for each update. Simultaneously, while arriving the shrinking condition (i.e., k%MaxShrinkIter = 0, where MaxShrinkIter is predefined by user as the maximum number of iterations for doing the shrinking), the inner loop (see lines 16-18 of Algorithm 2) do the shrinking (i.e., Algorithm 3), to reduce the size of active set A. Algorithm 3 describes the procedures of shrinking. All these computations can be run in parallel. Line 5 of Algorithm 2 is used to check whether the solution to the subproblem is an approximate optimal solution.

In addition, it should be noted that the training samples are partitioned into p independent group (see line 2 of Algorithm 2) similar to AsyGCD, where p is the number of cores in a multi-core shared memory machine, and each group is associated to a core. Thus, we update the active set A for each group to ensure that, the active set A is partitioned into the corresponding groups (see line 18 of Algorithm 2).

Outer loop: The outer loop is mainly used for checking the termination condition and reconstructing the active set. Specifically, we use the global variables Mˆ S andmˆ S to check the quality of the solution (see lines 3 and 22-24 of Algorithm 2). Given a tolerance ε, if the condition M -m ≤ ε is satisfied, the ε-approximate solution of (1) is obtained, and the AsyAGCD algorithm terminate. If the condition M -m ≤ ε is not satisfied, we need to reconstruct the active set and then run the inner loop.

3.3 Comparison with AsyGCD

In this section, we compare our AsyAGCD with AsyGCD with respect to the computational complexity. Table 1 provides the computational complexities of each iteration for AsyGCD and AsyAGCD, and the total computational complexities of AsyGCD and AsyAGCD (see the last row of Table 1).

For the AsyGCD algorithm (i.e., Algorithm 2), You et al. [You et al., 2016] pointed out the computational complexity of greedy coordinate selection is O(ℓ). The computational complexity of updating the gradients is also O(ℓ). Assume the size of total iterations of AsyGCD (or AsyAGCD) is IterTotal, the total computational complexity of AsyGCD is IterTotal ×O(ℓ), please see the last row of Table 1.

For the AsyAGCD algorithm (i.e., Algorithm 2), we first consider the computational complexity for the inner loop. Specifically, the computational complexities of greedy coordinate selection, updating GA, calculating MA and mA and shrinking are O(|A|). The computational complexity of updating the values of ¯G is O(ℓ). For the outer loop of Algorithm 2, The computational complexity of reconstructing gradients is O(ℓ × (ℓ - |A|)). The computational complexity of calculating Mˆ S and mˆ S is O(ℓ). Thus, the total computational complexity of AsyAGCD can be described as (I terT otal+I terInner)×O(|A|)+#GbarCount×O(ℓ)+ I terOuter ×O(ℓ × (ℓ - |A|)), where #GbarCount denotes the times of updating ¯G. Note that the above |A| represents the average size of the active set A. To have a deep insight to the computational complexity of AsyAGCD, we give the following two remarks.

Remark 1. Generally, we have the relationship (IterT otal + IterInner) ≫ #GbarCount > IterOuter. Thus, the computational complexity of our AsyAGCD is dominated by the average size of the active sets and the iteration number IterT otal + IterInner.

Remark 2. Generally, the active set size may become smaller and smaller, as the number of iterations increases. Then the average active set size |A| of Table 1 is smaller accordingly and shrinking technique can shorten the training time.

According to the above Remarks 1 and 2, if the number of iterations is large, our AsyAGCD would have a speedup compared to AsyGCD, because of using the active set technique. This is because the subproblem (8) would become very lightweight (i.e., |A| ≪ ℓ in Table 1) and the cache hit ratio improves simultaneously (the cache strategy is widely used in [Chang and Lin, 2011; You et al., 2015]). In addition, we can update ¯G quickly because the size of A would become small enough. Thus, our AsyAGCD can converge much faster than AsyGCD for the second half of iterations.

On the contrary, if the number of iterations is too small, the active set size would be relatively large during the whole procedure. Thus, the active set technique may not be so useful. This is because that, the time saved by the active set technique cannot make up the time introduced by the active set technique (i.e., IterInner ×O(|A|) +#GbarCount × O(ℓ) + IterOuter × O(ℓ × (ℓ - |A|)) as analyzed in Table 1). In this case, the AsyGCD algorithmwithout the active set technique may be faster than our AsyAGCD.

4 Experiments

In this section, we first give the experimental setup, then show our experimental results and discussions.

4.1 Experimental Setup

Design of experiments: In the experiments, we mainly compare the efficiency of our AsyAGCD with the existing state-of-art solver of SVMs. The compared algorithms used in our experiments include:

1. AsyAGCD 1 : Our AsyAGCD algorithm which can solve C-SVC and ϵ-SVR. As mentioned before, AsyAGCD is defaultly with shrinking (i.e., the active set technique).
2. AsyGCD 2 : The AsyGCD algorithm which was proposed by You et al. [You et al., 2016]. Note that AsyGCD can only solve C-SVC.
3. AsyAGCD without shrinking: A new version of AsyGCD without shrinking which can solve ϵ-SVR.
4. LIBSVM(OMP) with shrinking 3 : The parallel version of LIBSVM with shrinking, where a column of kernel matrix is computed parallelly.

For C-SVC, the AsyAGCD, AsyGCD and LIBSVM(OMP) with shrinking methods are compared. Specifically, we compare the accuracies and objective values over the running time among these methods. Note that, because of the different formulation of C-SVC solved by LIBSVM, we only compare the objective values over the running time between AsyAGCD and AsyGCD. For ϵ-SVR, the AsyAGCD, AsyAGCD without shrinking and LIBSVM(OMP) with shrinking methods are compared. We compare the mean absolute errors (MAE) and the mean squared errors (MSE) over the running time among these methods.

1 The code is available in https://sites.google.com/site/jsgubin/ our-software-codes/AsyAGCD code.zip.

2 The code is available in https://github.com/cjhsieh/asyn kernel svm

3 The code is available in http://www.csie.ntu.edu.tw/∼cjlin/ libsvm/

Implementation: We implement our proposed AsyAGCD algorithm by C++ based on AsyGCD, where the shared memory parallel computation is handled via OpenMP [Chandra, 2001].We perform experiments on a 36-core two-socket Intel Xeon E5-2696 machine with 48 GB RAM, where each socket has 18 cores. All the experiments are running with the Gaussian kernel K(x1, x2) = exp(-κ‖x1 - x2‖2 ), the parameter κ, C and ϵ are the default values at LIBSVM. MaxShrinkIter in AsyAGCD is set as 10, 000. All the experiments are running at least 10 times with 20 cores, and using 24GB memory space for kernel caching.

Datasets: Table 2 summarizes the eight datasets used in experiments which divided into two parts according to the tasks, i.e., the binary classification (BC) and regression (R). They are from https://archive.ics.uci.edu/ml/datasets.html and https://www.csie.ntu.edu.tw/∼cjlin/libsvmtools/datasets/. All the inputs of these datasets are mapped into the interval [0, 1].

4.2 Experimental Results and Discussion

C-SVC: Fig. 1 shows the accuracies(objective values) v.s. the running time of different algorithms. From the top graphs, it is easy to find that both AsyAGCD and AsyGCD converge faster than the LIBSVM(OMP) in the most cases. It should be noted that, the advantages of AsyAGCD are not obvious on the Ijcnn1 dataset. This phenomenon shows our AsyAGCD is more applicable to handle large datasets. For the Covtype dataset, the performance of AsyAGCD and AsyGCD is similar. This is because only a small number of variables are shrunk in each shrinking procedure. Thus, the running time saved from the active set technique cannot make up the time consuming introduced by the shrinking technology, as discussed in Section 3.3.

From the bottom graphs in Fig. 1, it is easy to find that our AsyAGCD is faster than AsyGCD in terms of the objective values. This follows the comparison of the computational complexities of AsyGCD and AsyAGCD as discussed in Section 3.3. Note that, in ijcnn1 and skinnoskin dataset, it looks like that the convergence points of AsyGCD and AsyAGCD are different, which contradicts to the theoretical result that, AsyGCD and AsyAGCD should converge to a same value of the objective function. Actually, if running AsyGCD and AsyAGCD for enough time, they can converge to a same objective values. However, the running time shown in figures is limited, which leads that the convergence points of AsyGCD and AsyAGCD look different. This phenomenon also confirms that our AsyAGCD is faster than AsyGCD.

Remark 3. In fact, we observe that most of computational time is wasted on selecting the active set and calculating the kernel vector. During the second half of iterations, selecting the active set is fairly rapid because the size of the active set becomes relatively small. Thus, our AsyAGCD can converge much faster than AsyGCD for the second half of iterations.

ϵ-SVR: Fig.2 shows the MAE values v.s. the running time of different algorithms. It is easy to find that both AsyAGCD

with and without shrinking converge faster than the LIBSVM(OMP). This is because AsyAGCD is an asynchronous parallel algorithm, and LIBSVM(OMP) is a synchronous parallel algorithm. Normally, asynchronous parallel computation is much faster than synchronous parallel computation. In addition, [Steinwart et al., 2011] pointed out that solving an SVM without offset can be significantly faster.

5 Conclusion

In this paper, we propose an asynchronous accelerated greedy coordinate descent algorithm (AsyAGCD) for SVMs. Compared with the existing state-of-art solver AsyGCD, our AsyAGCD has the following two-fold advantages: 1) AsyAGCD is an accelerated version of AsyGCD and converge much faster than AsyGCD for the second half of iterations. 2) Our AsyAGCD can handle more SVM formulations than AsyGCD. We provide the comparison of computational complexity of AsyGCD and our AsyAGCD. Experiment results confirm that our AsyAGCD is much faster than the existing SVM solvers (including AsyGCD). In future, we plan to extend our AsyAGCD to other learning problems [Gu et al., 2015a; 2015b; 2017; Gu and Ling, 2015; Gu et al., 2018c; Huo et al., 2018].

Acknowledgments

Dr. Heng Huang made scientific contributions to this paper. Unfortunately the conference chair doesn't allow us to change the author list. We acknowledge Dr. Heng Huang's contributions. This work was partially supported by the Natural Science Foundation of Jiangsu Province of China (No. BK20161534), Six talent peaks project (No. XYDXX-042) and the 333 Project (No. BRA2017455) in Jiangsu Province and the National Natural Science Foundation of China (No 61573191). BG and HH were partially supported by U.S. NSF-IIS 1302675, NSF-IIS 1344152, NSF-DBI 1356628, NSF-IIS 1619308, NSF-IIS 1633753, NIH R01 AG049371.