Scalable Kernel Methods via Doubly Stochastic Gradients

2014·arXiv

Abstract

1 Introduction

The general perception is that kernel methods are not scalable. When it comes to large-scale nonlinear learning problems, the methods of choice so far are neural nets where theoretical understanding remains incomplete. Are kernel methods really not scalable? Or is it simply because we have not tried hard enough, while neural nets have exploited sophisticated design of feature architectures, virtual example generation for dealing with invariance, stochastic gradient descent for efficient training, and GPUs for further speedup?

A bottleneck in scaling up kernel methods is the storage and computation of the kernel matrix, K, which is usually dense. Storing the matrix requires ) space, and computing it takes ) operations, where n is the number of data points and d is the dimension. There have been many great attempts to scale up kernel methods, including efforts from numerical linear algebra, functional analysis, and numerical optimization perspectives.

A common numerical linear algebra approach is to approximate the kernel matrix using low-rank factors, , with and rank . This low-rank approximation usually requires + nrd) operations, and then subsequent kernel algorithms can directly operate on A. Many works, such as Greedy basis selection techniques [1], Nystr¨om approximation [2] and incomplete Cholesky decomposition [3], all followed this strategy. In practice, one observes that kernel methods with approximated kernel matrices often result in a few percentage of losses in performance. In fact, without further assumption on the regularity of the kernel matrix, the generalization ability after low-rank approximation is typically of the order + 1) [4, 5], which implies that the rank needs to be nearly linear in the number of data points! Thus, in order for kernel methods to achieve the best generalization ability, the low-rank approximation based approaches quickly become impractical for big datasets due to their + ) preprocessing time and ) memory requirement.

Random feature approximation is another popular approach for scaling up kernel methods [6, 7]. Instead of approximating the kernel matrix, the method directly approximates the kernel function using explicit feature maps. The advantage of this approach is that the random feature matrix for n data points can be computed in time O(nrd) using O(nr) memory, where r is the number of random features. Subsequent algorithms then only operate on an O(nr) matrix. Similar to low-rank kernel matrix approximation approach, the generalization ability of random feature approach is of the order ) [8, 9], which implies that the number of random features also needs to be O(n). Another common drawback of these two approaches is that it is not easy to adapt the solution from a small r to a large . Often one is interested in increasing the kernel matrix approximation rank or the number of random features to obtain a better generalization ability. Then special procedures need to be designed to reuse the solution obtained from a small r, which is not straightforward.

Another approach that addresses the scalability issue rises from optimization perspective. One general strategy is to solve the dual forms of kernel methods using coordinate or block-coordinate descent (e.g., [10, 11, 12]). By doing so, each iteration of the algorithm only incurs O(nrd) computation and O(nr) memory, where r is the size of the parameter block. A second strategy is to perform functional gradient descent by looking at a batch of data points at a time (e.g., [13, 15]). Thus, the computation and memory requirements are also O(nrd) and O(nr) respectively in each iteration, where r is the batch size. These approaches can easily change to a different r without restarting the optimization and has no loss in generalization ability since they do not approximate the kernel matrix or function. However, a serious drawback of these approaches is that, without further approximation, all support vectors need to be kept for testing, which can be as big as the entire training set! (e.g., kernel ridge regression and non-separable nonlinear classification problems.)

In summary, there exists a delicate trade-off between computation, memory and statistics if one wants to scale up kernel methods. Inspired by various previous efforts, we propose a simple yet general strategy to scale up many kernel methods using a novel concept called “doubly stochastic functional gradients”. Our method relies on the fact that most kernel methods can be expressed as convex optimization problems over functions in reproducing kernel Hilbert spaces (RKHS) and solved via functional gradient descent. Our algorithm proceeds by making two unbiased stochastic approximations to the functional gradient, one using random training points and the other one using random features associated with the kernel, and then descending using this noisy functional gradient. The key intuitions behind our algorithm originate from

(i) the property of stochastic gradient descent algorithm that as long as the stochastic gradient is unbiased, the convergence of the algorithm is guaranteed [16]; and

(ii) the property of pseudo-random number generators that the random samples can in fact be completely determined by an initial value (a seed).

We exploit these properties and enable kernel methods to achieve better balances between computation, memory and statistics. Our method interestingly combines kernel methods, functional analysis, stochastic optimization and algorithmic trick, and it possesses a number of desiderata:

Generality and simplicity. Our approach applies to many kernel methods, such as kernel ridge regression, support vector machines, logistic regression, two-sample test, and many different types of kernels, such as shift-invariant kernels, polynomial kernels, general inner product kernels, and so on. The algorithm can be summarized in just a few lines of code (Algorithm 1 and 2). For a different problem and kernel, we just need

to adapt the loss function and the random feature generator.

Flexibility. Different from previous uses of random features which typically prefix the number of features and then optimize over the feature weightings, our approach allows the number of random features, and hence the flexibility of the function class, to grow with the number of data points. This allows our method to be applicable to data streaming setting, which is not possible for previous random feature approach, and achieve the full potential of nonparametric methods.

Efficient computation. The key computation of our method is evaluating the doubly stochastic functional gradient, which involves the generation of the random features with specific random seeds and the evaluation of these random features on the small batch of data points. For iteration t, the computational complexity is O(td).

Small memory. The doubly stochasticity also allows us to avoid keeping the support vectors which becomes prohibitive in large-scale streaming setting. Instead, we just need to keep a small program for regenerating the random features, and sample previously used random feature according to pre-specified random seeds. For iteration t, the memory needed is O(t) independent of the dimension of the data.

Theoretical guarantees. We provide a novel and nontrivial analysis involving Hilbert space martingale and a newly proved recurrence relation, and show that the estimator produced by our algorithm, which might be outside of the RKHS, converges to the optimal RKHS function. More specifically, both in expectation and with high probability, our algorithm can estimate the optimal function in the RKHS in the rate of O(1/t), which are indeed optimal [16], and achieve a generalization bound of ). The variance of the random features, introduced during our second approximation to the functional gradient, only contributes additively to the constant in the final convergence rate. These results are the first of the kind in kernel method literature, which can be of independent interest.

Strong empirical performance. Our algorithm can readily scale kernel methods up to the regimes which are previously dominated by neural nets. We show that our method compares favorably to other scalable kernel methods in medium scale datasets, and to neural nets in big datasets such as 8 million handwritten digits from MNIST, 2.3 million materials from MolecularSpace, and 1 million photos from ImageNet using convolution features. Our results suggest that kernel methods, theoretically well-grounded methods, can potentially replace neural nets in many large scale real-world problems where nonparametric estimation are needed.

In the remainder, we will first introduce preliminaries on kernel methods and functional gradients. We will then describe our algorithm and provide both theoretical and empirical supports.

2 Duality between Kernels and Random Processes

Kernel methods owe their name to the use of kernel functions, ) : , which are symmetric positive definite (PD), meaning that for all n > 1, and , and , we have 0. There is an intriguing duality between kernels and stochastic processes which will play a crucial role in our later algorithm design. More specifically,

Theorem 1 (e.g.,[17]; [18]) If ) is a PD kernel, then there exists a set Ω, a measure P on Ω, and random feature ) : from ), such that

Essentially, the above integral representation relates the kernel function to a random process with measure ). Note that the integral representation may not be unique. For instance, the random process can be a Gaussian process on X with the sample function ), and ) is simply the covariance function between two point x and . If the kernel is also continuous and shift invariant, ) for , then the integral representation specializes into a form characterized by inverse Fourier transformation (e.g., [19, Theorem 6.6]),

Theorem 2 (Bochner) A continuous, real-valued, symmetric and shift-invariant function ) on is a PD kernel if and only if there is a finite non-negative measure ) on , such that

2 cos() cos(tribution on [02 cos(

For Gaussian RBF kernel, ), this yields a Gaussian distribution ) with density proportional to exp(2); for the Laplace kernel, this yields a Cauchy distribution; and for the Martern kernel, this yields the convolutions of the unit ball [20].

Similar representation where the explicit form of ) and ) are known can also be derived for rotation invariant kernel, ), using Fourier transformation on sphere [20]. For polynomial kernels, , a random tensor sketching approach can also be used [21]. Explicit random features have been designed for many other kernels, such as dot product kernel [33], additive/multiplicative class of homogeneous kernels [34], e.g., Hellinger’s, , Jensen-Shannon’s and Intersection kernel, as well as kernels on Abelian semigroups [35]. We summarized these kernels with their explicit features and associated densities in Table 1.

Instead of finding the random process ) and function ) given a kernel, one can go the reverse direction, and construct kernels from random processes and functions (e.g., [19]).

Theorem 3 If ) for a nonnegative measure ) on Ω and ) : , each component from ), then ) is a PD kernel.

For instance, ), where ) can be a random convolution of the input x parametrized by , or )], where ) denote the random feature for kernel ). The former random features define a hierachical kernel [45], and the latter random features induce a linear combination of multiple kernels. It is worth to note that the Hellinger’s, , Jensen-Shannon’s and Intersection kernels in [34] are special cases of multiple kernels combination. For simplicity, we assume following, and our algorithm is still applicable to .

Another important concept is the reproducing kernel Hilbert space (RKHS). An RKHS H on X is a Hilbert space of functions from X to R. H is an RKHS if and only if there exists a ) : such that and If such a ) exist, it is unique and it is a PD kernel. A function if and only if , and its norm is dominated by RKHS norm

3 Doubly Stochastic Functional Gradients

Many kernel methods can be written as convex optimizations over functions in the RKHS and solved using the functional gradient methods [13, 15]. Inspired by these previous works, we will introduce a novel concept called “doubly stochastic functional gradients” to address the scalability issue. Let l(u, y) be a scalar (potentially non-smooth) loss function convex of . Let the subgradient of l(u, y) with respect to u be ). Given a PD kernel ) and the associated RKHS H, many kernel methods try to find a function which solves the optimization problem

where 0 is a regularization parameter, ) is a non-increasing function of , and the data (x, y) follow a distribution P(x, y). The functional gradient ) is defined as the linear term in the change of the objective after we perturb f by in the direction of g, i.e.,

For instance, applying the above definition, we have ), and .

Stochastic functional gradient. Given a data point () and , the stochastic functional gradient of )] with respect to is

which is essentially a single data point approximation to the true functional gradient. Furthermore, for any , we have ). Inspired by the duality between kernel functions and random processes, we can make an additional approximation to the stochastic functional gradient using a random feature ) sampled according to ). More specifically,

Doubly stochastic functional gradient. Let ), then the doubly stochastic gradient of )] with respect to is

Figure 1: is the angle betwen ) and ) may be outside of H.

Note that the stochastic functional gradient ) is in RKHS H but ) may be outside H, since ) may be outside the RKHS. For instance, for the Gaussian RBF kernel, the random feature 2 cos() is outside the RKHS associated with the kernel function.

We emphasize that the source of randomness associated with the random feature is not present in the data, but artificially introduced by us. This is crucial for the development of our scalable algorithm in the next section. Meanwhile, it also creates additional challenges in the analysis of the algorithm which we will deal with carefully.

4 Doubly Stochastic Kernel Machines

The first key intuition behind our algorithm originates from the property of stochastic gradient descent algorithm that as long as the stochastic gradient is unbiased, the convergence of the algorithm is guaranteed [16]. In our algorithm, we will exploit this property and introduce two sources of randomness, one from data and another artificial, to scale up kernel methods.

The second key intuition behind our algorithm is that the random features used in the doubly stochastic functional gradients will be sampled according to pseudo-random number generators, where the sequences of apparently random samples can in fact be completely determined by an initial value (a seed). Although these random samples are not the “true” random sample in the purest sense of the word, however they suffice for our task in practice.

More specifically, our algorithm proceeds by making two unbiased stochastic approximation to the functional gradient in each iteration, and then descending using this noisy functional gradient. The overall algorithms for training and prediction is summarized in Algorithm 1 and 2. The training algorithm essentially just performs random feature sampling and doubly stochastic gradient evaluation, and maintains a collection of real number , which is computationally efficient and memory friendly. A crucial step in the algorithm is to sample the random features with “seed i”. The seeds have to be aligned between training and prediction, and with the corresponding obtained from each iteration. The learning rate in the algorithm needs to be chosen as O(1/t), as shown by our later analysis to achieve the best rate of convergence. For now, we assume that we have access to the data generating distribution P(x, y). This can be modified to sample uniformly randomly from a fixed dataset, without affecting the algorithm and the later convergence analysis. Let the sampled data and random feature parameters be and respectively after t iteration, the function obtained by Algorithm 1 is a simple additive form of the doubly stochastic functional gradients

where ) are deterministic values depending on the step sizes ) and regularization parameter . This simple form makes it easy for us to analyze its convergence.

We note that our algorithm can also take a mini-batch of points and random features at each step, and estimate an empirical covariance for preconditioning to achieve potentially better performance.

Our algorithm is general and can be applied to most of the kernel machines which are formulated in the convex optimization (1) in a RKHS H associated with given kernel ). We will instantiate the doubly stochastic gradients algorithms for a few commonly used kernel machines for different tasks and loss functions, e.g., regression, classification, quantile regression, novelty detection and estimating divergence functionals/likelihood ratio. Interestingly, the Gaussian process regression, which is a Bayesian model, can also be reformulated as the solution to particular convex optimizations in RKHS, and therefore, be approximated by the proposed algorithm.

Kernel Support Vector Machine (SVM). Hinge loss is used in kernel SVM where

with . We have

we have

Kernel Logistic Regression. Log loss is used in kernel logistic regression for binary classification where )) with . We have and the step 5 in Algorithm. 1. becomes

C is the number of categories, and , otherwise In such scenario, we denote )], and therefore, the corresponding ]. The update rule for in Algorithm. 1. is

Kernel Ridge Regression. Square loss is used in kernel ridge regression where . We

have ) and the step 5 in Algorithm. 1. becomes

Kernel Robust Regression. Huber’s loss is used for robust regression [22] where

We have

-insensitive loss function is used in kernel SVR where

. We have

becomes

Remark: Note that if we set -intensitive loss function will become absolute deviatin, i.e., . Therefore, we have the updates for kernel least absolute deviatin regression. Kernel Quantile Regression. The loss function for quantile regression is (1

. We have

and the gradient of ) is

The step 5 in Algorithm. 1. becomes

Kernel Density Ratio Estimation. Based on the variational form of Ali-Silvey divergence, ( , where is a convex function with r(1) = 0, [24] proposed a nonparametric estimator for the logarithm of the density ratio, log , which is the solution of following convex optimization,

where denotes the Fenchel-Legendre dual of ). In Kullback-Leibler (KL) divergence, the log(). Its Fenchel-Legendre dual is

where Bernoulli(0.5). Denote ) exp(, we have

and the the step 5 in Algorithm. 1. becomes

In particular, the and are not sampled in pair, they are sampled independently from P(x) and Q(x) respectively.

[24] proposed another convex optimization based on ) whose solution is a nonparametric estimator for the density ratio. [25] designed log ) for novelty detection. Similarly, the doubly stochastic gradients algorithm is also applicable to these loss functions. Gaussian Process Regression. The doubly stochastic gradients can be used for approximating the posterior of Gaussian process regression by reformulating the mean and variance of the predictive distribution as the solutions to the convex optimizations with particular loss functions. Let y = f(x) + where ) and )), given the dataset , the posterior distribution of the function at the test point can be derived as

where and is the identity matrix.

Obviously, the posterior mean of the Gaussian process for regression can be thought as the solution to optimization problem (1) with square loss and setting . Therefore, the update rule for approximating the posterior mean will be the same as kernel ridge regression.

To compute the predictive variance, we need to evaluate the . Following, we will introduce two different optimizations whose solutions can be used for evaluating the quantity.

where is the Hilbert-Schmidt norm of the operator. We can achieve the optimum by which is equivalent to Eq. 10.

2. Assume that the testing points, , are given beforehand, instead of approximating the op- erator A, we target on functions where (, [)] and . Estimating () can be accomplished by solving the optimization problem (1) with square loss and setting , leading to the same update rule as kernel ridge regression.

After we obtain these estimators, we can calculate the predictive variance on by either ) ) or )(). We conduct experiments to justify the novel formulations for approximating both the mean and variance of posterior of Gaussian processes for regression, and the doubly stochastic update rule in Section.(7).

Note that, to approximate the operator A, doubly stochastic gradient requires ) memory. Although we do not need to save the whole training dataset, which saves O(dt) memory cost, this is still computationally expensive. When the m testing data are given, we estimate m functions and each of them requires O(t) memory cost, the total cost will be O(tm) by the second algorithm.

5 Theoretical Guarantees

In this section, we will show that, both in expectation and with high probability, our algorithm can estimate the optimal function in the RKHS with rate O(1/t), and achieve a generalization bound of ). The analysis for our algorithm has a new twist compared to previous analysis of stochastic gradient descent algorithms, since the random feature approximation results in an estimator which is outside the RKHS. Besides the analysis for stochastic functional gradient descent, we need to use martingales and the corresponding concentration inequalities to prove that the sequence of estimators, , outside the RKHS converge to the optimal function, , in the RKHS. We make the following standard assumptions ahead for later references: A. There exists an optimal solution, denoted as , to the problem of our interest (1). B. Loss function ) : and its first-order derivative is L-Lipschitz continous in terms of the first argument.

Theorem 4 (Convergence in expectation) (1, 2) 2+ 2

+ (21)(1 + 1)+

φ)LMθκ + φ)M

Figure 2: stands the error due to random features, and stands for the error due to random data.

+ (1 + + 16(8+

(

Proof sketch: We focus on the convergence in expectation; the high probability bound can be established in a similar fashion. The main technical difficulty is that may not be in the RKHS H. The key of the proof is then to construct an intermediate function , such that the difference between and and the difference between and can be bounded. More specifically,

For the error term due to random features, is constructed such that is a martingale, and the stepsizes are chosen such that , which allows us to bound the martingale. In other words, the choices of the stepsizes keep close to the RKHS. For the error term due to random data, since , we can now apply the standard arguments for stochastic approximation in the RKHS. Due to the additional randomness, the recursion is slightly more complicated, + + where ], and and depends on the related parameters. Solving this recursion then leads to a bound for the second error term.

Proof By the Lipschitz continuity of ) and Jensen’s Inequality, we have

Again, can be decomposed as two terms and ), which can be bounded similarly as in Theorem 5 (see Corollary 12 in the appendix).

Remarks. The overall rate of convergence in expectation, which is O(1/t), is indeed optimal. Classical complexity theory (see, e.g. reference in [16]) shows that to obtain -accuracy solution, the number of iterations needed for the stochastic approximation is Ω(1) for strongly convex case and Ω(1) for general convex case. Different from the classical setting of stochastic approximation, our case imposes not one but two

Table 2: Comparison of Computation and Memory Requirements

sources of randomness/stochasticity in the gradient, which intuitively speaking, might require higher order number of iterations for general convex case. However, the variance of the random features only contributes additively to the constant in the final convergence rate. Therefore, our method is still able to achieve the same rate as in the classical setting. Notice that these bounds are achieved by adopting the classical stochastic gradient algorithm, and they may be further refined with more sophisticated techniques and analysis. For example, techniques for reducing variance of SGD proposed in [37], mini-batch and preconditioning [41, 42] can be used to reduce the constant factors in the bound significantly. Theorem 4 also reveals bounds in and sense as in Appendix B. The choices of stepsizes and the tuning parameters given in these bounds are only for sufficient conditions and simple analysis; other choices can also lead to bounds in the same order.

6 Computation, Memory and Statistics Trade-oﬀ

To investigate computation, memory and statistics trade-off, we will fix the desired error in the function estimation to , and work out the dependency of other quantities on . These other quantities include the preprocessing time, the number of samples and random features (or rank), the number of iterations of each algorithm, and the computational cost and memory requirement for learning and prediction. We assume that the number of samples, n, needed to achieve the prescribed error is of the order ), the same for all methods. Furthermore, we make no other regularity assumption about margin properties or the kernel matrix such as fast spectrum decay. Thus the required number of random feature (or ranks), r, will be of the order ) [4, 5, 8, 9].

We will pick a few representative algorithms for comparison, namely, (i) NORMA [13]: kernel methods trained with stochastic functional gradients; (ii) k-SDCA [12]: kernelized version of stochastic dual coordinate ascend; (iii) r-SDCA: first approximate the kernel function with random features, and then run stochastic dual coordinate ascend; (iv) n-SDCA: first approximate the kernel matrix using Nystr¨om’s method, and then run stochastic dual coordinate ascend; similarly we will combine Pegasos algorithm [26], stochastic block mirror descent (SBMD) [38], and random block coordinate descent (RBCD) [39] with random features and Nystr¨om’s method, and obtain (v) r-Pegasos, (vi) n-Pegasos, (vii) r-SBMD, (viii) n-SBMD, (ix) r-RBCD, and (x) n-RBCD, respectively. The comparisons are summarized below in Table. 2

From Table 2, one can see that our method, r-SDCA, r-Pegasos, r-SBMD and r-RBCD achieve the best dependency on the dimension, d, of the data up to a log factor. However, often one is interested in increasing the number of random features as more data points are observed to obtain a better generalization ability, e.g., in streaming setting. Then special procedures need to be designed for updating the r-SDCA, rPegasos, r-SBMD and r-RBCD solutions, which is not clear how to do easily and efficiently with theoretical

Table 3: Comparison of Computation and Memory Requirement Per Iteration. b denotes the block size in algorithms SBMD and RBCD.

Table 4: Datasets

guarantees. As a more refined comparison, our algorithm is also the cheapest in terms of per training iteration computation and memory requirement. We list the computational and memory requirements at a particular iteration t < n for these five algorithms to achieve error in Table 3.

7 Experiments

We show that our method compares favorably to other scalable kernel methods in medium scale datasets, and neural nets in large scale datasets. Below is a summary of the datasets used. A “yes” for the last column means that virtual examples (random cropping and mirror imaging of the original pictures) are generated for training. K-ridge stands for kernel ridge regression; GPR stands for Gaussian processes regression; K-SVM stands for kernel SVM; K-logistic stands for kernel logistic regression.

Experiment settings. We first justify the doubly stochastic algorithm for Gaussian processes regression on dataset (1), comparing with NORMA. The dataset is medium size, so that the closed-form for posterior is tractable. For the large-scale datasets (2) — (5), we compare with the first seven algorithms for solving kernel methods discussed in Table 2. For the algorithms based on low rank kernel matrix approximation and random features, i.e., pegasos and SDCA, we set the rank r or number of random features r to be 2. We use the same batch size for both our algorithms and the competitors. We adopted two stopping criteria for different purposes. We first stopped the algorithms when they pass through the entire dataset once (SC1). This stopping criterion is designed for justifying our motivation. By investigating the performances of these algorithms with different levels of random feature approximations but the same number of training samples, we could identify that the bottleneck of the performances of the vanilla methods with explicit feature will be their approximation ability. To further demonstrate the advantages of the proposed algorithm in computational cost, we also conduct experiments on datasets (3) – (5) running the competitors within the same time budget as the proposed algorithm (SC2). We do not count the preprocessing time of Nystr¨om’s method for n-Pegasos and n-SDCA, though it takes substantial amount of time. The algorithms are executed on the machine with AMD 16 2.4GHz Opteron CPUs and 200G memory. It should be noticed that this gives advantage to NORMA and k-SDCA which could save all the data in the memory. For fairness, we also record

Figure 3: Illustration of the neural nets structure in our experiments. The first several red layers are convolutions with max pooling layers. The following blue layers are fully connected layes. The green layer is the output layer which is multiclass logistic regression model.

as many random features as the memory allowed. For datasets (6) — (8), we compare with neural nets for images (“jointly-trained”). In order to directly compare the performance of nonlinear classifiers rather than feature learning abilities, we also use the convolution layers of a trained neural net to extract features, then apply our algorithm and a nonlinear neural net on top to learn classifiers (“fixed”). The structures of these neural nets in Figure 3. For datasets (9) and (10), we compare with the neural net described in [30] and use exactly the same input. In all the experiments, we select the batch size so that for each update, the computation resources can be utilized efficiently.

7.1 Kernel Ridge Regression

In this section, we compare our approach with alternative algorithms for kernel ridge regression on 2D synthetic dataset. The data are generated by

where 5]and 1). We use Gaussian RBF kernel with kernel bandwidth chosen to be 0.1 times the median of pairwise distances between data points (median trick). The regularization parameter is set to be 10. The batch size and feature block are set to be 2.

Figure 4: Experimental results for kernel ridge regression on synthetic dataset.

Figure 5: Experimental results for Gaussian Processes regression.

The results are shown in Figure 4. In Figure 4(1), we plot the optimal functions generating the data. We justify our proof of the convergence rate in Figure 4(2). The blue dotted line is a convergence rate of 1/t as a guide. ˆdenotes the average solution after t-iteration, i.e., ˆ). It could be seen that our algorithm indeed converges in the rate of O(1/t). In Figure 4 (3), we compare the first seven algorithms listed in the Table 2 for solving the kernel ridge regression.

The comparison on synthetic dataset demonstrates the advantages of our algorithm clearly. Our algorithm achieves comparable performance with NORMA, which uses full kernel, in similar time but less memory cost. The pegasos and SDCA using 2random or Nystr¨om features perform worse.

7.2 Gaussian Processes Regression

As we introduced in Section. (4), the mean and variance of posterior of Gaussian processes for regression problem can be formulated as solutions to some convex optimization problems. We conduct experiments on synthetic dataset for justification. Since the task is computing the posterior, we evaluate the performances by comparing the solutions to the posterior mean and variance, denoted as and , obtained by closed- form (9). We select 2data from the same model in previous section for training and 2data for testing, so that the closed-form of posterior is tractable. We use Gaussian RBF kernel with kernel bandwidth chosen by median trick. The noise level is set to be 0.1. The batch size is set to be 64 and feature block is set to be 512.

We compared the doubly stochastic algorithm with NORMA. The results are shown in Figure 5. Both the doubly stochastic algorithm and NORMA converge to the posterior, and our algorithm achieves comparable performance with NORMA in approximating both the mean and variance.

Figure 6: Comparison with other kernel SVM solvers on datasets (3) – (5) with two different stopping criteria.

7.3 Kernel Support Vector Machine

We evaluate our algorithm solving kernel SVM on three datasets (3)–(5) comparing with other several algorithms listed in Table 2 using stopping criteria SC1 and SC2.

Adult. We use Gaussian RBF kernel with kernel bandwidth obtained by median trick. The regularization parameter is set to be 1/(100n) where n is the number of training samples. We set the batch size to be 2and feature block to be 2. After going through the whole dataset one pass, the best error rate is achieved by NORMA and k-SDCA which is 15% while our algorithm achieves comparable result 15.3%. The performances are illustrated in Figure 6(1). Under the same time budget, all the algorithms perform similarly in Figure 6(4). The reason of flat region of r-pegasos, NORMA and the proposed method on this dataset is that Adult dataset is unbalanced. There are about 24% positive samples while 76% negative samples.

MNIST 8M 8 vs. 6. We first reduce the dimension to 50 by PCA and use Gaussian RBF kernel with kernel bandwidth 03 obtained by median trick. The regularization parameter is set to be 1/n where n is the number of training samples. We set the batch size to be 2and feature block to be 2. The results are shown in Figure 6(2) and (5) under SC1 and SC2 respectively. Under both these two stopping criteria, our algorithm achieves the best test error 0.26% using similar training time.

Forest. We use Gaussian RBF kernel with kernel bandwidth obtained by median trick. The regularization parameter is set to be 1/n where n is the number of training samples. We set the batch size to be 2and feature block to be 2. In Figure 6(3), we shows the performances of all algorithms using SC1. NORMA and k-SDCA achieve the best error rate, which is 10%, while our algorithm achieves around 15%, but still much better than the pegasos and SDCA with 2features. In the same time budget, the proposed algorithm performs better than all the alternatives except NORMA in Figure 6(6).

As seen from the performance of pegasos and SDCA on Adult and MNIST, using fewer features does not deteriorate the classification error. This might be because there are cluster structures in these two binary classification datasets. Thus, they prefer low rank approximation rather than full kernel. Different from

Figure 7: Comparison with Neural Networks on datasets (6) – (10).

these two datasets, in the forest dataset, algorithms with full kernel, i.e., NORMA and k-SDCA, achieve best performance. With more random features, our algorithm performs much better than pegasos and SDCA under both SC1 and SC2. Our algorithm is preferable for this scenario, i.e., huge dataset with sophisticated decision boundary. Although utilizing full kernel could achieve better performance, the computation and memory requirement for the kernel on huge dataset are costly. To learn the sophisticated boundary while still considering the computational and memory cost, we need to efficiently approximate the kernel in O( ) with O(n) random features at least. Our algorithm could handle so many random features efficiently in both computation and memory cost, while for pegasos and SDCA such operation is prohibitive.

7.4 Classification Comparisons to Convolution Neural Networks

We also compare our algorithm with the state-of-the-art neural network. In these experiments, the block size is set to be O(10). Compared to the number of samples, O(10), this block size is reasonable.

MNIST 8M. In this experiment, we compare to a variant of LeNet-5 [32], where all tanh units are replaced with rectified linear units. We also use more convolution filters and a larger fully connected layer. Specifically, the first two convolutions layers have 16 and 32 filters, respectively, and the fully connected layer contains 128 neurons. We use kernel logistic regression for the task. We extract features from the last max-pooling layer with dimension 1568, and use Gaussian RBF kernel with kernel bandwidth equaling to four times the median pairwise distance. The regularization parameter is set to be 0.0005.

The result is shown in Figure 7(1). As expected, the neural net with pre-learned features is faster to train than the jointly-trained one. However, our method is much faster compared to both methods. In addition, it achieves a lower error rate (0.5%) compared to the 0.6% error provided by the neural nets.

From molecular structure to molecular property

Figure 8: The computational procedure for predicting molecular property from molecular structure.

CIFAR 10. In this experiment, we compare to a neural net with two convolution layers (after contrast normalization and max-pooling layers) and two local layers that achieves 11% test erroron CIFAR 10 [28]. The features are extracted from the top max-pooling layer from a trained neural net with 2304 dimension. We use kernel logistic regression for this problem. The kernel bandwidth for Gaussian RBF kernel is again four times the median pairwise distance. The regularization parameter is set to be 0.0005. We also perform a PCA (without centering) to reduce the dimension to 256 before feeding to our method.

The result is shown in Figure 7(2). The test error for our method drops significantly faster in the earlier phase, then gradually converges to that achieved by the neural nets. Our method is able to produce the same performance within a much restricted time budget.

ImageNet. In this experiment, we compare our algorithm with the neural nets on the ImageNet 2012 dataset, which contains 1.3 million color images from 1000 classes. Each image is of size 256 256, and we randomly crop a 240 240 region with random horizontal flipping. The jointly-trained neural net is Alex-net [29]. The 9216 dimension features for our classifier and fixed neural net are from the last pooling layer of the jointly-trained neural net. The kernel bandwidth for Gaussian RBF kernel is again four times the median pairwise distance. The regularization parameter is set to be 0.0005.

Test error comparisons are shown in Figure 7(3). Our method achieves a test error of 44.5% by further max-voting of 10 transformations of the test set while the jointly-trained neural net arrives at 42% (without variations in color and illumination). At the same time, fixed neural net can only produce an error rate of 46% with max-voting. There may be some advantages to train the network jointly such that the layers work together to achieve a better performance. Although there is still a gap to the best performance by the jointly-trained neural net, our method comes very close with much faster convergence rate. Moreover, it achieves superior performance than the neural net with pre-learned features, both in accuracy and speed.

7.5 Regression Comparisons to Neural Networks

We test our algorithm for kernel ridge regression with neural network proposed in [30] on two large-scale real-world regression datasets, (9) and (10) in Table 4. To our best knowledge, this is the first comparison between kernel ridge regression and neural network on the dataset MolecularSpace.

QuantumMachine. In this experiment, we use the same binary representations converted based on random Coulomb matrices as in [30]. We first generate a set of randomly sorted coulomb matrices for each molecule. And then, we break each dimension of the Coulomb matrix apart into steps and convert them to the binary predicates. Predictions are made by taking average of all prediction made on various Coulomb matrices of the same molecule. The procedure is illustrated in Figure. 8. For this experiment, 40 sets of randomly permuted matrices are generated for each training example and 20 for each test example. We use Gaussian kernel with kernel bandwidth and the feature block is 2. The total dimension of random features is 2.

The results are shown in Figure 7(4). In QuantumMachine dataset, our method achieves Mean Absolute Error (MAE) of 2.97 kcal/mole, outperforming neural nets results, 3.51 kcal/mole. Note that this result is already close to the 1 kcal/mole required for chemical accuracy.

MolecularSpace. In this experiment, the task is to predict the power conversion efficiency (PCE) of the molecule. This dataset of 2.3 million molecular motifs is obtained from the Clean Energy Project Database. We use the same feature representation as for “QuantumMachine” dataset [30]. We set the kernel bandwidth of Gaussian RBF kernel to be 290 by median trick. The batch size is set to be 25000 and the feature block is 2. The total dimension of random features is 2.

The results are shown in Figure 7(5). It can be seen that our method is comparable with neural network on this 2.3 million dataset.

8 Discussion

Our work contributes towards making kernel methods scalable for large-scale datasets. Specifically, by introducing artificial randomness associated with kernels besides the random data samples, we propose doubly stochastic functional gradient for kernel machines which makes the kernel machines efficient in both computation and memory requirement. Our algorithm successfully reduces the memory requirement of kernel machines from O(dn) to O(n). Meanwhile, we also show that our algorithm achieves the optimal rate of convergence, O(1/t), for strongly convex stochastic optimization. We compare our algorithm on both classification and regression problems with the state-of-the-art neural networks as well as some other competing algorithms for kernel methods on several large-scale datasets. With our efficient algorithm, kernel methods could perform comparable to sophisticated-designed neural network empirically.

The theoretical analysis, which provides the rate of convergence independent to the dimension, is also highly non-trivial. It twists martingale techniques and the vanilla analysis for stochastic gradient descent and provides a new perspective for analyzing optimization in infinite-dimensional spaces, which could be of independent interest. It should be pointed out that although we applied the algorithm to many kernel machines even with non-smooth loss functions, our current proof relies on the Lipschitz smoothness of the loss function. Extending the guarantee to non-smooth loss function will be one interesting future work.

Another key property of our method is its simplicity and ease of implementation which makes it versatile and easy to be extened in various aspects. It is straightforward to replace the sampling strategy for random features with Fastfood [7] which enjoys the efficient computational cost, or Quasi-Monte Carlo sampling [43], data-dependent sampling [47] which enjoys faster convergence rate with fewer generated features. Meanwhile, by back-propogation trick, we could refine the random features by adapting their weights for better performance [36].

Acknowledgement

M.B. is supported in part by NSF grant CCF-1101283, AFOSR grant FA9550-09-1-0538, a Microsoft Faculty Fellowship, and a Raytheon Faculty Fellowship. L.S. is supported in part by NSF IIS-1116886, NSF/NIH BIGDATA 1R01GM108341, NSF CAREER IIS-1350983, and a Raytheon Faculty Fellowship.

References

[1] A. J. Smola and B. Sch¨olkopf. Sparse greedy matrix approximation for machine learning. In ICML, pages 911–918, San Francisco, 2000. Morgan Kaufmann Publishers.

[2] C. K. I. Williams and M. Seeger. Using the Nystrom method to speed up kernel machines. In T. G. Dietterich, S. Becker, and Z. Ghahramani, editors, NIPS, 2000.

[3] S. Fine and K. Scheinberg. Efficient SVM training using low-rank kernel representations. JMLR, 2:243–264, 2001.

[4] P. Drineas and M. Mahoney. On the nystr om method for approximating a gram matrix for improved kernel-based learning. JMLR, 6:2153–2175, 2005.

[5] Corinna Cortes, Mehryar Mohri, and Ameet Talwalkar. On the impact of kernel approximation on learning accuracy. In AISTATS, pages 113–120, 2010.

[6] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In J.C. Platt, D. Koller, Y. Singer, and S. Roweis, editors, NIPS. MIT Press, Cambridge, MA, 2008.

[7] Q.V. Le, T. Sarlos, and A. J. Smola. Fastfood — computing hilbert space expansions in loglinear time. In ICML, 2013.

[8] Ali Rahimi and Benjamin Recht. Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning. In NIPS, 2009.

[9] David Lopez-Paz, Suvrit Sra, A. J. Smola, Zoubin Ghahramani, and Bernhard Schlkopf. Randomized nonlinear component analysis. In ICML, 2014.

[10] John C. Platt. Sequential minimal optimization: A fast algorithm for training support vector machines. Technical Report MSR-TR-98-14, Microsoft Research, 1998.

[11] T. Joachims. Making large-scale SVM learning practical. In B. Sch¨olkopf, C. J. C. Burges, and A. J. Smola, editors, Advances in Kernel Methods — Support Vector Learning, pages 169–184, Cambridge, MA, 1999. MIT Press.

[12] Shai Shalev-Shwartz and Tong Zhang. Stochastic dual coordinate ascent methods for regularized loss. JMLR, 14(1):567–599, 2013.

[13] J. Kivinen, A. J. Smola, and R. C. Williamson. Online learning with kernels. IEEE Transactions on Signal Processing, 52(8), Aug 2004.

[14] S. S. Keerthi and D. DeCoste. A modified finite Newton method for fast solution of large scale linear SVMs. J. Mach. Learn. Res., 6:341–361, 2005.

[15] N. Ratliff and J. Bagnell. Kernel conjugate gradient for fast kernel machines. In IJCAI, volume 20, January 2007.

[16] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro. Robust stochastic approximation approach to stochastic programming. SIAM J. on Optimization, 19(4):1574–1609, January 2009.

[17] A. Devinatz. Integral representation of pd functions. Trans. AMS, 74(1):56–77, 1953.

[18] M. Hein and O. Bousquet. Kernels, associated structures, and generalizations. Technical Report 127, Max Planck Institute for Biological Cybernetics, 2004.

[19] H. Wendland. Scattered Data Approximation. Cambridge University Press, Cambridge, UK, 2005.

[20] Bernhard Sch¨olkopf and A. J. Smola. Learning with Kernels. MIT Press, Cambridge, MA, 2002.

[21] N. Pham and R. Pagh. Fast and scalable polynomial kernels via explicit feature maps. In KDD. ACM, 2013.

[22] K.-R. M¨uller, A. J. Smola, G. R¨atsch, B. Sch¨olkopf, J. Kohlmorgen, and V. Vapnik. Predicting time series with support vector machines. In W. Gerstner, A. Germond, M. Hasler, and J.-D. Nicoud, editors, , volume 1327 of Lecture Notes in Comput. Sci., pages 999–1004, Berlin, 1997. Springer-Verlag.

[23] B. Sch¨olkopf, J. Platt, J. Shawe-Taylor, A. J. Smola, and R. C. Williamson. Estimating the support of a high-dimensional distribution. Neural Computation, 13(7):1443–1471, 2001.

[24] X.L. Nguyen, M. Wainwright, and M. Jordan. Estimating divergence functionals and the likelihood ratio by penalized convex risk minimization. In Advances in Neural Information Processing Systems 20, pages 1089–1096. MIT Press, Cambridge, MA, 2008.

[25] Alex J Smola, Le Song, and Choon H Teo. Relative novelty detection. In International Conference on Artificial Intelligence and Statistics, pages 536–543, 2009.

[26] Shai Shalev-Shwartz, Yoram Singer, and Nathan Srebro. Pegasos: Primal estimated sub-gradient solver for SVM. In ICML, 2007.

[27] G. Loosli, S. Canu, and L. Bottou. Training invariant support vector machines with selective sampling. In L. Bottou, O. Chapelle, D. DeCoste, and J. Weston, editors, Large Scale Kernel Machines, pages 301–320. MIT Press, 2007.

[28] A. Krizhevsky. Learning multiple layers of features from tiny images. Technical report, University of Toronto, 2009.

[29] A. Krizhevsky, I. Sutskever, and G. Hinton. Imagenet classification with deep convolutional neural networks. In NIPS, 2012.

[30] Gr´egoire Montavon, Katja Hansen, Siamac Fazli, Matthias Rupp, Franziska Biegler, Andreas Ziehe, Alexandre Tkatchenko, Anatole von Lilienfeld, and Klaus-Robert M¨uller. Learning invariant representations of molecules for atomization energy prediction. In NIPS, pages 449–457, 2012.

[31] Alexander Rakhlin, Ohad Shamir, and Karthik Sridharan. Making gradient descent optimal for strongly convex stochastic optimization. In ICML, pages 449–456, 2012.

[32] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied to document recog- nition. Proceedings of the IEEE, 86(11):2278–2324, November 1998.

[33] Purushottam Kar and Harish Karnick. Random feature maps for dot product kernels. In Neil D. Lawrence and Mark A. Girolami, editors, AISTATS-12, volume 22, pages 583–591, 2012.

[34] Andrea Vedaldi and Andrew Zisserman. Efficient additive kernels via explicit feature maps. IEEE Trans. Pattern Anal. Mach. Intell., 34(3):480–492, 2012.

[35] Jiyan Yang, Vikas Sindhwani, Quanfu Fan, Haim Avron, and Michael W. Mahoney. Random laplace feature maps for semigroup kernels on histograms. In CVPR, 2014.

[36] Zichao Yang, Marcin Moczulski, Misha Denil, Nando de Freitas, Alexander J. Smola, Le Song, and Ziyu Wang. Deep fried convnets. CoRR, abs/1412.7149, 2014d. URL http://arxiv.org/abs/1412.7149.

[37] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduc- tion. In NIPS, pages 315–323, 2013.

[38] Cong D. Dang and Guanghui Lan. Stochastic block mirror descent methods for nonsmooth and stochas- tic optimization. Technical report, University of Florida, 2013.

[39] Yurii Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization, 22(2):341–362, 2012.

[40] Andrew Cotter, Shai Shalev-Shwartz, and Nati Srebro. Learning optimally sparse support vector ma- chines. In Proceedings of the 30th International Conference on Machine Learning, ICML 2013, Atlanta, GA, USA, 16-21 June 2013, pages 266–274, 2013.

[41] A. Agarwal, S. Kakade, N. Karampatziakis, L. Song, and G. Valiant. Least squares revisited: Scalable approaches for multi-class prediction. In International Conference on Machine Learning (ICML), 2014.

[42] Tianbao Yang, Rong Jin, and Shenghuo Zhu. On data preconditioning for regularized loss minimization. CoRR, 2014.

[43] Jiyan Yang, Vikas Sindhwani, Haim Avron, and Michael W. Mahoney. Quasi-monte carlo feature maps for shift-invariant kernels. In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014, pages 485–493, 2014.

[44] Gaurav Pandey and Ambedkar Dukkipati. Learning by stretching deep networks. In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014, pages 485–493, 2014.

[45] Youngmin Cho and Lawrence K. Saul. Kernel methods for deep learning. In Y. Bengio, D. Schuurmans, J.D. Lafferty, C.K.I. Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems 22, pages 342–350, 2009.

[46] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA, 2006.

[47] Francis R. Bach. On the equivalence between quadrature rules and random features. CoRR, abs/1502.06800, 2015. URL http://arxiv.org/abs/1502.06800.

Appendix A Convergence Rate

We first provide specific bounds and detailed proofs for the two error terms appeared in Theorem 4 and Theorem 5.

A.1 Error due to random features

Lemma 7 We have

Then by Azuma’s Inequality, for any 0,

which is equivalent as

Moreover,

Lemma 8 Suppose (1 ) and (1, 2) . Then we have

Proof (1) follows by induction on is trivially true. We have

When (1, 2), + 1 for any 1, so . When , if 1, then ; if 1, then [1, 2),

When and 2

A.2 Error due to random data

Lemma 9 Assume ) is L-Lipschitz continous in terms of . Let be the optimal solution to our target problem. Then

Proof For the sake of simple notations, let us first denote the following three different gradient terms, which are

Note that by our previous definition, we have 1.

Hence, we have

Proof for (i): Let us denote . We first show that are bounded. Specifically, we have for 1,

We prove these results separately in Lemma 10 below. Let us denote ], given the above bounds, we arrive at the following recursion,

When with such that (1, 2) , from Lemma 8, we have . Consequently, and . Applying these bounds leads to the refined recursion as follows

that can be further written as

where and (1 + . Invoking Lemma 14 with 1, we obtain

Proof for (ii): Cumulating equations (12) with i = 1, . . . t, we end up with the following inequality

We first show that

Again, the proofs of these results are given separately in Lemma 10. Applying the above bounds leads to the refined recursion as follows,

with probability 1. When with such that , with similar reasons in Lemma 8, we have and also we have (1 ) )(1 . Therefore, we can rewrite the above recursion as

where (1+(1+. Invoking Lemma 15, we obtain

with the specified .

Lemma 10 In this lemma, we prove the inequalities (1)–(5) in Lemma 9.

Proof Given the definitions of in Lemma 9, we have

where the first and third inequalities are due to Cauchy–Schwarz Inequality and the second inequality is due to L-Lipschitz continuity of ) in the first parameter, and the last step is due to Lemma 7 and the definition of .

is martingale difference sequence since max, with .

• V ar(d4κM t.

where ) and , we immediately obtain the above inequality as desired.

Applying these lemmas immediately gives us Theorem 4 and Theorem 5, which implies pointwise distance between the solution ) and ). Now we prove similar bounds in the sense of and distance.

B L∞ distance, L2 distance, and generalization bound

Theorem 4 also implies a bound in sense, namely, + 2

This is because , where always exists since X is closed and bounded. Note that the result for average solution can be improved without log factor using more sophisticated analysis (see also reference in [31]).

With the choices of in Lemma 9, we have

Proof (i) follows directly from Theorem 4. (ii) can be proved as follows. First, we have

From Lemma 9, with probability at least 1 , we have

From Lemma 7, for any , we have

Since , the above inequality can be writen as

which leads to

By Fubini’s theorem and Markov’s inequality, we have

From the analysis in Lemma 7, we also have that . Therefore, with probability at least 1 over (), we have

Summing up equation (17) and (16), we have

From the bound on distance, we can immediately get the generalization bound.

Then the theorem follows from Corollary 12.

C Suboptimality

For comprehensive purposes, we also provide the O(1/t) bound for suboptimality.

Corollary 13 If we set with , then the average solution ˆsatisfies

where + 2, with defined as in Lemma 9.

Proof From the anallysis in Lemma 9,we have

Invoking strongly convexity of R(f), we have )+ . Taking expectaion on both size and use the bounds in last lemma, we have

which can be further bounded by

By convexity, we have (ˆ. The corollary then follows from the fact that [ ˆ[ˆ] and [ˆ(ˆ)].

C.1 Technical lemma for recursion bounds

Lemma 14 Suppose the sequence satisfies Γ0, and 1

where η > 1, β, β> 0

Proof The proof follows by induction. When t = 1, it always holds true by the definition of R. Assume the conclusion holds true for t with 1, i.e., Γ, then we have

where the last step follows from the defintion of .

Lemma 15 Suppose the sequence satisfies

Proof The proof follows by induction. When 1, therefore,

Since+2+ 2+ , we have (2+2+ + 2. Hence,

D Doubly Stochastic Gradient Algorithm for Posterior Variance Operator in Gaussian Process Regression

As we show in Section 4, the estimation of the variance of the predictive distribution of Gaussian process for regression problem could be recast as estimating the operator A defined in (10). We first demonstrate that the operator A is the solution to the following optimization problem

where is the Hilbert-Schmidt norm of the operator. The gradient of R(A) with respect to A is

Set + , exactly the same as (10).

To derive the doubly stochastic gradient update for A, we start with stochastic functional gradient of R(A). Given ), the stochastic functional gradient of R(A) is

With such update rule, we could show that ) by induction. Let ). Assume at t-th iteration, ), and notice that

Recall