Efficient Kernel Learning from Side Information Using ADMM

1 Introduction

Kernel methods have been highly successful in various aspects of machine learning, such as classification, regression, clustering, ranking, and dimensionality reduction [Sch¨olkopf and Smola, 2002]. Because of the central role of the kernel, it is important to identify an appropriate kernel function or matrix for the task at hand. Over the past decade, there have been a large body of literature on this kernel learning problem [Lanckriet et al., 2004; Bach et al., 2004; Cortes et al., 2009]. While a parametric form of the kernel [Chapelle et al., 2003] or a combination of multiple kernels [Sonnenburg et al., 2006] are often assumed, nonparametric kernel learning, which takes no such assumptions, is more flexible and has received significant interest in recent years [Li et al., 2008; Wu et al., 2009; Kulis et al., 2009; Hu et al., 2011; Zhuang et al., 2011; Shang et al., 2012].

To facilitate kernel learning, obviously one has to utilize information from the data. The most straightforward approach is to use class labels. However, obtaining label information may sometimes be expensive and time-consuming. In this paper, we focus on a weaker form of supervisory information, namely, the so-called must-link and cannot-link pairwise constraints [Wagstaff et al., 2001]. These pairwise constraints, or side information, define whether the two patterns involved should belong to the same class or not. Another useful source of information, which is commonly used in semi-supervised learning, is the data manifold [Belkin et al., 2006]. This encourages patterns that are locally nearby on the manifold to have similar predicted labels.

A notable nonparametric kernel learning approach that nicely integrates these two sources of information is the pairwise constraint propagation (PCP) algorithm [Li et al., 2008]. Let n be the number of patterns, andM, C be the sets of mustlink and cannot-link pairs, respectively. PCP is formulated as the following optimization problem

where K = [Kij] ∈ Rn×n is the kernel matrix to be learned (which has to be symmetric and positive semidefinite (psd), denoted K ⪰ 0), L is the graph Laplacian matrix of the data manifold, T = [Tij ] with

T =M∪C∪{(i, i)}ni=1 , and γ is a regularization parameter. In words, the tr(KL) term in (1) encourages K to be aligned with L and thus smooth on the manifold, while (Kij - Tij)2 encourages Kij to be close to 1 for two must-link patterns i, j, and 0 for cannot-link patterns. Note from (2) that all the patterns are also encouraged to be normalized in the kernel representation (i.e., Kii ≃ 1).

As in other nonparametric kernel learning algorithms, the optimization problem in (1) is a semidefinite program (SDP) [Vandenberghe and Boyd, 1996]. Though off-theshelf solvers such as CSDP can be used [Li et al., 2008], they are slow, particularly when n is large or the pairwise constraints are abundant. Recently, Laue [2012] proposed a more efficient hybrid solver that is based on gradually increasing the rank of the K solution. However, each time the rank is increased, it has to perform an approximate eigen-decomposition and nonlinear optimization. As will be demonstrated in Section 4.2, this can be inefficient on large kernel learning problems.

Recently, there have been a lot of interest in scaling up this SDP. In many machine learning problems, a prominent tool that often drastically reduces the computational effort is lowrank approximation [Fine and Scheinberg, 2001]. In the context of kernel learning, this involves replacing the n×n matrix K by V⊤V, where V ∈ Rr×n and r is the rank. The solver in [Burer and Choi, 2006] can then be used, though the line search and computation of the coefficients in each iteration can be expensive. Alternatively, one can optimize V columnby-column using coordinate descent [Hu et al., 2011]. The resultant algorithm is scalable, but the decoupling of variables as required in coordinate descent makes it unable to handle the (soft) constraints of Kii ≃ 1 in (2).

Another popular low-rank approximation option is to replace K by QUQ⊤, as in [Li and Liu, 2009; Wu et al., 2009] and the modified fixed point continuation (MFPC) algorithm [Shang et al., 2012]. Here, Q ∈ Rn×r is a fixed matrix with orthogonal columns (usually the trailing eigenvectors of the graph Laplacian matrix L), and U ∈ Rr×r is the psd matrix to be optimized. Typically, r ≪ n, and thus optimization of U is much more manageable. Moreover, Wu et al. [2009] showed that problem (1) belongs to a special class of SDP called quadratic semidefinite program (QSDP). By exploiting this structure, they reformulated (1) as a semidefinitequadratic-linear program (SQLP), which reduces the complexity from O(n9 ) for the standard QSDP formulation to O(n6.5 ). However, empirically, this still does not scale well. As reported in [Wu et al., 2009], kernel learning on a data set with 2000 instances and 100 pairwise constraints takes about 4000 seconds. Moreover, the fixed Q limits the ability of K to fit the data, leading to possibly deteriorated performance of the learned model.

Instead of solving problem (1) directly, Zhuang et al. [2011] noted that this, together with other similar machine learning problems, can be reformulated as optimizing the inner product of the matrix variable with another fixed matrix. By adding an additional constraint tr(KK) ≤ B, where B > 0 is a user-defined constant, this optimization problem has a closed-form solution of K. The resultant algorithm, called simple nonparametric kernel learning (SimpleNPKL), is significantly more efficient than existing algorithms. However, computing the closed-form K solution requires eigendecomposition at each iteration, which can be expensive on large kernel learning problems with a large number of instances or constraints. Moreover, its performance is sometimes sensitive to the setting of B.

In this paper, our key idea is to decompose the matrix K as V⊤U, where V,U ∈ Rr×n , together with the constraint that V = U. Obviously, this is essentially the same as the commonly-used low-rank approximation V⊤V discussed above [Burer and Choi, 2006; Hu et al., 2011; Li and Liu, 2009], but with the important difference that V and U are now decoupled. Together with the QSDP structure in (1), this allows simple and efficient optimization based on the alternating direction method of multipliers (ADMM) [Boyd et al., 2011], which has been popularly used in diverse fields such as machine learning, data mining and image processing. Experimental results show that the proposed method is as accurate as solving the SDP directly, but is significantly more efficient and scalable.

The rest of this paper is organized as follows. Section 2 first gives a brief introduction to ADMM. Section 3 then describes the proposed solver and its properties. Experimental results on a number of real-world data sets are presented in Section 4, and the last section gives some concluding remarks.

Notations: In the sequel, matrices and vectors are denoted in bold, with upper-case letters for matrices and lower-case for vectors. The transpose of a vector/matrix is denoted by the superscript ⊤, tr(A) denotes the trace of matrix A, A⊙B is the element-wise multiplication of matrices A and B (i.e., [A⊙B]ij = AijBij), and A· B = tr(A⊤B) is their inner product. Moreover, I is the identity matrix, and 0 (resp. 1) is the vector (or matrix) of all zeros (resp. ones).

2 Alternating DirectionMethod of Multipliers (ADMM)

ADMM is a simple but powerful algorithm first introduced in [Glowinski and Marrocco, 1975], and has recently been successfully used in diverse fields such as machine learning, data mining and image processing. The standard ADMM is for solving convex problems. Here, we consider the more general bi-convex problem in [Boyd et al., 2011]:

where F (·, ·) is bi-convex and G(·, ·) is bi-affine 1 . As in the method of multipliers, ADMM considers the augmented Lagrangian of (3): L(x, y,Λ) = F (x, y) + Λ⊤G(x, y) + ρ2 ‖G(x, y)‖2 , where Λ is the vector of Lagrangian multipliers, and ρ > 0 is a penalty parameter. At the kth iteration, the values of x, y and Λ (denoted xk , yk and Λk ) are updated as

Note that while the method of multipliers minimizes L (x, y, Λk ) w.r.t. x and y jointly, ADMM allows easier decomposition of the optimization problem by minimizing them in an alternating manner.

3 Kernel Learning by ADMM

In this section, we present a novel solver for problem (1) based on the ADMM. Section 3.1 shows that the update steps can be efficiently computed. A necessary condition for convergence is provided in Section 3.2, and Section 3.3 discusses the time complexity in each iteration.

As in [Wu et al., 2009; Hu et al., 2011; Li and Liu, 2009], we use a low-rank approximation on the kernel matrix K. In other words, K is approximated as

where V ∈ Rr×n and rank r < n. Problem (1) can then be rewritten as

1 In other words, for any fixed x, y, F(·, y) and F(x, ·) are convex; while G(·, y) and G(x, ·) are affine.

Recall that the graph Laplacian matrix L is symmetric [Luxburg, 2007], this can be further written as

where U ∈ Rr×n , and vi's (resp. ui's) are the columns of V (resp. U). Note that the objective is bi-convex, and the constraint function V - U is bi-affine.

3.1 ADMMUpdate Steps

Let Λ be the Lagrangian multiplier for the constraint V - U = 0. As in Section 2, the augmented Lagrangian of (5) is

and the update steps at the kth iteration of ADMM are

It can be easily shown that

where

Λki is the ith column of Λk , and Ti = {j | (i, j) ∈ T }.

On setting (9), (10) to zero, the solutions of (7) and (8) can be obtained as

for i = 1, · · · , n. Note from (11), (12) that Aki and Bki are always invertible.

To monitor convergence of the algorithm, we require the primal and dual residuals at iteration k + 1

Δprimal = ∥∥Vk+1 -Uk+1∥∥ , Δdual = ρ ∥∥Vk+1 -Vk ∥∥ (14) to be small [Boyd et al., 2011]. Moreover, to improve convergence, it is common to vary the penalty parameter ρ in each ADMM iteration. Specifically, following [Boyd et al., 2011], we update ρ as

The whole procedure is shown in Algorithm 1.

Algorithm 1 KLADMM: Kernel learning by ADMM.

• 1: Input: rank r, V0 = [v0 1 , . . . , v0n ],U0 = [u01 , . . . , u0n ], parameters ρ, ε and IterMax.
• 2: Output: K = Vk⊤Vk (or Uk⊤Uk ).
• 3: k ← 0;
• 4: repeat
• 5: for i = 1 to n do

• 7: end for

• 9: compute the primal and dual residuals in (14);
• 10: update ρ using (15);
• 11: k ← k +1;
• 12: until max(Δprimal,Δdual) < ε or k > IterMax;

3.2 Convergence

With the low-rank decomposition in (4), problem (5) is nonconvex in V and so we can only consider local convergence [Boyd et al., 2011]. As in [Xu et al., 2011], we show below a necessary condition for local convergence.

Lemma 1 Let {(Uk , Vk , Λk )} be a sequence generated by Algorithm 1. If the sequence {Λk } is bounded and satisfies

then Vk - Vk+1 → 0 and Uk - Uk+1 → 0.

Proof 1 First, observe that L(U,V,Λ) is bounded, which follows from the boundedness of Λ and

Moreover, L is strongly convex w.r.t. each coordinate variable of vi and ui. For any vi, it holds that

for any △vi, where L(vi, ·) denotes that all the variables in (7) except vi are fixed. Moreover, v⋆ i is a minimizer of L (vi, ·) if

Note that vk+1 i is a minimizer of L (vi, ·) at the kth iteration. Then, on combining (17) and (18), we have

Similarly, one can show that

Using (17) - (20) and Λk+1 - Λk = ρ(Vk+1 - Uk+1 ), we have

Taking summation of the above inequality from 1 to ∞, and recall that L is bounded, then

Since the second term on the left is bounded, it follows that ∞∑ k=1 (∥∥Vk - Vk+1∥∥2 + ∥∥Uk -Uk+1∥∥2) <∞, fromwhich we immediately have Vk -Vk+1 → 0 and Uk -Uk+1 → 0.

■Proposition 1 Let {(Uk , Vk )} be a sequence generated by Algorithm 1. Then any accumulation point of {(Uk , Vk )} satisfies the Karush-Kuhn-Tucker (KKT) condition for problem (5).

Proof 2 From (13), we have

Lemma 1 implies that the left- and right-hand sides of (21), (22), (23) all go to zero, i.e.,

and thus the KKT conditions for problem (5), namely,

are satisfied.

■

In particular, note from (24) that although the algorithm does not guarantee Uk = Vk at each iteration k, the difference goes to zero as the iteration proceeds. This will also be experimentally demonstrated in Section 4.1.

3.3 Time Complexity

As Aki ,Bki are of size r × r, computing their inverses in (13) takes O(r3 ) time. In practice, we often have r > Ti. For example, please refer to the experimental setup in Section 4. By rewriting Aki as γ ( ρ γ I + PiP⊤i ), where Pi = [ukj1 , . . . , ukj|Ti| ] ∈ Rr×|Ti| , and using the Sherman-MorrisonWoodbury (SMW) identity [Golub and Van Loan, 1996], we have

where ρ γ I + P⊤i Pi is of size |Ti| × |Ti|. Hence, computing (Aki )-1 takes only O(|Ti|3 + r|Ti|) time. The same applies to the computing of (Bki )-1 . The time complexity for each iteration of Algorithm 1 is then

In practice, |Ti| is usually very small, and (26) essentially scales linearly with n.

4 Experiments

As in [Hoi et al., 2007; Li et al., 2008], we study the performance of the proposed approach in the context of data clustering. Specifically, a kernel matrix is learned from the pairwise constraints, which is then used for clustering by the kernel k-means algorithm. Experiments are performed on a number of benchmark data sets (Table 1) that have been commonly used for nonparametric kernel learning [Hoi et al., 2007; Hu et al., 2011; Zhuang et al., 2011].

The constructions of the data manifold and pairwise constraints follow that in [Hoi et al., 2007; Hu et al., 2011;

Zhuang et al., 2011]. Specifically, the similarity between patterns xi, xj is defined as exp(-‖xi -xj‖2 /2σ2 ), where σ is equal to half of the average distance between each pattern and its top-10 nearest neighbor. The 50-nearest-neighbor graph is then used for the manifold of the adult data sets, and the 5nearest-neighbor graph for the others. As for the pairwise constraints, we randomly select 0.6n pairs of patterns that belong to the same class to construct the must-link pairs, another 0.6n pairs of patterns that belong to different classes to construct the cannot-link pairs, and another n must-link (i, i) pairs for each pattern i. On all the data sets, the average number of constraints (i.e., |Ti|) associated with a random pattern i is always around 2.

As for the rank of the kernel matrix, we follow [Zhuang et al., 2011; Hu et al., 2011] and set its value to be the largest r satisfying r(r + 1)/2 ≤ m, where m is the total number of constraints in T . Table 1 shows the rank values (i.e., r) used. As can be seen, they are all much larger than |Ti| ≃ 2, i.e., the average number of constraints associated with a random pattern i. This justifies the use of the SMW identity to compute the matrix inverse in Section 3.3.

The proposed KLADMM algorithm 2 is compared with the following solvers discussed in Section 1. Similar to KLADMM, all are based on a rank-r approximation of K.

1. The hybrid algorithm [Laue, 2012], with the rank of the K solution increasing up to r;
2. SimpleNPKL 3 [Zhuang et al., 2011], which is extended for the quadratic loss in (1);
3. MFPC [Shang et al., 2012], which is based on the lowrank representation QUQ⊤, where the columns of Q are the r trailing eigenvectors of the graph Laplacian L;
4. The SQLP solver [Wu et al., 2009], which is also based on the same low-rank representation as MFPC.

All the algorithms are implemented in MATLAB and run on a PC with i7-3370 3.4GHz CPU and 32G RAM.

The following performance evaluation measures are used:

1. percentage of pattern pairs that are correctly clustered together [Hoi et al., 2007]:

where ci (resp. ˆ ci) is the true (resp. predicted) cluster label of pattern xi, and I (·) is the indicator function that returns 1 when the argument is true, and 0 otherwise;

1. CPU time for learning the kernel matrix K.

4.1 Accuracy and Speed

Results on the clustering accuracy and CPU time are shown in Table 2. SQLP, which shares the same low-rank model as MFPC, has comparable accuracy as MFPC but is much slower (often by one to two orders of magnitude). Hence, its results are not reported here because of the lack of space. Moreover, some solvers become very slow on the larger data sets and so are not run to completion.

2V0 and U0 is initialized randomly, ρ = 100, and IterMax set to 500 (and to 100 on the adult datasets).

3 Recall that an additional B parameter is needed to constrain tr(KK). Preliminary results show that the solution quality can be sensitive to B. Here, we set B = 100n, which performs better in our experiments than the value suggested in [Zhuang et al., 2011].

As can be seen, KLADMM is much faster than SimpleNPKL, which in turn is faster than the hybrid algorithm. In particular, on the largest data set (adult-a4a) that SimpleNPKL can handle, KLADMM is about 30 times faster. MFPC, though is almost as fast as KLADMM, has much inferior accuracy, and cannot be improved even by increasing the size of U. As discussed in Section 1, for efficient optimization, MFPC has to fix the Q matrix in the low-rank representation (QUQ⊤) of K. This severely restricts the ability of K to adapt to the data.

On the other hand, KLADMM, though involves non-convex optimization, is as accurate as the hybrid algorithm which converges to the globally optimal solution of the SDP problem [Laue, 2012]. Figure 1 shows the speed in which the objective value obtained by KLADMM converges towards the globally optimal objective of the hybrid algorithm. As can be seen, KLADMM often converges in fewer than 25 iterations.

Finally, Figure 2 shows the progress of KLADMM's primal and dual residuals with the number of iterations. Recall from (14) that the primal residual Δprimal measures the difference between Vk and Uk . As can be seen, this difference quickly goes to zero, as is discussed in Section 3.2.

4.2 Scalability

In this section, we compare the scalability of the various algorithms, except for MFPC and SQLP which have been shown to have much inferior accuracy (Section 4.1). Experiments are performed on a synthetic 10-dimensional binary classification data set. Patterns from the first class are generated from the normal distribution N(1, I), while those from the second class are from N(-1, I).

In the first experiment, we fix the number of constraints at 1,000. Figure 3(a) shows the CPU time as the number of patterns n increases. As in Section 4.1, KLADMM scales much better than the others. The curves are redrawn as a log-log plot in Figure 3(b). A line corresponding to a linear scaling between time and n is also shown. As can be seen, KLADMM scales almost linearly in n, which agrees with the analysis in Section 3.3. In the second experiment, we fix n at 1,000 and gradually increase the number of constraints (recall that the rank r also increases according to our experimental setup). Results are shown in Figure 4. As can be seen, KLADMM is again much more efficient than the other methods.

5 Conclusion

In this paper, we proposed an efficient solver for nonparametric low-rank kernel learning. Using ADMM, the proposed KLADMM decouples the original optimization problem into computationally inexpensive subproblems. Unlike existing solvers which perform a large eigen-decomposition in each iteration, KLADMM only requires operations on very small matrices. Experimental results on a number of real-world data sets demonstrate that KLADMM is as accurate as directly solving the SDP, but is much faster than existing solvers. Its CPU time scales linearly with the number of patterns, and can be used on much larger kernel learning problems.

In the future, the proposed solver will be extended to optimize other nonparametric kernel learning models.

Acknowledgment

This research was supported in part by the Research Grants Council of the Hong Kong Special Administrative Region (Grant No. 614012), the National Natural Science Foundation of China (No. 61165012) , the Natural Science Foundations Fund (Grant No. 2011FZ074) and the Department of Education Fund (No. 2011Y305) of Yunnan Province.