Scalable Kernel Methods via Doubly Stochastic Gradients

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 $O(n^{2})$ space, and computing it takes $O(n^{2}d)$ 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.

Random feature approximation is another popular approach for scaling up kernel methods . 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 $O(1/\sqrt{r}+1/\sqrt{n})$ , 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 $r^{\prime}$. 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., ). 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., ). 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

the property of stochastic gradient descent algorithm that as long as the stochastic gradient is unbiased, the convergence of the algorithm is guaranteed ; and

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 , and achieve a generalization bound of $O(1/\sqrt{t})$. 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.

Duality between Kernels and Random Processes

Doubly Stochastic Functional Gradients

For instance, applying the above definition, we have $\nabla f(x)=\nabla\left\langle f,k(x,\cdot)\right\rangle_{\mathcal{H}}=k(x,\cdot)$, and $\nabla\left\|f\right\|_{\mathcal{H}}^{2}=\nabla\left\langle f,f\right\rangle_{\mathcal{H}}=2f$.

Note that the stochastic functional gradient $\xi(\cdot)$ is in RKHS $\mathcal{H}$ but $\zeta(\cdot)$ may be outside $\mathcal{H}$, since $\phi_{\omega}(\cdot)$ may be outside the RKHS. For instance, for the Gaussian RBF kernel, the random feature $\phi_{\omega}(x)=\sqrt{2}\cos(\omega^{\top}x+b)$ 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.

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 . 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.

where $a_{t}^{i}=-\gamma_{i}\prod_{j=i+1}^{t}(1-\gamma_{j}\nu)$ are deterministic values depending on the step sizes $\gamma_{j}(i\leqslant j\leqslant t)$ and regularization parameter $\nu$. 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 $\mathcal{H}$ associated with given kernel $k(x,x^{\prime})$. 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 $l(u,y)=\max\{0,1-uy\}$ with $y\in\{-1,1\}$. We have $l^{\prime}(u,y)=\begin{cases}0&\,\text{if }yu\geqslant 1\\ -y&\,\text{if }yu<1\\ \end{cases}$ and the step 5 in Algorithm. 1. becomes

Kernel Logistic Regression. Log loss is used in kernel logistic regression for binary classification where $l(u,y)=\log(1+\exp(-yu))$ with $y\in\{-1,1\}$. We have $l^{\prime}(u,y)=-\frac{y\exp(-yu)}{1+\exp(-yu)}$ and the step 5 in Algorithm. 1. becomes

Kernel Ridge Regression. Square loss is used in kernel ridge regression where $l(u,y)=\frac{1}{2}(u-y)^{2}$. We have $l^{\prime}(u,y)=(u-y)$ and the step 5 in Algorithm. 1. becomes

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

Remark: Note that if we set $\epsilon=0$, the $\epsilon$-intensitive loss function will become absolute deviatin, i.e., $l(u,y)=|u-y|$. Therefore, we have the updates for kernel least absolute deviatin regression.

Kernel Quantile Regression. The loss function for quantile regression is $l(u,y)=\max\{\tau(y-u),(1-\tau)(u-y)\}$. We have $l^{\prime}(u,y)=\begin{cases}1-\tau&\,\text{if }u\geqslant y\\ -\tau&\,\text{if }u<y\\ \end{cases}$ and the step 5 in Algorithm. 1. becomes

Kernel Novelty Detection. The loss function $l(u,\tau)=\max\{0,\tau-u\}$ is proposed for novelty detection. Since $\tau$ is also a variable which needs to be optimized, the optimization problem is formulated as

where $r^{*}$ denotes the Fenchel-Legendre dual of $r$, $r(\tau):=\sup_{\chi}\chi\tau-r^{*}(\chi)$. In Kullback-Leibler (KL) divergence, the $r_{KL}(\tau)=-\log(\tau)$. Its Fenchel-Legendre dual is

where $z\sim\text{Bernoulli}(0.5)$. Denote $l(u_{x},u_{y},z)=\delta_{1}(z)\exp(u_{y})-\delta_{0}(z)u_{x}$, we have

and the the step 5 in Algorithm. 1. becomes

proposed another convex optimization based on $r_{KL}(\tau)$ whose solution is a nonparametric estimator for the density ratio. designed $r_{nv}(\tau)=\max(0,\rho-\log\tau)$ 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)+\epsilon$ where $\epsilon\sim\mathcal{N}(0,\sigma^{2})$ and $f(x)\sim\mathcal{G}\mathcal{P}(0,k(x,x^{\prime}))$, given the dataset $\{x_{i},y_{i}\}_{i=1}^{n}$, the posterior distribution of the function at the test point $x_{*}$ can be derived as

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 $\nu=2\sigma^{2}$. 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 ${k^{*}}^{\top}\left(K+\sigma^{2}I\right)^{-1}k^{*}$. Following, we will introduce two different optimizations whose solutions can be used for evaluating the quantity.

Denote $\phi=[k(x_{1},\cdot),\ldots,k(x_{n},\cdot)]$, then

where the second equation based on identity $\left(\phi\phi^{\top}+\sigma^{2}I\right)\phi=\phi\left(\phi^{\top}\phi+\sigma^{2}I\right)$. Therefore, we just need to estimate the operator:

We can express $\mathcal{A}$ as the solution to the following convex optimization problem

where $\|\cdot\|_{HS}$ is the Hilbert-Schmidt norm of the operator. We can achieve the optimum by $\nabla R=0$, which is equivalent to Eq. 10.

Based on this optimization, we approximate the $\mathcal{A}_{t}$ using $\sum_{i\leqslant j,i=1}^{t}\theta_{ij}\phi_{\omega_{i}}(\cdot)\otimes\phi_{\omega_{j}}(\cdot)$ by doubly stochastic functional gradients. The update rule for $\theta$ is

Please refer to Appendix D for the details of the derivation.

Assume that the testing points, $\{x_{i}^{*}\}_{i=1}^{m}$, are given beforehand, instead of approximating the operator $\mathcal{A}$, we target on functions $F^{*}=[f^{*}_{1},\ldots,f^{*}_{m}]^{\top}$ where $f^{*}_{i}(\cdot)=k(\cdot)^{\top}\left(K+\sigma^{2}I\right)^{-1}k_{i}^{*}$, $k(\cdot)=[k(x_{1},\cdot),\ldots,k(x_{2},\cdot)]$ and $k_{i}^{*}=[k(x_{i}^{*},x_{1}),\ldots,k(x_{i}^{*},x_{n})]^{\top}$. Estimating $f^{*}_{i}(\cdot)$ can be accomplished by solving the optimization problem (1) with square loss and setting $y_{j}=k(x_{i}^{*},x_{j}),\forall j=1,\ldots,n$, $\nu=2\sigma^{2}$, leading to the same update rule as kernel ridge regression.

After we obtain these estimators, we can calculate the predictive variance on $x_{i}^{*}$ by either $k(x_{i}^{*},x_{i}^{*})-\mathcal{A}(x_{i}^{*},x_{i}^{*})$ or $k(x_{i}^{*},x_{i}^{*})-f_{i}^{*}(x_{i}^{*})$. 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 $\mathcal{A}$, doubly stochastic gradient requires $O(t^{2})$ 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.

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 $O(1/\sqrt{t})$. 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, $f_{t+1}$, outside the RKHS converge to the optimal function, $f_{\ast}$, in the RKHS. We make the following standard assumptions ahead for later references:

There exists an optimal solution, denoted as $f_{*}$, to the problem of our interest (1).

There exists $\kappa>0$ and $\phi>0$, such that $k(x,x^{\prime})\leqslant\kappa,\,|\phi_{\omega}(x)\phi_{\omega}(x^{\prime})|\leqslant\phi,\forall x,x^{\prime}\in\mathcal{X},\omega\in\Omega.$ For example, when $k(\cdot,\cdot)$ is the Gaussian RBF kernel, we have $\kappa=1$, $\phi=2$.

We now present our main theorems as below. Due to the space restrictions, we will only provide a short sketch of proofs here. The full proofs for the these theorems are given in the Appendix A-C.

where $Q_{1}=\max\left\{\left\|f_{\ast}\right\|_{\mathcal{H}},(Q_{0}+\sqrt{Q_{0}^{2}+(2\theta\nu-1)(1+\theta\nu)^{2}\theta^{2}\kappa M^{2}})/(2\nu\theta-1)\right\}$, with $Q_{0}=2\sqrt{2}\kappa^{1/2}(\kappa+\phi)LM\theta^{2}$, and $C^{2}=4(\kappa+\phi)^{2}M^{2}\theta^{2}$.

where $C$ is as above and $Q_{2}=\max\left\{\left\|f_{\ast}\right\|_{\mathcal{H}},Q_{0}+\sqrt{Q_{0}^{2}+\kappa M^{2}(1+\theta\nu)^{2}(\theta^{2}+16\theta/\nu)}\right\}$, with $Q_{0}=4\sqrt{2}\kappa^{1/2}M\theta(8+(\kappa+\phi)\theta L)$.

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 $f_{t+1}$ may not be in the RKHS $\mathcal{H}$. The key of the proof is then to construct an intermediate function $h_{t+1}$, such that the difference between $f_{t+1}$ and $h_{t+1}$ and the difference between $h_{t+1}$ and $f_{*}$ can be bounded. More specifically,

Proof By the Lipschitz continuity of $l(\cdot,y)$ and Jensen’s Inequality, we have

Again, $\|f_{t+1}-f_{*}\|_{2}$ can be decomposed as two terms $O\left(\|f_{t+1}-h_{t+1}\|_{2}^{2}\right)$ and $O(\left\|h_{t+1}-f_{\ast}\right\|_{\mathcal{H}}^{2})$, 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 ) shows that to obtain $\epsilon$-accuracy solution, the number of iterations needed for the stochastic approximation is $\Omega(1/\epsilon)$ for strongly convex case and $\Omega(1/\epsilon^{2})$ for general convex case. Different from the classical setting of stochastic approximation, our case imposes not one but two 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 , mini-batch and preconditioning can be used to reduce the constant factors in the bound significantly. Theorem 4 also reveals bounds in $L_{\infty}$ and $L_{2}$ sense as in Appendix B. The choices of stepsizes $\gamma_{t}$ 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.

Computation, Memory and Statistics Trade-off

To investigate computation, memory and statistics trade-off, we will fix the desired $L_{2}$ error in the function estimation to $\epsilon$, i.e., $\|f-f_{*}\|_{2}^{2}\leqslant\epsilon$, and work out the dependency of other quantities on $\epsilon$. 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 $\epsilon$ is of the order $O(1/\epsilon)$, 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 $O(n)=O(1/\epsilon)$ .

We will pick a few representative algorithms for comparison, namely, (i) NORMA : kernel methods trained with stochastic functional gradients; (ii) k-SDCA : 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öm’s method, and then run stochastic dual coordinate ascend; similarly we will combine Pegasos algorithm , stochastic block mirror descent (SBMD) , and random block coordinate descent (RBCD) with random features and Nyström’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. 2We only considered general kernel algorithms in this section. For some specific loss functions, e.g., hinge-loss, there are algorithms proposed to achieve better memory saving with extra training cost, such as support vector reduction technique .

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, r-Pegasos, r-SBMD and r-RBCD solutions, which is not clear how to do easily and efficiently with theoretical 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 $\epsilon$ error in Table 3.

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.

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^{8}$. 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öm’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 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 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.

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 $x\in^{2}$ and $e\sim\mathcal{N}(0,1)$. We use Gaussian RBF kernel with kernel bandwidth $\sigma$ chosen to be $0.1$ times the median of pairwise distances between data points (median trick). The regularization parameter $\nu$ is set to be $10^{-6}$. The batch size and feature block are set to be $2^{10}$.

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. $\hat{f}_{t}$ denotes the average solution after $t$-iteration, i.e., $\hat{f}_{t}(x)=\frac{1}{t}\sum_{i=1}^{t}f_{i}(x)$. 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 $2^{8}$ random or Nyström features perform worse.

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 $f_{gp}$ and $\sigma_{gp}^{2}$, obtained by closed-form (9). We select $2^{11}$ data from the same model in previous section for training and $2^{10}$ data for testing, so that the closed-form of posterior is tractable. We use Gaussian RBF kernel with kernel bandwidth $\sigma$ chosen by median trick. The noise level $\sigma^{2}$ 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.

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.

We use Gaussian RBF kernel with kernel bandwidth obtained by median trick. The regularization parameter $\nu$ is set to be $1/(100n)$ where $n$ is the number of training samples. We set the batch size to be $2^{6}$ and feature block to be $2^{5}$. 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 $\sigma=9.03$ obtained by median trick. The regularization parameter $\nu$ is set to be $1/n$ where $n$ is the number of training samples. We set the batch size to be $2^{10}$ and feature block to be $2^{8}$. 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 $\nu$ is set to be $1/n$ where $n$ is the number of training samples. We set the batch size to be $2^{10}$ and feature block to be $2^{8}$. 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 $2^{8}$ features. 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 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(\frac{1}{\epsilon})$ 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.

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^{4})$. Compared to the number of samples, $O(10^{8})$, this block size is reasonable.

In this experiment, we compare to a variant of LeNet-5 , 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 $\sigma$ equaling to four times the median pairwise distance. The regularization parameter $\nu$ 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.

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 errorThe specification is at https://code.google.com/p/cuda-convnet/ on CIFAR 10 . 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 $\sigma$ for Gaussian RBF kernel is again four times the median pairwise distance. The regularization parameter $\nu$ 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 $\times$ 256, and we randomly crop a 240 $\times$ 240 region with random horizontal flipping. The jointly-trained neural net is Alex-net . 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 $\sigma$ for Gaussian RBF kernel is again four times the median pairwise distance. The regularization parameter $\nu$ 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.

5 Regression Comparisons to Neural Networks

We test our algorithm for kernel ridge regression with neural network proposed in 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.

In this experiment, we use the same binary representations converted based on random Coulomb matrices as in . 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 $\sigma=60$ obtained by median trick. The batch size is set to be $50000$ and the feature block is $2^{11}$. The total dimension of random features is $2^{20}$.

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 . 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^{11}$. The total dimension of random features is $2^{20}$.

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.

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 which enjoys the efficient computational cost, or Quasi-Monte Carlo sampling , data-dependent sampling 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 .

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

Appendix A Convergence Rate

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

For any $x\in\mathcal{X}$, with probability at least $1-\delta$ over $(\mathcal{D}^{t},\bm{\omega}^{t})$,

Proof Let $V_{i}(x)=V_{i}(x;\mathcal{D}^{i},\bm{\omega}^{i}):=a_{t}^{i}\left(\zeta_{i}(x)-\xi_{i}(x)\right)$. Since $V_{i}(x)$ is a function of $(\mathcal{D}^{i},\bm{\omega}^{i})$ and

we have that $\left\{V_{i}(x)\right\}$ is a martingal difference sequence. Further note that

Then by Azuma’s Inequality, for any $\epsilon>0$,

Since $f_{t+1}(x)-h_{t+1}(x)=\sum_{i=1}^{t}V_{i}(x)$, we immediately obtain the two parts of the lemma.

$|a^{i}_{t}|\leqslant\frac{\theta}{t}$. Consequently, $\sum_{i=1}^{t}(a^{i}_{t})^{2}\leqslant\frac{\theta^{2}}{t}.$

Proof $(1)$ follows by induction on $i$. $|a_{t}^{t}|\leqslant\frac{\theta}{t}$ is trivially true. We have

A.2 Error due to random data

Particularly, if $\theta\nu=1$, we have $Q_{1}\leqslant\max\left\{\left\|f_{\ast}\right\|_{\mathcal{H}},4\sqrt{2}((\kappa+\phi)L+\nu)\cdot\frac{\kappa^{1/2}M}{\nu^{2}}\right\}$.

Particularly, if $\theta\nu=1$, we have $Q_{2}\leqslant\max\left\{\left\|f_{\ast}\right\|_{\mathcal{H}},8\sqrt{2}((\kappa+\phi)L+9\nu)\cdot\frac{\kappa^{1/2}M}{\nu^{2}}\right\}$.

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 $h_{t+1}=h_{t}-\gamma_{t}g_{t},\forall t\geqslant 1$.

Denote $A_{t}=\left\|h_{t}-f_{\ast}\right\|^{2}_{\mathcal{H}}$. Then we have

Because of the strongly convexity of (1) and optimality condition, we have

Proof for $(i)$: Let us denote $\mathcal{M}_{t}=\left\|g_{t}\right\|^{2}_{\mathcal{H}}$, $\mathcal{N}_{t}=\langle h_{t}-f_{\ast},\bar{g}_{t}-\hat{g}_{t}\rangle_{\mathcal{H}}$, $\mathcal{R}_{t}=\langle h_{t}-f_{\ast},\hat{g}_{t}-g_{t}\rangle_{\mathcal{H}}$. We first show that $\mathcal{M}_{t},\mathcal{N}_{t},\mathcal{R}_{t}$ are bounded. Specifically, we have for $t\geqslant 1$,

$\mathcal{M}_{t}\leqslant\kappa M^{2}(1+\nu c_{t})^{2}$, where $c_{t}=\sqrt{\sum_{i,j=1}^{t-1}|a_{t-1}^{i}|\cdot|a_{t-1}^{j}|}$ for $t\geqslant 2$ and $c_{1}=0$;

where $\beta_{1}=4\sqrt{2}\kappa^{1/2}LM(k+\phi)\theta^{2}$ and $\beta_{2}=\kappa M^{2}(1+\nu\theta)^{2}\theta^{2}$. Invoking Lemma 14 with $\eta=2\theta\nu>1$, we obtain

where $Q_{1}=\max\left\{\left\|f_{\ast}\right\|_{\mathcal{H}},\frac{Q_{0}+\sqrt{Q_{0}^{2}+(2\theta\nu-1)(1+\theta\nu)^{2}\theta^{2}\kappa M^{2}}}{2\nu\theta-1}\right\}$, and $Q_{0}=2\sqrt{2}\kappa^{1/2}(\kappa+\phi)LM\theta^{2}.$

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

Let us denote $b_{t}^{i}=\gamma_{i}\prod_{j={i+1}}^{t}(1-2\nu\gamma_{j}),1\leqslant i\leqslant t$, the above inequality is equivalent as

for any $0<\delta<1/e$ and $t\geqslant 4$, with probability $1-\delta$ over $(\mathcal{D}^{t},\bm{\omega}^{t})$,

where $C_{0}=\frac{4\max_{1\leqslant i\leqslant t}\mathcal{M}_{i}}{\nu}$.

for any $\delta>0$, with probability $1-\delta$ over $(\mathcal{D}^{t},\bm{\omega}^{t})$,

where $\hat{B}^{2}_{2,i}=2M^{2}(\kappa+\phi)^{2}\ln\left(\frac{2t}{\delta}\right)\sum_{j=1}^{i-1}|a^{j}_{i-1}|^{2}$.

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

where $\beta_{1}=\kappa M^{2}(1+\nu\theta)^{2}\theta^{2}$, $\beta_{2}=2\sqrt{2}\kappa^{1/2}LM(\kappa+\phi)\theta^{2}$, $\beta_{3}=16\kappa^{1/2}M\theta$, $\beta_{4}=16\kappa M^{2}(1+\theta\nu)^{2}\theta/\nu$. Invoking Lemma 15, we obtain

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

Proof Given the definitions of $\mathcal{M}_{t},\mathcal{N}_{t},\mathcal{R}_{t}$ in Lemma 9, we have

$\mathcal{M}_{t}\leqslant\kappa M^{2}(1+\nu\sqrt{\sum_{i,j=1}^{t-1}|a_{t-1}^{i}|\cdot|a_{t-1}^{j}|})^{2}$; This is because

where the first and third inequalities are due to Cauchy–Schwarz Inequality and the second inequality is due to $L$-Lipschitz continuity of $l^{\prime}(\cdot,\cdot)$ in the first parameter, and the last step is due to Lemma 7 and the definition of $A_{t}$.

for any $0<\delta<1/e$ and $t\geqslant 4$, with probability at least $1-\delta$ over $(\mathcal{D}^{t},\bm{\omega}^{t})$,

where $C_{0}=\frac{4\max_{1\leqslant i\leqslant t}\mathcal{M}_{i}}{\nu}$. This result follows directly from Lemma 3 in . Let us define $d_{i}=d_{i}(\mathcal{D}^{i},\bm{\omega}^{i}):=b_{t}^{i}\mathcal{N}_{i}=b_{t}^{i}\langle h_{i}-f_{\ast},\bar{g}_{i}-\hat{g}_{i}\rangle_{\mathcal{H}},1\leqslant i\leqslant t$, we have

$|d_{i}|\leqslant\max_{i}|b_{t}^{i}|\cdot C_{0}$, with $C_{0}=\frac{4\max_{1\leqslant i\leqslant t}\mathcal{M}_{i}}{\nu}$, $\forall 1\leqslant i\leqslant t$.

$Var(d_{i}|\mathcal{D}^{i-1},\bm{\omega}^{i-1})\leqslant 4\kappa M^{2}|b_{t}^{i}|^{2}A_{i},\forall 1\leqslant i\leqslant t$.

Plugging in these specific bounds in Lemma 3 in [Alexander et.al., 2012], which is,

where $\sigma_{t}^{2}=\sum_{i=1}^{t}Var_{i-1}(d_{i})$ and $d_{max}=\max_{1\leqslant i\leqslant t}|d_{i}|$, we immediately obtain the above inequality as desired.

for any $\delta>0$, with probability at least $1-\delta$ over $(\mathcal{D}^{t},\bm{\omega}^{t})$,

where $\hat{B}^{2}_{2,i}=2M^{2}(\kappa+\phi)^{2}\ln\left(\frac{2t}{\delta}\right)\sum_{j=1}^{i-1}|a^{j}_{i-1}|^{2}$.

This is because, for any $1\leqslant i\leqslant t$, recall that from analysis in (3), we have $\mathcal{R}_{i}\leqslant\kappa^{1/2}L|f_{t}(x_{t})-h_{t}(x_{t})|\cdot\|h_{t}-f_{*}\|_{\mathcal{H}}$, therefore from Lemma 9,

Taking the sum over $i$, we therefore get

Applying these lemmas immediately gives us Theorem 4 and Theorem 5, which implies pointwise distance between the solution $f_{t+1}(\cdot)$ and $f_{*}(\cdot)$. Now we prove similar bounds in the sense of $L_{\infty}$ and $L_{2}$ distance.

Theorem 4 also implies a bound in $L_{\infty}$ sense, namely,

Consequently, for the average solution $\hat{f}_{t+1}(\cdot):=\frac{1}{t}\sum_{i=1}^{t}f_{i}(\cdot)$, we also have

This is because $\left\|f_{t+1}-f_{\ast}\right\|_{\infty}=\max_{x\in\mathcal{X}}|f_{t+1}(x)-f_{\ast}(x)|=|f_{t+1}(x_{\ast})-f_{\ast}(x_{\ast})|$, where $x_{*}\in\mathcal{X}$ always exists since $\mathcal{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 ).

With the choices of $\gamma_{t}$ in Lemma 9, we have

$\|f_{t+1}-f_{*}\|_{2}^{2}\leqslant\frac{C^{2}\ln(8\sqrt{e}t/\delta)+2\kappa Q_{2}^{2}\ln(2t/\delta)\ln^{2}(t)}{t},$ with probability at least $1-3\delta$ over $(\mathcal{D}^{t},\bm{\omega}^{t})$.

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

From Lemma 9, with probability at least $1-2\delta$, we have

From Lemma 7, for any $x\in\mathcal{X}$, we have

Since $C^{2}=4(\kappa+\phi)^{2}M^{2}\theta^{2}$, the above inequality can be writen as

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

From the analysis in Lemma 7, we also have that $|f_{t+1}(x)-h_{t+1}(x)|\leqslant C^{2}$. Therefore, with probability at least $1-\delta$ over $(\mathcal{D}^{t},\bm{\omega}^{t})$, we have

Let $\epsilon=\frac{\delta}{4t}$, we have

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

Proof By the Lipschitz continuity of $l(\cdot,y)$ and Jensen’s Inequality, we have

Then the theorem follows from Corollary 12.

Appendix C Suboptimality

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

If we set $\gamma_{t}=\frac{\theta}{t}$ with $\theta\nu=1$, then the average solution $\hat{f}_{t+1}:=\frac{1}{t}\sum_{i=1}^{t}f_{i}$ satisfies

where $Q=(4\kappa M^{2}+2\sqrt{2}\kappa^{1/2}LM(\kappa+\phi)Q_{1})/\nu$, with $Q_{1}$ defined as in Lemma 9.

Proof From the anallysis in Lemma 9,we have

Invoking strongly convexity of $R(f)$, we have $\langle h_{t}-f_{*},\bar{g}_{t}\rangle\geqslant R(h_{t})-R(f_{*})+\frac{\nu}{2}\|h_{t}-f_{*}\|_{\mathcal{H}}^{2}$. Taking expectaion on both size and use the bounds in last lemma, we have

Assume $\gamma_{t}=\frac{\theta}{t}$ with $\theta=\frac{1}{\nu}$, then cumulating the above inequalities leads to

Suppose the sequence $\{\Gamma_{t}\}_{t=1}^{\infty}$ satisfies $\Gamma_{1}\geqslant 0$, and $\forall t\geqslant 1$

where $\eta>1,\beta_{1},\beta_{2}>0$. Then $\forall t\geqslant 1$,

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 $t\geqslant 1$, i.e., $\Gamma_{t}\leqslant\frac{R}{t}$, then we have

where the last step can be verified as follows.

where the last step follows from the defintion of $R_{0}$.

Suppose the sequence $\{\Gamma_{t}\}_{t=1}^{\infty}$ satisfies

where $\beta_{1},\beta_{2},\beta_{3},\beta_{4}>0$ and $\delta\in(0,1/e)$. Then $\forall 1\leqslant j\leqslant t(t\geqslant 4)$,

Proof The proof follows by induction. When $j=1$ it is trivial. Let us assume it holds true for $1\leqslant j\leqslant t-1$, therefore,

Since $\sqrt{R}\geqslant 2\beta_{2}+2\sqrt{2}\beta_{3}+\sqrt{(2\beta_{2}+2\sqrt{2}\beta_{3})^{2}+\beta_{1}+\beta_{4}}$, we have $(2\beta_{2}+2\sqrt{2}\beta_{3})\sqrt{R}+\frac{\beta_{1}}{2}+\frac{\beta_{4}}{2}\leqslant R/2$. Hence, $\Gamma_{j+1}\leqslant\frac{R\ln(2t/\delta)\ln^{2}(t)}{j+1}$.

Appendix 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 $\mathcal{A}$ defined in (10). We first demonstrate that the operator $\mathcal{A}$ is the solution to the following optimization problem

where $\|\cdot\|_{HS}$ is the Hilbert-Schmidt norm of the operator. The gradient of $R(\mathcal{A})$ with respect to $\mathcal{A}$ is

Set $\nabla R(\mathcal{A})=0$, we could obtain the optimal solution, \mathcal{C}\big{(}\mathcal{C}+\frac{\sigma^{2}}{n}I\big{)}^{-1}, exactly the same as (10).

where $\widehat{\mathcal{C}}=k(x_{i},\cdot)\otimes k(x_{i},\cdot)$ which leads to update

With such update rule, we could show that $\mathcal{A}_{t+1}=\sum_{i=1,j\geqslant i}^{t}\beta_{ij}^{t+1}k(x_{i},\cdot)\otimes k(x_{j},\cdot)$ by induction. Let $\mathcal{A}_{1}=0$, then, $\mathcal{A}_{2}=\gamma_{1}k(x_{1},\cdot)\otimes k(x_{1},\cdot)$. Assume at $t$-th iteration, $\mathcal{A}_{t}=\sum_{i=1,j\geqslant i}^{t-1}\beta_{ij}^{t}k(x_{i},\cdot)\otimes k(x_{j},\cdot)$, and notice that

we have $\mathcal{A}_{t+1}=\sum_{i=1,j\geqslant i}^{t}\beta_{ij}^{t+1}k(x_{i},\cdot)\otimes k(x_{j},\cdot)$ where

Therefore, inductively, we could approximate $\mathcal{A}_{t+1}$ by