Learning the Optimal Neighborhood Kernel for Classification

1 Introduction

Kernel methods [Scholk¨opf and Smola, 2002] work by embedding the data into a high-dimensional (possibly infinitedimensional) feature space, where the embedding is implicitly specified by a kernel. They have been applied successfully in many areas such as computer vision and bioinformatics [Scholk¨opf and Smola, 2002; Lanckriet et al., 2004a; 2004b]. The so-called kernel matrix, specified by the inner product of data points in the feature space is symmetric and positive semidefinite. The kernel matrix specifies the geometry of the feature space, and induces a notion of similarity in the sample space. The learning of a good kernel matrix is essential for the successful application of kernel methods. However, it is commonly the case in practice that the "ideal" kernel matrix is unknown and needs to be learnt.

The problem of learning the optimal kernel matrix has been addressed by many researchers. The pioneer work by [Chapelle et al., 2002] learned the parameters of some pre-specified kernel function using gradient descent method. More recent works focus on learning the kernel itself in a nonparametric manner. Lanckriet et al. [2004a; 2004b] proposed the framework of Multiple Kernel Learning (MKL) by representing the learnt kernel matrix as a linear combination of a set of pre-specified kernel matrices, and casted the problem as a semidefinite program or a quadratically-constrained quadratic program. Several efficient methods for accelerating the computation of MKL have been proposed in [Bach et al., 2004; Sonnenburg et al., 2006]. Bousquet and Herrmann [2002] proposed to learn the spectral classes of kernel matrices, and analyzed the complexities of the learnt kernel matrices. Kulis et al. [2006] learned low rank kernel matrices from the given distance and similarity constraints on the data. Luss and Aspremont [2007] proposed to learn a proxy kernel matrix by treating the indefinite kernel matrix as a noisy observation of the positive semidefinite kernel matrix, with the attractive property that the support vectors and the proxy kernel can be simultaneously computed.

In many applications, the kernel matrix generated can be treated as a noisy observation of the ideal one. In this paper, we focus on learning an optimal kernel matrix, called the optimal neighborhood kernel (ONK) matrix, from the pre-specified positive semidefinite kernel matrix for improved classification performance. The ONK matrix is shown to be the summation of the pre-specified kernel matrix and a rankone matrix. We formulate the problem of learning ONK as a convex quartic problem. Our proposed ONK formulation is different from the neighborhood mismatch kernel proposed in [Weston et al., 2003], where the kernel matrix was constructed by utilizing the neighborhood samples.

We propose two iterative methods for solving the resulting convex quartic problem: level method and constrained gradient descent. Level method is a memory-based method from the cutting plane family [Nemirovski, 1994]. Recently, the cutting plane methods have been successfully applied to a number of convex optimization problems in machine learning, e.g., the bundle method employed in [Teo et al., 2007] for solving regularized risk minimization, and the analytic center cutting plane method [Luss and Aspremont, 2007] for solving the robust classification problem with indefinite kernels. Nemirovski [1994] showed that, level method has a better convergence rate than the Kelley bundle method for solving the general nonsmooth convex optimizations, and achieves better practical performance than methods such as subgradient descent. However, since the level method keeps all the "prehistory" information, its computational cost can be high. In contrast, gradient descent is a non-memory-based method. To solve the constrained optimizations, we can extend gradient descent to its constrained case using gradient mapping [Nemirovski, 1994]. We show that one of the key steps in constrained gradient descent is the Euclidean projection onto a particular domain, which is a special case of the singly linearly constrained quadratic programs [Dai and Fletcher, 2006]. We cast the Euclidean projection as a zero finding problem, which can be solved in linear time. Empirical results on several benchmark data sets demonstrate the efficiency and effectiveness of the proposed algorithms.

The rest of the paper is organized as follows: we formulate the learning of the optimal neighborhood kernel in Section 2, present the algorithms for solving the convex quartic problem in Sections 3 & 4, report empirical results in Section 5, and conclude the paper in Section 6.

Notations Vectors are represented in bold notation, e.g., x ∈ Rn, and x ≤ 0 is understood componentwise. Matrices are represented in italic script, e.g., X ∈ Rm×n, and X ⪰ 0 means that X is square, symmetric, and positive semidefinite. ‖X‖F denotes the Frobenius norm of the matrix X, and ‖x‖2 denotes the Euclidean norm of the vector x.

2 Learning the Optimal Neighborhood Kernel

Assume that we are given a collection of training data {(xi, yi)}ni=1 , where xi ∈ Rd , and yi ∈ {±1}. We consider the SVM formulation in the primal form:

where ξ = [ξ1, ξ2, . . . , ξn]T denotes the vector of slack variables, φ(.) is a kernel feature mapping, C is the regularization parameter, and (w, b) determines the hyperplane in the feature space. The dual problem of (1) is given by

where e is a vector containing n 1's, y = [y1, y2, . . . , yn]T , Y = diag(y), K is the kernel matrix specified by the kernel κ(., .) as Kij = κ(xi, xj) = φ(xi)T φ(xj), and α = [α1, α2, . . . , αn]T ∈ Rn can be solved by quadratic programming. For an unknown test sample x, we compute

and x is classified as +1 if g(x) is positive, and -1 otherwise.

2.1 The Optimal Neighborhood Kernel Matrix

We assume that the pre-specified kernel matrix K is the noisy observation of the "ideal" kernel matrix G. We thus propose to learn the kernel matrix G in the neighborhood of K by solving the following optimization problem:

where the parameter ρ > 0 controls the magnitude of the penalty on the square distance between the optimal neighborhood kernel matrix G and the pre-specified kernel matrix K.

It is easy to verify that, the objective function in (3) is convex in G and concave in α. Therefore, this is a convex concave optimization problem, and the existence of a saddle point is given by the well-known von Neumann Lemma [Nemirovski, 1994]. Moreover, since the optimization is strictly feasible, the minimization and maximization can be exchanged, resulting in the following optimization problem:

We first show that the optimal neighborhood kernel matrix G that solves the optimization problem in (4) is the summation of the pre-specified kernel matrix K and a rank-one matrix, as summarized in the following theorem:

Theorem 1 The optimal G∗ to the inner minimization problem of (4) can be analytically solved as

Proof: The inner minimization problem of (4) can be written as

which is convex in G. We first remove the constraint G ⪰ 0, and compute the optimal G∗ for (6). It is easy to show that the optimal G∗ satisfies

which leads to (5). Meanwhile, it is clear that G∗ ⪰ 0. Thus, the resulting G∗ is exactly the optimal solution to (6). □

Based on Theorem 1, the optimization problem in (4) can be reformulated as

which is a convex quartic problem. We propose to solve this convex quartic problem by level method in Section 3, and constrained gradient descent in Section 4.

2.2 Discussion

It follows from (5) that when ρ tends to +∞, G∗ approaches K, that is, the learnt optimal neighborhood kernel matrix reduces to the pre-specified kernel matrix.

The difference between the optimal neighborhood kernel matrix and the pre-specified kernel matrix is a rank-one matrix E = 1 2ρ (Yα)(Yα)T . Unlike K that is determined by the sample features and the pre-specified kernel κ(., .), E is specified by the computed α and the labels of the samples. In addition, denoting ˜ φ(xi) = (φ(xi), yiαi √2ρ ) for the training samples and ˜ φ(x) = (φ(x), 0) for the sample with unknown label, we can write the learnt optimal neighborhood kernel as ˜κ(x, z) = ˜ φ(x)T ˜ φ(z).

3 Optimization by the Level Method

The problem in (7) is convex in α so that a global optimum can be found. In this section, we show how this problem can be solved by the level method, which is a cutting plane method from the bundle family [Nemirovski, 1994].

Given a convex function f(α), the basic idea behind a cutting plane method is that, f(α) can be lower-bounded by its first-order Taylor approximation, i.e.,

where f ′ ( ˜α) is the derivative of f(α) at point ˜α. Assume we have obtained a set of α's in the set Ω = {α1,α2, . . . ,αi}. We define the following model:

which is a piecewise linear convex function, and underestimates the objective as fi(α) ≤ f(α).

Given the computed α's in Ω = {α1, . . . ,αi}, in the i-th step, the level method calculates αi+1 by sequentially performing the following three steps:

1. Compute the upper-bound fi = min1≤j≤i f(αj), and the lower-bound f i = fi(ui), where

can be solved by linear programming. Obviously, fi and f i respectively constitute the upper- and lower- bounds of minα∈P f(α), with the gap given by Δi = f i - f i .

2) Calculate the level

whose value is between fi and f i , with 0 < τ < 1 (usually set to 0.5) controlling the tradeoff. Moreover, we can form the level set Gi = {α|α ∈ P, fi(α) ≤ li}.

3) Project αi to the level set Gi to get the new solution:

which can be solved by quadratic programming [Boyd and Vandenberghe, 2004].

Algorithm 1 Level Method

Input: f(.), ε, τ, α0 ∈ P

Output: α

• 1: Initialize α1 = α0
• 2: for i = 1 to . . . do
• 3: Compute the lower-bound f i , the upper-bound fi, and the gap Δi
• 4: Calculate the level li, and form the level set Gi
• 5: Compute the new solution αi+1 = πGi (αi)
• 6: if Δi ≤ ε then
• 7: α = argmin1≤j≤i+1 f(αj), break
• 8: end if
• 9: end for

The pseudo code for solving the convex quartic problem by the level method is summarized in Algorithm 1. We have the following convergence property for the level method:

Theorem 2 [Nemirovski, 1994, Chapter 8, pages 135-137] Let L(f) < ∞ be the Lipschitz constant of the convex and Lipschitz continuous function f(α) on the convex domain P , i.e., |f(α) - f(β)| ≤ L(f)‖α-β‖2,α, β ∈ P. Let D(P) be the diameter of the bounded convex domain P , and c(τ ) = 1 (1-τ)2 τ(2-τ) . Then the gaps Δi's converge to zero, and for any ε > 0, if N ≥ c(τ)L2 (f )D2 (P )ε-2 , we have ΔN ≤ ε.

Theorem 2 shows that the level method has a theoretical complexity estimate O(ε-2 ), which is better than Kelly bundle method's O(ε-3 ). Moreover, as pointed out in [Nemirovski, 1994], the level method has a good practical performance, with an empirical O(n log( 1ε )) complexity.

4 Constrained Gradient Descent with Efficient Euclidean Projections

Although the level method has a good practical convergence property, it needs to solve at each iteration the linear programming problem (8) and the quadratic programming problem (10), whose computational costs increase with an increasing number of constraints during the iterations. In this section, we propose a non-memory-based method, "constrained gradient descent", equipped with an efficient Euclidean projection onto the domain P, with the favorable property that the Euclidean projection can be solved in linear time. Specifically, we introduce the gradient mapping for solving the constrained optimization problem in Section 4.1, cast the Euclidean projection onto the domain P as a zero finding problem in Section 4.2, and present the constrained gradient descent method for solving the problem (7) in Section 4.3.

4.1 Gradient Mapping

To deal with constrained optimization problems, we construct the gradient mapping, which acts a similar role as the gradient in unconstrained optimization. Let γ > 0. We define

which is the tangent line of f(.) at x, regularized by the squared distance between y and x. Minimizing fγ,x(y) in the domain P is the problem of Euclidean projections onto P :

We shall discuss the efficient computation of (11) in the next subsection. We call

the "gradient mapping" of f(.) on P . From (12), we have

which shows that, πP(x - 1 γ f ′ (x)) can be viewed as the result of the "gradient" step in the anti-direction of the gradient mapping p(γ, x) with stepsize 1 γ .

We say that γ is "appropriate" for x, if

It has been shown in [Nemirovski, 1994] that, γ is always appropriate for any x ∈ P , if γ ≥ Lf , the Lipschitz gradient of f(.). Moreover, when γ is appropriate for x, ∀y ∈ P , we have [Nemirovski, 1994, chapter 11, page 174]:

which is key to the convergence rate proof.

4.2 Casting Euclidean Projection onto the Domain P as a Zero Finding Problem

In this subsection, we discuss the efficient computation of the Euclidean projection onto P = {x|0 ≤ x ≤ Ce, xT y = 0}:

Our methodology is to cast (15) as a zero finding problem.

Introducing the Lagrangian variables λ, μ and ν respectively for the constraints xT y = 0, 0 ≤ x and x ≤ Ce, we can write the Lagrangian of (15) as

Let x⋆ be the primal optimal point, and λ⋆ , μ⋆ and ν ⋆ be the dual optimal points. We have the following results according to the KKT conditions [Boyd and Vandenberghe, 2004]:

From (16)-(19), we have

An analysis is given as follows. When 0 < vi-λ⋆yi < C, we have μ⋆i = ν ⋆ i = 0 (which leads to x⋆i = vi - λ⋆ yi), because if μ⋆i ̸= 0, we have x⋆i = 0 from (19), ν ⋆ i > 0 from (17) and x⋆i = C from (19), leading to a contradiction; similar contradiction can be constructed by assuming ν ⋆ i ̸= 0. When vi - λ⋆yi ≥ C, we have x⋆i = C, because if x⋆i < C, we have ν ⋆ i > 0 from (17), and x⋆i = C from (19). Similarly, when vi - λ⋆ yi ≤ 0, we have x⋆i = 0.

Combining (20) and (21), we have

which shows that λ⋆ is the root of the function h(.). The following theorem shows that λ⋆ exists and is unique.

Theorem 3 Let yi ∈ {±1}, i = 1, 2, . . . , n and y has both negative and positive elements. The auxiliary function h(.) is a continuous and monotonically decreasing function, and it has a unique root.

Proof: It is easy to verify that, ∀i, hi(λ) = yi max(min(vi - λyi,C) is continuous and monotonically decreasing in the whole domain (-∞,+∞). ∀i ∈ {i|yi = 1}, hi(.) is strictly decreasing in the interval [vi - C, vi]; and ∀i ∈ {i|yi = -1}, hi(.) is strictly decreasing in the interval [-vi,C - vi]. Denote λ11 = mini∈{i|yi=1} vi - C, λ12 = maxi∈{i|yi=1} vi, λ21 = mini∈{i|yi=-1} -vi, λ22 = maxi∈{i|yi=-1} -vi + C, λ1 = min(λ11, λ21), and λ2 = max(λ12, λ22). It is easy to get that h(.) is strictly decreasing in the interval [λ1, λ2]. Moreover, ∀i ∈ {i|yi = 1}, we have that, if λ ≤ λ11, hi(λ) = C; and if λ ≥ λ12, hi(λ12) = 0. Similarly, ∀i ∈ {i|yi = -1}, we have that, if λ ≤ λ21, hi(λ) = 0; and if λ ≥ λ22, hi(λ12) = -C. Therefore, we have h(λ1) > 0 and h(λ2) < 0. According to the Intermediate Value Theorem, h(.) has a unique root lying in the interval [λ1, λ2]. □

To solve the zero finding problem (22), we can make use of bisection, which produces a sequence of intervals of uncertainty with strictly decreasing lengths (decreased by a factor of two). For given v, y and C, λ1 and λ2 in Theorem 3 are fixed, and the number of bisection iterations is at most log2 ( λ2-λ1 δ ) under the machine precision δ. Therefore, we have a linear complexity for solving the Euclidean projection problem (15) by bisection. However, the efficiency of bisection is independent of the function h(.), and it cannot be improved even when h(.) is "well-behaved", since bisection only utilizes the signs of h(.) at the two boundary points, but not their values. To improve the efficiency, we follow the work of [Dekker, 1969] in using the interpolation methods that have a better local convergence rate than bisection.

4.3 Constrained Gradient Descent

Constrained Gradient Descent operates in a similar way to the classical gradient descent, and we present the pseudo code in Algorithm 2. In Step 3, we perform a line search for finding γk that is appropriate for xk, according to the ArmijoGoldstein rule [Nemirovski, 1994]. In Step 4, we compute the new approximate solution xk+1 = πP (xk1 γk f ′ (xk)) by the zero finding method proposed in Section 4.2. In Step 5, we terminate the algorithm, once the relative change in the approximate solution xk is not larger than the pre-specified parameter tol. Making use of (14) and following the similar technique as the one given in [Nemirovski, 1994, Chapter 10], we can establish the following convergence rate for constrained gradient descent.

Theorem 4 Let Lf be the Lipschitz gradient of f(.), i.e. ‖f ′(x) - f ′ (y)‖ ≤ Lf ‖x-y‖, ∀x, y ∈ P , and let Rf be the Euclidean distance from the starting point x0 to the optimal solution set of f(.). After N iterations, we have

5 Experiments

We conduct experiments on the following UCI data sets 1 : Liver Disorder, Ionosphere, Pima Indians Diabetes, Hill Valley, a2a, and a4a. Table 1 summarizes the characteris-

1 http://archive.ics.uci.edu/ml/

Algorithm 2 Constrained Gradient Descent Input: f(.), x0 ∈ P, tol Output: x 1: Initialize x1 = x0 2: for k = 1 to . . . do 3: Find γk that is appropriate for xk, i.e., satisfying (13) 4: Compute xk+1 = πP (xk - 1 γk f ′ (xk )) 5: if ‖xk+1 - xk‖ ≤ tol × max(‖xk‖, 1) then 6: x = xk+1, break 7: end if 8: end for

tics of the data sets. In this study, we generate the prespecified kernel by the radius basis kernel: κ(xi, xj) = exp(-‖xi - xj‖22/σ), and try three different values for σ: σk = 4-1+k ∑ni=1 ∑nj=1 ‖xi - xj‖22/n2 , k = 1, 2, 3. We implement the level method, by using the general solver MOSEK 2 for solving (8) and (10); and implement the constrained gradient descent with the efficient Euclidean projection proposed in Section 4.2. All experiments were carried out on an Intel (R) Core (TM)2 Duo 3.00GHZ processor.

Convergence Evaluation We set C = 1, ρ = 100, σ = σ2, and explore the convergence of the proposed level method (Algorithm 1) and constrained gradient descent (Algorithm 2) for learning the optimal neighborhood kernel matrix on data sets: Liver Disorder and a4a. We use all the samples in this task, and thus the learnt kernel matrices are of sizes 345 × 345 and 4781 × 4781, respectively. We report the results in Fig. 1, fromwhich we can observe that: (1) both level method and constrained gradient descent converge fast; and (2) level method consumes fewer iterations than constrained gradient descent, due to the usage of the whole "prehistory" for generating the approximate solution.

Time Efficiency We report the computational time for solving (7) on the six data sets (including all the samples) in the left figure of Fig. 2, from which we can observe that: (1) level method works relatively well; and (2) constrained gradient descent is quite efficient. It takes about 0.5 seconds to learn a kernel matrix of size 1212 × 1212. Compared to the level method, the constrained gradient descent is much more efficient, since each of its iteration takes less time. To explore the efficiency of the Euclidean projection method proposed in Section 4.2, we compare it with the quadratic programming provided in MOSEK, and report the results in the right figure

2 http://www.mosek.com/

n

Figure 2: Computational time. The left figure reports the computational time for learning the kernel matrix on the six data sets; and the right figure shows the total computational time for solving 100 independent Euclidean projection problems by the method proposed in Section 4.2 and the quadratic programming provided by MOSEK.

of Fig. 2, where the x-axis corresponds to n, the problem size, and the y-axis denotes the total computational time for solving 100 independent problems. We can observe from the figure that, the proposed Euclidean projection method is about 100 times faster than the quadratic programming of MOSEK for this task.

Classification Performance In exploring the classification performance of the proposed method, we randomly choose 20% from each data set for training, and the remaining for testing, and report the classification performance by averaging over 20 runs. We choose the values of the parameters C and ρ through cross-validation, and report the recognition results in Table 2. We can observe from the table that the learnt optimal neighborhood kernel can yield better classification performance than the pre-specified kernel. In many cases, the improvement is significant. This demonstrates the effectiveness of the proposed optimal neighborhood kernel.

Sensitivity Study We perform the sensitivity study in terms of ρ. Fig. 3 shows the classification performance of the learnt

Liver Ionosphere Diabetes Hill Valley a2a a4a

500

1000

2000

4000

8000

16000

optimal neighborhood kernel under different ρ's on two data sets: Liver Disorder and Pima Indians Diabetes. We can observe from this figure that: (1) when ρ is quite small, the classification performance is in general not very good, since the rank-one matrix E = 1 2ρ (Yα)(Yα)T dominates the learnt kernel matrix G = K + E; (2) for a relatively wide range of values for ρ, the learnt optimal neighborhood kernel achieves better classification performance than the pre-specified kernel; and (3) when ρ is large, the classification performance of the optimal neighborhood kernel is almost identical to that of the pre-specified kernel, since the effect of E is diminished with a large value of ρ.

6 Conclusion

In this paper, we propose to learn the optimal neighborhood kernel matrix from the pre-specified positive semidefinite kernel matrix for improved classification performance. The optimal neighborhood kernel matrix is shown to be the summation of the pre-specified kernel matrix and a rank-one matrix. We formulate the problem of learning the optimal neighborhood kernel as a convex quartic problem, and propose to solve it by the level method and the constrained gradient descent. To solve the constrained problem, we propose an efficient Euclidean projection onto the domain P . Empirical results on a collection of benchmark data sets demonstrate the efficiency and effectiveness of the proposed algorithms.

We have focused on the learning from a single kernel matrix in this paper, and we plan to extend the current work to multiple kernel learning, where an optimal kernel matrix is learnt from a given collection of kernel matrices. We plan to accelerate the proposed constrained gradient descent by using the Nesterov's method [Nemirovski, 1994], and apply the proposed efficient Euclidean projection to solving other related constrained optimization problems. We also plan to extend the proposed algorithm to the transductive setting.

Acknowledgments

This work was supported by NSF IIS-0612069, IIS-0812551, CCF-0811790, NIH R01-HG002516, NGA HM1582-08-10016, NSF of China 60773061, and NSF of Jiangsu Province BK2008381.