Kernel Conjugate Gradient for Fast Kernel Machines

1 Introduction

Kernel methods, in their various incarnations (e.g. Gaussian Processes (GPs), Support Vector Machines (SVMs), Kernel Logistic Regression (KLR)) have recently become a preferred approach to non-parametric machine learning and statistics. They enjoy this status because of their conceptual clarity, strong empirical performance, and theoretical foundations.

The primary drawback to kernel methods is their computational complexity. GPs require the inversion of an n × n (covariance/kernel) matrix, implying a running time of O(n3 ), where n is the size of the training set. SVMs require similar computation to solve the convex program, although intense research has gone into fast, specialized approximations [Sch¨olkopf & Smola, 2002].

State-of-the-art approaches to kernel learning revolve largely around two techniques: iterative optimization algorithms, and learning by representing the solution with only a subset of the original data points. Our algorithm applies most directly to the former line of research, although we address the latter in the conclusions.

We propose a novel variant of the conjugate gradient algorithm, Kernel Conjugate Gradient (KCG), designed to speed up learning for kernel machines with differentiable loss functions (e.g. Gaussian Process mean inference, Kernel Logistic Regression). This algorithm is motivated by the understanding that all gradient-based methods rely, at least implicitly, on a particular metric or inner product [Sch¨olkopf & Smola, 2002]. It is natural in kernel learning problems for the algorithm to inherit the metric on functions that the kernel provides. In Section 2.4 we show that such an approach can be interpreted as a Riemannian metric method similar to Amari's Natural Gradient [Amari & Nagaoka, 2000], although the kernel metric is much less expensive to compute.

In Section 5, we establish an upper bound on the number of iterations for KCG that indicates it should require fewer than the square root of the number of iterations that standard conjugate gradient requires. The algorithm is simple to implement, requires no more computation per iteration than standard approaches, and is well motivated by Reproducing Kernel Hilbert Space (RKHS) theory. In practice, for various differentiable kernel learning problems, we find KCG consistently, and significantly, outperforms existing techniques.

Recent research has demonstrated the benefits of spacepartitioning data-structures to speed up certain kernel machines [Shen et al., 2006; Gray & Moore, 2001]. We show that these techniques work well with KCG as they reduce the dominant computational burdens in training: RKHS function evaluations and inner product computations.

2 Preliminaries

Below we briefly review the theory of kernel machines in terms of Reproducing Kernel Hilbert Spaces. We then describe the regularized risk functionals we are interested in optimizing during learning, and finish by reviewing the functional gradient, a generalization of the standard gradient to inner product spaces of functions.

2.1 Reproducing Kernel Hilbert Spaces

An RKHS of functions Hk is a complete inner product space, known as a Hilbert space, that arises from the completion of a set of basis functions BX k = {k(x, .) | x ∈ X}, where k : X ×X → R is a symmetric, positive-definite kernel function, and X is a continuous domain sometimes known as an index set. A common kernel, and one used exclusively in our experiments, is the exponential Radial Basis Function (RBF)

k(x, x′ ) = e -‖x-x′ ‖2 2σ2 . The RKHS inner product between two functions f = ∑i αik(xi, .) and g = ∑j βjk(xj, .) is defined as

Central to the idea of an RKHS is the reproducing property which follows directly from the above definition. It states that the basis functions k(x, .) ∈ BX k are representers of evaluation. Formally, for all f ∈ Hk and x ∈ X, 〈f, k(x, .)〉Hk = f(x). When the basis functions are normalized, this means the evaluation of f at x is the scalar projection of f onto k(x, .). Note that there exists a simple mapping φ : X → BX k between the domain X and the RKHS basis BX k defined by the kernel as φ(x) = k(x, .). It follows that for any x, x′ ∈ X

A complete overview of these concepts can be found in [Aronszajn, 1950].

A fundamental result first proven in [Kimeldorf &Wahba, 1971] and then generalized in [Sch¨olkopf et al., 2001] is the Representer Theorem, which makes possible the direct minimization of a particular class of functionals. It states that given a subset D = {(xi, yi)}ni=1 ⊂ X × R (i.e. the dataset), the minimizer of a functional with the form F [f ] = c((x1, y1, f(x1)), . . . , (xn, yn, f(xn)))+g(‖f‖2Hk ), where c : X × R2 → R is arbitrary and g : [0,∞] → R is strictly monotonically increasing, must have the form ˜ f = ∑xi∈D αik(xi, .).

2.2 Regularized Risk Functionals

An important class of functionals, common in machine learning, for which the Representer Theorem holds is the regularized risk functional [Sch¨olkopf & Smola, 2002]:

These functionals combine a data-dependent risk term, with a "prior" (or regularization) term that controls the complexity of the solution by penalizing the norm of the function f in the RKHS. Our work focuses on the case where l is differentiable in its third argument; many important kernel machines have this property. For instance, we can write Kernel Logistic Regression [Zhu & Hastie, 2001] in this form. Let y ∈ {-1, 1}. Then

Similarly, we may write the popular Gaussian Process Regression (for the mean function) and the Regularized Least Squares Classification algorithm in this form 1 :

1 These two differently motivated algorithms optimize the same functional.

Numerous other important examples occur in the literature [Sch¨olkopf & Smola, 2002]. We focus on these, two of the best known kernel algorithms other than the nondifferentiable Support Vector Machine, as they are particularly important tools in machine learning for classification and regression.

2.3 The Functional Gradient

For large-scale problems and problems for which direct solution is inapplicable, the most popular approach to optimizing these functionals is through the use of iterative, gradientbased techniques. Gradient-based algorithms, such as steepest descent, are traditionally defined with respect to the gradient arising from the Euclidean inner product between parameter vectors. There are many ways in which a given regularized risk functional can be parameterized, though, and each gives rise to a different parameter gradient. It is intuitively natural to follow the gradient defined uniquely by the RKHS inner product. We review some basic concepts behind the functional gradient below.

The functional gradient may be defined implicitly as the linear term of the change in a function due to a small perturbation ϵ in its input [Mason et al., 1999]

Following this definition using the RKHS inner product 〈., .〉Hk allows us to define the kernel gradient of a regularized risk functional [Sch¨olkopf & Smola, 2002]. We use three basic formulas that can be easily derived from this definition:

1. Gradient of the evaluation functional. Define Fx : Hk → R as Fx[f] = f(x) = ∑i αik(xi, x). Then ∇kFx[f] = k(x, .).
2. Gradient of the square RKHS norm. Define F〈.,.〉 : Hk → R as F〈.,.〉[f] = ‖f‖2Hk = 〈f, f〉Hk . Then ∇kF〈.,.〉[f ] = 2f .
3. Chain rule. Let g : R → R be a differentiable function, and let F : Hk → R be an arbitrary differentiable functional. Then ∇kg(F [f ]) = g′ (F [f ])∇kF [f ].

A straight forward application of the above formulas brings us to the following result. The kernel gradient of a regularized risk functional (Equation 2) is

where we again use αi to denote the parameters of the expansion of f in terms of the kernels k(xi, .).

Equation 5 allows us to easily find the kernel gradient of the functionals we described above and use in the remainder of this paper:

Kernel Logistic Regression (KLR):

Regularized Least Squares (RLS), Gaussian Process:

2.4 The Kernel Gradient and RiemannianMetrics

Above we described the kernel gradient as a function; that is as a linear combination of basis functions. The kernel gradient as in Equation 5 demonstrates a property similar to the Representer Theorem: namely that ∇kF[f] = ∑ni=1 γik(xi, .) for appropriate γi ∈ R. In other words, the kernel gradient of F is represented in the finite-dimensional subspace SD = span{BD k } of Hk. A gradient descent type method through SD then amounts to modifying the coefficients αi of the current function f by γi. That is, we can understand the kernel gradient as modifying parameters,

just as a standard gradient descent algorithm would. The difference is that the coefficients γi for the kernel gradient are not the same as those of the parameter gradient. One can verify that they differ by ∇αF[f] = Kγ where K is the kernel matrix and γ is the vector of coefficients.

We can derive this relation in another way. Starting from the parameter gradient, we define a Riemannian metric [Hassani, 1998] on the space of parameters. This defines our notion of size in the space of parameters. We then consider an alternate definition of the gradient, as the direction of steepest ascent for a "small" change in coefficients α:

It can be shown that taking ‖γ‖2α as γ T γ gives the vanilla parameter gradient ∇αF , while defining the norm with respect to the RKHS inner product ‖γ‖2Hk = ∑i,j γiγjk(xi, xj) = γ T Kγ gives the functional gradient coefficients ∇kF = K-1∇αF. [Hassani, 1998]

This interpretation of the kernel gradient makes connections with other metric methods more clear. For instance, Amari [Amari & Nagaoka, 2000] considers the use of a metric derived from information geometry that leads to the "natural gradient". Such algorithms are applicable here as well since we can compute the metric for the probabilistic models given by Gaussian Processes or KLR. Unfortunately, computing the natural gradient in these cases is very expensive: for instance, in Gaussian Processes it is as expensive as inverting the kernel matrix, the very computational difficulty we are striving to avoid. By contrast, computing the kernel gradient is very cheap: cheaper in fact then the standard parameter gradient. 2

2 Furthermore, in practice we often find deriving the kernel gradient using the functional gradient rules easier than deriving the parameter gradient.

Algorithm 1 Kernel Conjugate Gradient

• 12: end procedure

3 The Kernel Conjugate Gradient Algorithm

Both in theory and in practice it is understood that conjugate gradient (CG) methods outperform standard steepest descent procedures [Ashby et al., 1990]. These techniques have been used profusely throughout machine learning, in particular, for regularized risk minimization and kernel matrix inversion [Gibbs, 1997; Sch¨olkopf & Smola, 2002].

In this section, we present an algorithm we term Kernel Conjugate Gradient (KCG) that takes advantage of conjugate direction search while utilizing the RKHS inner product 〈f, g〉Hk = αT Kβ. Algorithm 1 gives the general (nonlinear) KCG algorithm in Polak-Ribi`ere form [Ashby et al., 1990]. In essence, Algorithm 1 comes directly from conjugate gradient by replacing all gradients by their functional equivalents and replacing Euclidean inner products with an RKHS inner product.

Note that the computational complexity per iteration of KCG is essentially identical to that of the more conventional Parameter Conjugate Gradient (PCG) algorithm. Intuitively, while the kernel inner product takes time O(n2 ) compared to the O(n) vanilla inner product used by PCG, KCG is correspondingly more efficient in the gradient computation since ∇αF[f] = Kγ, where ∇kF[f] = ∑i γik(xi, .). It is possible in the case of RLS to step through an iteration of each algorithm and show that the number of operations is equivalent.

We emphasize that despite its somewhat involved derivation, the implementation of this algorithm is just a simple extension of PCG. The differences amount to only a change of inner product αT β → ∑i,j αiβjk(xi, xj) = αT Kβ, and a different, though in some ways simpler, gradient computation. We also point out that the line optimization (step 6) can be solved in closed-form in the case of quadratic risk functionals (e.g. RLS). For starting point f = ∑i αik(xi, .) and search direction h = ∑i γik(xi, .) we have

where A is the Hessian of the quadratic functional when parameterized by α. Note that this formula differs from that derived under the parameter gradient (-αT γ/γ T Aγ) only in the numerator's inner product, as is a common theme throughout this algorithm. The theoretical and experimental results given below suggest that there is little reason why one should prefer PCG to KCG in most (differentiable) kernel algorithms.

4 Experimental Results - Kernel Conjugate Gradient

We bench-marked KCG against PCG for both classification and regression tasks. In all cases, KCG significantly out performed PCG.

Our first test was performed using KLR on the USPS dataset (with a training/test size of 7291/2007) for the common one-vs-all task of recognizing the digit 4. We used a length scale hyperparameter σ = 5 as was used in [Rifkin et al., 2003] for RLS classification, and a regularization constant λ = 0. Figure 1 summarizes the results in log scale.

Second, we used RLS for both regression and classification using the Abalone and Soil datasets in addition to the USPS dataset. The Abalone dataset 3 consisted of 3133 training examples and 1044 test examples in 8 attributes. Our experimental setup was equivalent to that in [Smola & Sch¨olkopf, 2000]. The Soil dataset contained three-dimensional examples of soil pH levels in areas of Honduras partitioned into a training set of size 1709 and a test set of size 383. The latter dataset and corresponding hyperparameters (λ = 0.1514 and σ = 0.296283) were provided by [Gonzalez, 2005]. Again, the results are summarized in Figure 1.

For RLS, there exists a quadratic

that can provide a lower bound to the regularized risk [Gibbs, 1997; Sch¨olkopf & Smola, 2002]. As theory suggests (see below) the upper bound under KCG converges comparably with the lower bound, while the upper bound under PCG lags considerably behind. This implies faster convergence under a gap termination criterion [Sch¨olkopf & Smola, 2002].

The right-most plots of Figure 1 contrast (in iterations, equivalent to multiples of wall-clock time) the speeds of PCG and KCG for both RLS and KLR. We plot of the number of iterations of each to reach the same level of performance in terms of loss on the data-set. Plots are terminated when convergence is achieved as measured with the gap termination criterion[Sch¨olkopf & Smola, 2002], and hyper-parameters were chosen on a hold-out set. These figures confirm what analysis in the next section suggests: it takes more than the square of the amount of time to achieve the same level of performance using PCG as it does with KCG. Finally, in Table 2 we directly compare the number of iterations needed to achieve convergence on a number of data sets, all taken again from the UCI data repository, for both RLS and KLR. Averaged over the datasets, KCG is 54 times faster than the standard conjugate gradient approach.

3 See UCI repository: http://www.ics.uci.edu/

∼mlearn/MLRepository.html

5 KCGAnalysis

We derived the Kernel Conjugate Gradient algorithm from a normative point of view arguing that 〈f, g〉Hk defined the natural notion of inner product in the RKHS and hence for the optimization procedure as well. The strong empirical performance of KCG noted in the previous section, while in some sense not surprising given we are using the "correct" inner product, deserves analysis. We examine here the linear case (as in RLS) where the analysis is more transparent, although presumably similar results hold near the optima of non-linear risk functionals.

We note a classic bound on the error reduction of CG (see [Luenberger, 2003]),

where i is the iteration number, A is the Hessian of the quadratic form with condition number κ, ‖x‖A = xT Ax is a norm on x, and ei = xi - x∗ . Loosely speaking, this gives a running time complexity of O(√κ) for PCG.

We start by analyzing the effective condition number of KCG. As with essentially all variants of CG, the algorithm's dynamics can be described in terms of a preconditioning to the spectrum of the Hessian [Ashby et al., 1990]. It can be verified by inspection of the algorithm, that KCG is equivalent to an implicitly preconditioned conjugate gradient algorithm with preconditioner K[Ashby et al., 1990]. The following theorem relates the running time of PCG to KCG in light of the bound given above.

Theorem. Let κPCG be the condition number of Rrls (Equation 4), and let κK be the condition number of the kernel matrix K = [k(xi, xj)]i,j. Then the condition number κKCG resulting from preconditioning the RLS risk functional by K has the relation κPCG = κKκKCG.

Proof. Let σ1 ≥ σ2 ≥ . . . ≥ σn be the eigenvalues of K. The condition number of K is then κK = σ1/σn. The Hessian of Rrls is A = KT K + λK and has eigenvalues σ2 i + λσi = σi(σi + λ), given in terms of the eigenvalues of K. This implies

Since K is symmetric, positive-definite, the preconditioned Hessian becomes K-1A = K-1 (KT K + λK) = K + λI, with corresponding eigenvalues σi + λ. Thus, κPCG = κKκKCG. ✷

The condition number κK of K is typically very large. In particular, as the regularization constant decreases, the asymptotic bound on the convergence of PCG approaches the square of the bound on KCG. Alternatively, as the regularization constant increases, κKCG approaches 1 implying an O(1) convergence bound for KCG, while the convergence bound of PCG remains bounded below by O(κK). We would thus expect a number of iterations for KCG that is dramatically less than that of PCG.

It is informative to note that the decrease in computational complexity from PCG to KCG (O(κ1/2 ) to O(κ1/4 )) is at least the amount we see from steepest descent(O(κ) to PCG O(κ1/2 )) [Luenberger, 2003].

6 Tree-Augmented Algorithms

For many stationary kernels, it is often the case that the majority of the basis functions in BD k are nearly orthogonal. This often stems from a relationship between the degree of orthogonality of two basis functions k(x, .), k(x′ , .) ∈ BD k and the Euclidean distance between x and x′ . This is often the case when using the exponential RBF kernel given above, in which the orthogonality between two basis functions increases exponentially with the square Euclidean distance between the two points ‖x - x′ ‖2 .

Previous work has shown how the evaluation of RKHS functions f(x) = ∑i αik(xi, x) can be made fast when this holds using N-body type algorithms [Gray & Moore, 2001]. Intuitively, the idea is to store the training data in a spacepartitioning tree, such as a KD-tree [Moore, 1990] as was used in our experiments below, and recursively descend the tree during evaluation, pruning negligible contributions.

The maximum and minimum impact of each set of pruned points can be easily calculated resulting in upper and lower bounds on the evaluation error. We demonstrate how such data-structures and algorithms can be used to reduce the per iteration O(n2 ) computational cost of both KCG and PCG during learning as well as evaluation.

The inner loop computational bottleneck of KCG is in evaluating functions and calculating the kernel inner product. If we rewrite the RKHS inner product between f = ∑i αik(xi, .), g = ∑i βik(xi, .) as 〈f, g〉Hk = ∑i,j αiβjk(xi, xj) = ∑i αig(xi), then reducing the computational complexity of RKHS function evaluations will simultaneously encompass both of these bottlenecks. Similarly, the iteration complexity of PCG

is dominated by the computation of the parameter gradient. We can rewrite the parameter gradient as ∇αF[f] = [∇kF[f](x1),∇kF[f](x2), . . . ,∇kF[f](xn)]T (see 2.4), reducing the complexity to that of finding ∇kF[f] ∈ Hk and evaluating it n times. As was the case with the KCG algorithm without trees, the tradeoff still balances out so that the per iteration complexity is essentially equivalent between KCG and PCG using tree-augmented function evaluation.

Noting [f(x1), . . . , f(xn)]T = Kα suggests that the closed form quadratic line minimization for both the upper and lower bounds of RLS under either KCG or PCG can easily be augmented as well by expanding the Hessians Au = KT K+λK in the case of the upper bound and Al = K+λI in the case of the lower bound. This was used in the treeaugmented RLS experiments described below.

7 Experimental Results - Tree-Augmented Algorithms

For these experiments, in addition to using the Soil dataset described in section 4, we performed large scale RLS regressions using a tree-augmented KCG on variable sized subsets of a PugetSound elevation map 4 using hyperparameters λ = 0.1/n and σ2 = kN/(nπ), where n is the size of the training set, and N is the size of the entire height map. In this case, we chose N = 500 × 500 = 250, 000, and k = 15. The largest resulting datasets were on the order of 100,000 points. It should be noted that the n¨aive implementation in this case did not cache kernel evaluations in a kernel matrix as such matrices for datasets above O(15, 000) points proved intractable for the machine on which the experiments were performed.

4 http://www.cc.gatech.edu/projects/large models/

Figure 1 shows that tree-augmentation significantly outperforms the na¨ıve algorithm. Extrapolating from the plot on the right, trees make possible accurate kernel learning on very large datasets without requiring explicit subset selection techniques.

8 Conclusions and FutureWork

We have demonstrated that the novel gradient method, Kernel Conjugate Gradient, can dramatically improve learning speed for differentiable kernel machines. Furthermore, we have shown how this can be understood as a very efficient preconditioning that naturally derives from the inner product on functions defined by the kernel. In practice, for various differentiable kernel learning problems, we find KCG consistently, and significantly, outperforms existing techniques. We emphasize that the algorithm is simple to implement and requires no more computation per iteration than standard approaches.

Further, we demonstrated that space-partitioning datastructures, also developed by other authors [Shen et al., 2006], for optimizing Gaussian Processes extend naturally to other kernel methods. We find that this approach meshes well with the KCG algorithm, by significantly speeding up the inner loop computations of functions and inner products.

While conjugate gradient is a powerful optimization algorithm, there are other approaches, like Limited-Memory Quasi-Newton methods [Luenberger, 2003] that may also be derived in similar ways in terms of the kernel inner product. These algorithms have proved to be practical when using the Euclidean inner product; we expect they would also gain the benefits of preconditioning that KCG enjoys.

Finally, very large scale kernel applications, it seems, will invariably need to rely on sparse representation techniques that do not use kernels at all of the data points. Nearly all of these methods require the efficient solution of large kernel problems using an iterative approximation. It is natural to explore how Kernel Conjugate Gradient can speed up the expensive inner loop of these approximation procedures.

Acknowledgements

The authors gratefully acknowledge the partial support of this research by the DARPA Learning for Locomotion contract.