Data-driven Random Fourier Features using Stein Effect

1 Introduction

Kernel methods offer a comprehensive collection of nonparametric techniques for modeling a broad range of problems in machine learning. However, standard kernel machines such as kernel SVM or kernel ridge regression suffer from slow training and prediction which limits their use in real-world large-scale applications. Consider, for example, the training of linear regression over n labeled data points {(xi, yi)}ni=1 where xi ∈ Rd with n ≫ d and yi is a binary label. The time complexity of standard least square fit is O(nd2 ) and the memory demand is O(nd). Its kernel-based counterpart requires solving a linear system with an n-by-n kernel matrix, which takes O(n3 + n2d) time and O(n2 ) memory usage as typical. Such complexity is far from desirable in big-data applications.

In recent years, significant efforts have been devoted into low-rank kernel approximation for enhancing the scalability of kernel methods. Among existing approaches, random features [Rahimi and Recht, 2007; Rahimi and Recht, 2009; Le et al., 2013; Bach, 2015] have drawn considerable attention and yielded state-of-the-art results in classification accuracy on benchmark data sets [Huang et al., 2014; Dai et al., 2014;

Li and P´oczos, 2016]. Specifically, inspired fromBochner's theorem [Rudin, 2011], random Fourier features have been studied for evaluating the expectation of shift-invariant kernels (i.e., k(x, x′ ) = g(x-x′ ) for some function g). Rahimi and Recht [2007] proposed to use Monte-Carlo methods (MC) to estimate the expectation; Yang et al. [2014] leveraged the low-discrepancy properties of Quasi-Monte Carlo (QMC) sequences to reduce integration errors.

Both the MC and QMC methods rely on the (implicit) assumption that all the M random features are equally important, and hence assign a uniform weight 1M to the features in kernel estimation. Such a treatment, however, is arguably sub-optimal for minimizing the expected risk in kernel approximation. Avron et al. [2016b] presented a weighting strategy for minimizing a loose error bound which depends on the maximum norm of data points. On the other hand, Bayesian Quadrature (BQ) is a powerful framework that solves the integral approximation by E[f(x)] ≈ ∑m βmf(xm) where f(xm) is feature function and βm is with non-uniform weights. BQ leverages Gaussian Process (GP) to model the prior of function f(·), and the resulted estimator is Bayes optimal [Husz´ar and Duvenaud, 2012]. Thus, BQ could be considered a natural choice for kernel approximation. Nonetheless, a wellknown limitation (or potential weakness) is that the performance of Bayesian-based models heavily relies on extensive hyper-parameter tuning, e.g., on the condition of the covariance matrix in the Gaussian prior, which is not suitable for large real-world applications.

To address the fundamental limitations of the aforementioned work, we propose a novel approach to data-driven feature weighting in the approximation of shift-invariant kernels, which motivated by the by Stein Effect in the statistical literature (Section 4), and solve it using an efficient stochastic algorithm with a convex optimization objective. We also present a natural extension of BQ to the applications of kernel approximation (Section 3). The adapted BQ together with standard MC and QMC methods serve as representative baselines in our empirical evaluations on six benchmark data sets. The empirical results (Section 5) show that the proposed SteinEffect Shrinkage (SES) estimator consistently outperforms the baseline methods in terms of lowering the kernel approximation errors, and was competitive to or better than the best performer among the baseline methods in most task-oriented evaluations (on supervised learning of regression and classification tasks).

2 Preliminaries

Let us briefly outline the related background and key concepts in random-feature based kernel approximation.

2.1 Reproducing Kernel Hilbert Space

At the heart of kernel methods lies the idea that inner products in high-dimensional feature spaces can be computed in an implicit form via a kernel function k : X × X → R, which is defined on a input data domain X ⊂ Rd such that

where φ : X → H is a feature map that associates kernel k with an embedding of the input space into a high-dimensional Hilbert space H. A key to kernel methods is that as long as kernel algorithms have access to k, one need not to represent φ(x) explicitly. Most often, that means although φ(x) can be high-dimensional or even infinite-dimensional, their inner products, can be evaluated by k. This idea is known as the kernel trick.

The bottleneck of kernel methods lies in both training and prediction phase. The former requires to compute and store the kernel matrix K ∈ Rn×n for n data points, while the latter need to evaluate the decision function f(x) via kernel trick f(x) = ∑Ni=1 αik(xi, x), which sums over the nonzero support vectors αi. Unfortunately, Steinwart and Christmann [2008] show that the number of nonzero support vectors growth linearly in the size of training data n. This posts a great challenge for designing scalable algorithms for kernel methods in large-scale applications. Some traditional way to address this difficulty is by making trade-offs between time and space. [Kivinen et al., 2004; Chang and Lin, 2011].

2.2 RandomFeatures (RF)

Rahimi and Recht [2007] propose a low-dimensional feature map z(x) : X →R2M for the kernel function k

under the assumption that k fall in the family of shift-invariant kernel. The starting point is a celebrated result that characterizes the class of positive definite functions:

Theorem 1 (Bochner's theorem [Rudin, 2011]). A continuous, real valued, symmetric and shift-invariant function k on Rd is a positive definite kernel if and only if there is a positive finite measure P(w) such that

The most well-known kernel that belongs to the shift-invariant family is Gaussian kernel k(x, x′ ) = exp(- ‖x-x′ ‖2 2σ2 ), the associated density P(w) is again Gaussian, N(0, σ-2 Id).

Rahimi and Recht [2007] approximate the integral representation of the kernel (2) by Monte Carlo method as follows,

where

and w1, . . . ,wM are i.i.d. samples from P(w).

After obtaining the randomized feature map z, training data could be transformed into {(z(xi), yi}ni=1 . As long as M is sufficiently smaller than n, this leads to more scalable solutions, e.g., for regression we get back to O(nM 2 ) training and O(Md) prediction time, with O(nM) memory requirements. We can also apply random features to unsupervised learning problems, such as kernel clustering [Chitta et al., 2012].

2.3 Quasi-Monte Carlo Technique

Instead of using plain MC approximation, Yang et al. [2014] propose to use the low-discrepancy properties of Quasi-Monte Carlo (QMC) sequences to reduce the integration error in approximations of the form (3). Due to space limitation, we restrict our discussion to the background that is necessary for understanding subsequent sections. We refer interested readers to the comprehensive reviews [Caflisch, 1998; Dick et al., 2013] for more detailed exposition.

The QMC method is generally applicable to integrals over a unit cube. So the procedure is to first generate a discrepancy sequence t1, . . . , tM ∈ [0, 1]d , and transform it into a sequence w1, . . . ,wM in Rd . To convert the integral presentation of kernel (2) to an integral over the unit cube, a simple change of variables suffices. For t ∈ Rd , define Φ-1 (t) = [Φ-1 1 (t1), . . . ,Φ-1 d (td)] ∈ Rd , where we assume that the density function in (2) can be written as p(w) = ∏Mi=1 pj(wj) and φj is the cumulative distribution function (CDF) of pj , j = 1, . . . , d. By setting w = Φ-1 (t), the integral (2) is equivalent to

3 Bayesian Quadrature for RandomFeatures

Random features constructed by either MC [Rahimi and Recht, 2007] or QMC [Yang et al., 2014] approaches employ equal weights 1/M on integrand functions as shown in (3). An important question is how to construct different weights on the integrand functions to have a better kernel approximation.

Given a fixed QMC sequence, Avron et al. [2016b] solve the weights by optimizing a derived error bound based on the observed data. There are two potential weakness of [Avron et al., 2016b]. First, Avron et al. [2016b] does not fully utilize the data information because the optimization objective only depends on the upper bound information of |xi| instead of the distribution of xi. Second, the error bound used in Avron et al. [2016b] is potentially loose. Due to technical reasons, Avron et al. [2016b] relaxes hx,x′(w) to hx,x′(w)sinc(Tw) when they derive the error bound to optimize. It is not studied whether the relaxation results in a proper upper bound.

On the other hand, selecting weights by such way is closely connected to the Bayesian Quadrature (BQ), originally proposed by Ghahramani and Rasmussen [2002] and theoretically guaranteed under Bayes assumptions. BQ is a Bayesian approach that puts prior knowledge on functions to solve the integral estimation problem. We first introduce the work of Ghahramani and Rasmussen [2002], and then discuss how to apply BQ to kernel approximations.

Standard BQ considers a single integration problem with form

To utilize the information of f , BQ puts a Gaussian Process (GP) prior on f with mean 0 and covariance matrix KGP where

After conditioning on sample w1, · · · , wM from P, we obtain a closed-form GP posterior over f and use the posterior to get the optimal Bayes estimator ˆ f 1 as

where γm = ∫ KGP (w,wm)dP(w), βm BQ = (K-1 GP γ)m, and f(W) = [f(w1), · · · , f(wM)]. Clearly, the resulting estimator is simply a new weighted combination of the empirical observations f(wm) by β(m) BQ derived from the Bayes perspective.

In kernel approximation, we could treat it as a series of the single integration problem. Therefore, we could directly apply BQ on it. Here we discuss the technical issues of applying BQ on kernel approximation. First, we need to evaluate γ. For Gaussian kernel, we derive the closed-form expression

where γm is a random variable evaluated at wm.

3.1 Data-driven by Kernel Hyperparameters

The most import issue of applying BQ on kernel approximation is that it models the information of hx,x′ by the kernel KGP . In practice, we usually use Gaussian kernel for KGP , then the feature information of (x, x′ ) is modeling by the kernel bandwidth. As suggested by Rasmussen [2006], one could tune the best hyper-parameter either by maximizing the marginal likelihood or cross-validation based on the metrics in interested. Note that the existing BQ approximates the integral representation of a kernel function value, evaluated at a single pair of data (x, x′ ). However, it is unfeasible to tune individual σGP for n2 pairs of data. Therefore, we tune a global hyper-parameter σGP on a subset of pair data points.

1 For the square loss.

4 Stein Effect Shrinkage (SES) Estimator

In practice, it is well known that the performance of Gaussian Process based models heavily depend on hyper-parameter σGP , i.e. the bandwidth of kernel function kGP . Therefore, instead of using Bayesian way to obtain the optimal Bayes estimator from BQ, we start from the risk minimization perspective by considering the Stein effect to derive the shrinkage estimator for random features.

The standard estimator of Monte-Carlo methods is the unbiased empirical average using equal weights 1M for M samples which sum to 1. However, Stein effect (or James-Stein effect) suggests a family of biased estimators could result in a lower risk than the standard empirical average. Here we extend the analysis of non-parametric Stein effect to the kernel approximation with random features by considering the risk function R(k, k′ ) = Ex,x′ ,w(k(x - x′ ) - k′ (x - x′ ; w))2 .

Theorem2. (Stein effect on kernel approximation) Let ˆ k(x - x′ ,w) = 1M ∑Mm=1 hx,x′ (wm). For any estimator µ, which is independent of x and w, there exists 0 ≤ α < 1 such that ˜ k(·) = αµ + (1 - α)ˆ k(·) is an estimator with lower risk R(k, ˜ k) < R(k, ˆ k).

Proof. The proof closely follows Muandet et al. [2014] for a different problem under the non-parametric setting 2 . By bias-variance decomposition, we have

By elementary calculation, we have R(k, ˜ k) < R(k, ˆ k) when

and the optimal value happened at

Corollary 2.1. (Shrinkage estimator) If k(x-x′ ) > 0, there is a shrinkage estimator ˜ k(·) = (1 - α∗ )ˆ k(·) that has lower risk than the empirical mean estimator ˆ k(·) by setting µ = 0.

The standard shrinkage estimator derived from Stein effect uses the equal weights 1-α M to minimize the variance; however, a non-uniform weights usually results in better performance in practice [Husz´ar and Duvenaud, 2012]. Therefore, based on the analysis of Theorem 2, we proposed the estimator ˜ k(x - x′ , w) = ∑Mm=1 βmh(wm) obtained by minimizing the empirical risk with the constraint ∑Mm=1 β2m ≤ c to shrink the weights, where c is a small constant. With Lagrangian multiplier, we end up the following risk minimization formulation which is data-driven by directly taking k(x, x′ ) into account,

2 The original Stein effect has a Gaussian distribution assumption.

where ◦ denotes the Hadamard product, and λ is the regularization coefficient to shrink the weights.

The objective function (6) can be further simplified as a least squares regression (LSR) formulation,

where t ∈ Rn2 ,Z ∈ Rn2×M , such that for all i, j ∈ P = {1, . . . , n}, there exist a corresponding pair mapping function ξ : ξP (i, j) → k satisfying tk := φ(xi)T φ(xj) and Z(k, : ) := z(xi) ◦ z(xj).

We approximately solve optimization problem (7) via the matrix sketching techniques [Avron et al., 2013; Avron et al., 2016a; Woodruff and others, 2014]:

where S ∈ Rr×n2 is the randomized sketching matrix. Solving the sketched LSR problem (8) now costs only O(rd2 +Ts) where Ts is the time cost of sketching. From the optimization perspective, suppose β∗ be the optimal solution of LSR problem (7), and ˜ β be the solution of sketched LSR problem (8). [Wang et al., 2017] prove that if r = O(M/ϵ + poly(M)), the the objective function value of (8) at ˜ β is at most ϵ worse than the value of (7) at β∗ . In practice, we consider an efficient Ts = O(r) sampling-based sketching matrix where S is a diagonal matrix with Si,i = 1/pi if row i was include in the sample, and Si,i = 0 otherwise. That being said, we form a sub-matrix of Z by including each row of Z in the sub-matrix independently with probability pi. [Woodruff and others, 2014] suggests to set pi proportional to ‖Zi‖, the norm of ith row of Z, to have a meaningful error bound.

Algorithm 1 Data-driven random feature (SES)

• 1: procedure SES
• 2: Compute the Fourier transform p of the kernel
• 3: Generate a sequence w1, . . . ,wM by QMC or MC
• 4: Compute un-weighted random feature z(x) by (4)
• 5: Calculate shrinkage weight by (8)
• 6: Compute data-driven random features by

5 Experiments

In this section, we empirically demonstrate the benefits of our proposed method on six data sets listed in Table 1. The features in all the six data sets are scaled to zero mean and unit variance in a preprocessing step.

We consider two different tasks: i) the quality of kernel approximation and ii) the applicability in supervised learning. Both the relative kernel approximation errors and the regression/classification errors are reported on the test sets.

5.1 Methods for Comparison

• MC: Monte Carlo approximation with uniform weight.
• QMC: Quasi-Monte Carlo approximation with uniform weight. We adapt the implementation of Yang et al. [2014] on Github. 3
• BQ: QMC sequence with Bayesian Quadrature weights. We modify the source code from Husz´ar and Duvenaud [2012]. 4

Gaussian RBF kernel is used through all the experiments where the kernel bandwidth σ is tuned on the set {2-10 , 2-8 , . . . , 28 , 210 } that is in favor of Monte Carlo methods. For QMC method, we use scrambling and shifting techniques recommended in the QMC literature (See Dick et al. [2013] for more details). In BQ framework, we consider the GP prior with the covariance matrix kGP (w,w′ ) = exp(- ‖w-w′ ‖22 2σ2GP ) where σGP is tuned on the set {2-8 , 2-6 , . . . , 26 , 28 } over a subset of sampled pair data (x, x′ ). Likewise, regularization coefficient λβ of in the proposed objective function (6) is tuned on the set {2-8 , 2-6 , . . . , 26 , 28 } over the same subset of sampled pair data. For regression task, we solve ridge regression problems where the regularization coefficient λ is tuned by five folds cross-validation. For classification task, we solve L2-loss SVM where the regularization coefficient C is tuned by five folds cross-validation. All experiments code are written in MATLAB and run on a Intel(R) Xeon(R) CPU 2.40GHz Linux server with 32 cores and 190GB memory.

5.2 Quality of Kernel Approximation

In our setting, the most natural metric for comparison is the quality of approximation of the kernel matrix. We use the relative kernel approximation error ‖K- ˜K‖F /‖K‖F to measure the quality, as shown in Figure 1.

SES outperforms the other methods on cpu, census, adult, mnist, covtype, and is competitive with the best baseline methods over the years data sets. Notice that SES is particularly strong when the number of random features is relatively small. This phenomenon confirms our theoretical analysis, that is, the smaller the random features are, the larger the variance would be in general. Our shrinkage estimator with Stein effect enjoys better performance than the empirical mean estimators (MC and QMC) using uniformly weighted features without shrinkage. Although BQ does use weighted random features, its performance is still not as good as SES because tuning the covariance matrix in the Gaussian Process is rather difficult. Moreover, it is even harder to tune a global bandwidth for all integrand functions that approximate the kernel function.

3 https://github.com/chocjy/QMC-features

4 http://www.cs.toronto.edu/ duvenaud/

5.3 Supervised Learning

For regression task, we report the relative regression error ‖y-˜ y‖2/‖y‖2. For classification task, we report the classification error.

We investigate if better kernel approximations reflect better predictions in the supervised learning tasks. The results are shown in Table 2. We can see that among the same data set (cpu, census, adult, and covtype) our method performs very well in terms of kernel approximation, and it also yields lower errors in most of the cases of the supervised learning tasks. Note that we only present partial experiment results on some selected number of random features. However, we also observe that the generalization error of our proposed method may be higher than that of other state-of-the-art methods at some other number of random features that we did not present here. This situation is also observed in other works [Avron et al., 2016b]: better kernel approximation does not always result in better supervised learning performance. Under certain conditions, less number of features (worse kernel approximation) can play a similar role of regularization to avoid overfitting [Rudi et al., 2016]. Rudi et al. [2016] suggest to solve this issue by tuning number of random features.

5.4 Distributions of FeatureWeights

We are interested in the distributions of random feature weights derived from BQ and SES. For different numbers of random features (M = 32, 128, 512), we plotted the corresponding histograms and compared them with the equal weight of 1/M used in MC and QMC, as shown in Figure 2. When the number of random features is small, the difference between BQ and SES is big. As the number of random features increases, the weights of each method gradually converge to a stable distribution, and get increasingly closer to the uniform weight spike. This empirical observation agrees with the theoretical analysis we found, namely with the role of the weight estimator in the bias-variance trade-off.

5.5 Behavior of Randomized Optimization Solver

In theory, our randomized solver only needs to sample O(M) data in our proposed objective function (8) to have a sufficiently accurate solution. We verify this assertion by varying the sampled data size r in sketching matrix S and examine it affects the relative kernel approximation errors. The results are shown in Figure 3. For census data set, we observe that as the number of sampled data r exceeds the number of random features M, the kernel approximation of our proposed method outperforms the state-of-the-art, QMC. Furthermore, the number of sampled data required by our proposed randomized solver is still the same order of M even for a much larger data set covtype, with 0.5 millions data points. On the other hand, BQ does not have consistent behavior regarding the sampled data size, which is used for tuning the hyperparmaeter σGP . This again illustrates that BQ is highly sensitive to the covariance matrix and is difficult to tune σGP on only a small subset of sampled data. In practice, we found that we usually only need to sample up to 2 ∼ 4 times the number of random features M to get a relatively stable and good minimizer β∗ . This observation coincides with the complexity analysis of our proposed randomized solver in Section 4.

5.6 Uniform v.s. Non-uniform ShrinkageWeight

We conducted experiments with SES using a uniform shrinkage weight versus using non-uniform shrinkage weights on the cpu and census data sets. The results in Figures 4 confirm our conjecture that non-uniform shrinkage weights are more beneficial, although both enjoy the lower empirical risk from the Stein effect. We also observed similar behavior on other data sets, and we omit those results.

6 Discussion and Conclusion

Nystr¨om methods [Williams and Seeger, 2001; Drineas and Mahoney, 2005] is yet another line of research for kernel approximation. Specifically, eigendecomposition of the kernel matrix is achieved via low rank factorization with various sampling techniques [Wang and Zhang, 2013]. We refer interested readers to the comprehensive study [Yang et al., 2012] for more detailed comparisons. Due to the page limits, the

(c) M = 64, covtype data set (d) M = 256, covtype data set

Figure 4: Comparison of uniform and non-uniform shrinkage weight.

scope of this paper only focus on the study of random Fourier features.

Finally, we notice [Sinha and Duchi, 2016] proposes a weighting scheme for random features which matches the label correlation matrix with the goal to optimize the supervise objective (e.g. classification errors). In contrast, our work aims to optimize the kernel approximation, which does not require the label information to solve weights and can be applied to different tasks after solving the weights once.

To conclude, we propose a novel shrinkage estimator, which enjoys the benefits of Stein effect w.r.t. lowering the empirical risk compared to empirical mean estimators. Our estimator induces non-uniform weights on the random features with desirable data-driven properties. We further provide an efficient randomized solver to optimize this estimation problem for real-world large-scale applications. The extensive experimental results demonstrate the effectiveness of our proposed method.

Acknowledgments

We thank the reviewers for their helpful comments. This work is supported in part by the National Science Foundation (NSF) under grants IIS-1546329 and IIS-1563887.