Title: QUADRATURE-BASED FEATURES FOR KERNEL APPROXIMATION

Abstract: We consider the problem of improving kernel approximation via randomized feature maps. These maps arise as Monte Carlo approximation to integral representations of kernel functions and scale up kernel methods for larger datasets. We propose to use more efficient numerical integration technique to obtain better estimates of the integrals compared to the state-of-the-art methods. Our approach allows the use of information about the integrand to enhance approximation and facilitates fast computations. We derive the convergence behavior and conduct an extensive empirical study that supports our hypothesis.

1. INTRODUCTION:

Kernel methods proved to be an efficient technique in numerous real-world problems. The core idea of kernel methods is the kernel trick -compute an inner product in a high-dimensional (or even infinite-dimensional) feature space by means of a kernel function k: k(x, y) = ψ(x), ψ(y) ,

where ψ : X → F is a non-linear feature map transporting elements of input space X into a feature space F. It is a common knowledge that kernel methods incur space and time complexity infeasible to be used with large-scale datasets directly. For example, kernel regression has O(N 3 + N d 2 ) training time, O(N 2 ) memory, O(N d) prediction time complexity for N data points in original d-dimensional space X . One of the most successful techniques to handle this problem (Rahimi & Recht (2008) ) introduces a low-dimensional randomized approximation to feature maps:

EQUATION

This is essentially carried out by using Monte-Carlo sampling to approximate scalar product in (1) . A randomized D-dimensional mapping Ψ(•) applied to the original data input allows employing standard linear methods, i.e. reverting the kernel trick. In doing so one reduces the complexity to that of linear methods, e. 2016)) aims to improve the convergence of approximation so that a smaller D can be used to obtain the same quality of approximation.

This paper considers kernels that allow the following integral representation k(x, y) = E q(w) g xy (w) ≈ E p(w) f xy (w) = I(f xy ), p(w) = 1 (2π) d/2 e -w 2 2 ,

(3) where q(w) is a density associated with a kernel, e.g. the popular Gaussian kernel has q(w) = p(w), so the exact equality holds with g xy (w) = f xy (w) = φ(w x) φ(w y), where φ(•) = [cos(•), sin(•)] .

The class of kernels admitting the form in (3) covers shift-invariant kernels (e.g. radial basis function (RBF) kernels) and Pointwise Nonlinear Gaussian (PNG) kernels. They are widely used in practice and have interesting connections with neural networks (Cho & Saul (2009) , Williams (1997)).

The main challenge for the construction of low-dimensional feature maps is the approximation of the expectation in (3) which is d-dimensional integral with Gaussian weight. While standard Monte-Carlo rule is easy to implement, there are better quadrature rules for such kind of integrals. For example, Yang et al. (2014) apply quasi-Monte Carlo (QMC) rules and obtain better quality kernel matrix approximations compared to random Fourier features of Rahimi & Recht (2008) .

Unlike other research studies we refrain from using simple Monte Carlo estimate of the integral, instead, we propose to use specific quadrature rules. We now list our contributions:

1. We propose to use advanced quadrature rules to improve kernel approximation accuracy. We also provide an analytical estimate of the error for the used quadrature rules. 2. We note that for kernels with specific integrand f xy (w) in (3) one can improve on its properties. For example, for kernels with even function f xy (w) we derive the reduced quadrature rule which gives twice smaller embedded dimension D with the same accuracy. This applies, for example, to any RBF kernel. 3. We use structured orthogonal matrices (so-called butterfly matrices) when designing quadrature rule that allow fast matrix by vector multiplications. As a result, we speed up the approximation of the kernel function and reduce memory requirements. 4. We demonstrate our approach on a set of regression and classification problems. Empirical results show that the proposed approach has a better quality of approximation of kernel function as well as better quality of classification and regression when using different kernels.

2. QUADRATURE-BASED RANDOM FEATURES:

We start with rewriting the expectation in equation (3) as integral of f xy with respect to p(w):

I(f xy ) = (2π) -d 2 ∞ -∞ • • • ∞ -∞ e -w w 2 f xy (w)dw.

Integration can be performed by means of quadrature rules. The rules usually take a form of interpolating function that is easy to integrate. Given such a rule, one may sample points from the domain of integration and calculate the value of the rule at these points. Then, the sample average of the rule values would yield the approximation of the integral.

We use the average of sampled quadrature rules developed by Genz & Monahan (1998) to yield unbiased estimates of I(f xy ). A change of coordinates is the first step to facilitate stochastic sphericalradial rules. Now, let w = rz, with z z = 1, so that w w = r 2 for r ∈ [0, ∞], leaving us with

EQUATION

where

U d = {z : z z = 1, z ∈ R d }.

As mentioned earlier we are going to use a combination of radial R and spherical S rules. We now describe the logic behind the used quadratures.

Stochastic radial rules. Stochastic radial rule R(h) of degree 2l + 1 has the form of weighted symmetric sums:

R(h) = l i=0 w i h(ρ i ) + h(-ρ i ) 2 ,

where h is an integrand in infinite range integral

T (h) = ∞ -∞ e -r 2 2 |r| d-1 h(r)dr.

To get an unbiased estimate for T (h), points ρ i are sampled from specific distributions which depend on the degree of the rule. Weights w i are derived so that R has a polynomial degree 2l + 1, i.e. is exact for integrands h(r) = r p with p = 0, 1, . . . , 2l + 1. For radial rules of degree three R 3 the point ρ 0 = 0, while ρ 1 ∼ χ(d + 2) follows Chi-distribution with d + 2 degrees of freedom. Higher degrees require samples from more complex distributions which are hard to sample from.

Stochastic spherical rules. Spherical rule S(s) approximates an integral of a function s(z) over the surface of unit d-sphere U d and takes the following form:

S(s) = p j=1 w j s(z j ),

where z j are points on U d , i.e. z z = 1. If we set weight w j = |U d | 2(d+1) and sum function s values at original and reflected vertices v j of randomly rotated d-simplex V, we will end up with a degree three rule:

S 3 Q (f xy (ρQz)) = |U d | 2(d + 1) d+1 j=1 f xy (-Qv j ) + f xy (Qv j ) ,

where v j is the j'th vertex of d-simplex V with vertices on U d and Q is a random d × d orthogonal matrix. We justify the choice of the degree in Appendix A.

Since the value of the integral is approximated as the sample average, the key to unbiased estimate is proper randomization. In this case, randomization is attained through the matrix Q. It is crucial to generate uniformly random orthogonal matrices to achieve an unbiased estimate for spherical surface integrals. We consider various designs of such matrices further in Section 3.

Stochastic spherical-radial rules. Meanwhile, combining foregoing rules results in stochastic spherical-radial rule of degree three:

EQUATION

which we finally apply to the approximation of (4) by averaging the samples of SR 3,3 Q,ρ :

EQUATION

where n is the number of sampled SR rules. Speaking in terms of approximate feature maps, the new feature dimension D in case of quadrature based approximation equals 2n(d + 1) as we sample n rules and evaluate each of them at 2(d + 1) points. Surprisingly, empirical results (see Section 5) show that even a small number of rule samples n provides accurate approximations.

Properties of the integrand. We also note here that for specific functions f xy (w) we can derive better versions of SR rule by taking on advantage of the knowledge about the integrand. For example, the Gaussian kernel has f xy (w) = cos(w (xy)). Note that f is even, so we can discard an excessive term in the summation in (5), since f (w) = f (-w), i.e SR 3,3 rule reduces to

EQUATION

)

Variance of the error. We contribute the variance estimation for the stochastic spherical-radial rules when applied to kernel function. To the best of our knowledge, it has not been done before. In case of kernel functions the integrand f xy can be represented as

EQUATION

where z 1 , z 2 are scalar values. Using this representation of kernel function and its Taylor expansion we can obtain the following proposition (see Appendix B for detailed derivation of the result): Proposition 2.1. The quadrature rule ( 6) is an unbiased estimate of integral of any integrable function f . If function f can be represented in the form (8), i.e. f (w) = g(z 1 , z 2 ), z 1 = w x, z 2 = w y for some x, y ∈ R d , all 4-th order partial derivatives of g are bounded and D = 2n(d+1) is the number of generated features, then

V[ Î(f )] ≤ 2.66M 2 1 L 8 nd + 212M 1 M 2 L 6 nd 3 + (d + 95)M 2 2 L 4 4nd 2 (d -2)

,

where

M 1 = max sup z ∂ 4 g ∂z 4 1 , sup z ∂ 4 g ∂z 4 2 , sup z ∂ 4 g ∂z 2 1 ∂z 2 2 , M 2 = max j=0,1,2 ∂ 2 g ∂z j 1 ∂z 2-j 2 (0, 0)

and L = max{ x , y }.

Constants M 1 , M 2 in the proposition are upper bounds on the derivatives of function g and don't depend on the data set, while L plays the role of the scale of inputs. The proposition implies that the error of approximation is proportional to L -the less the scale, the better the accuracy (see Figure 1 ). However, scaling input vectors is equivalent to changing the parameters of the kernel function.

For example, decreasing the norm of input variables for RBF kernel is equivalent to increasing the kernel width σ. Therefore, the wide RBF kernels are approximated better than the narrow ones. This result also gives us the rate of convergence O( 1 /nd) for the quadrature rule. Figure 1 : Relative error of approximation of kernel matrix with 95% confidence interval depending on the scaling factor, Gaussian kernel was used. In this experiment for each scaling factor α we construct approximate kernel matrix K and the exact kernel matrix K using scaled input vectors x = x/α. To plot the confidence interval we run each experiment 10 times each time generating new weights. Experiment was conducted for 3 different input dimensions: d = 10, 100, 500.

The quadrature rule (5) grants us some freedom in the choice of random orthogonal matrix Q.

The next section discusses such matrices and suggests butterfly matrices for fast matrix by vector multiplication as the SR 3,3 rule implementation involves multiplication of the matrix QV by the data vector x.

3. GENERATING UNIFORMLY RANDOM ORTHOGONAL MATRICES:

Previously described stochastic spherical-radial rules require a random orthogonal matrix Q (see equation ( 5)). If Q follows Haar distribution on the set of all matrices in the orthogonal group O(d) in dimension d, then the averages of spherical rules S 3 Qi (s) provide unbiased degree three estimates for integrals over unit sphere. Essentially, Haar distribution means that all orthogonal matrices in the group are equiprobable, i.e. uniformly random. Methods for sampling such matrices vary in their complexity of generation and multiplication.

Techniques based on QR decomposition (Mezzadri (2006) ) have complexity cubic in d, and the resulting matrix does not allow fast matrix by vector multiplications. Another set of methods is based on a sequence of reflectors (Stewart (1980) ) or rotators (Anderson et al. (1987) ). The complexity is better (quadratic in d), however the resulting matrix is unstructured and, thus, implicates no fast matrix by vector multiplication. In Choromanski et al. (2017) random orthogonal matrices are considered. They are constructed as a product of random diagonal matrices and Hadamard matrices and therefore enable fast matrix by vector products. Unfortunately, they are not guaranteed to follow the Haar distribution.

To satisfy both our requirements, i.e low computational/space complexity and generation of Haar distributed orthogonal matrices, we propose to use so-called butterfly matrices.

Butterfly matrices. The method from Genz (1998) generates Haar distributed random orthogonal matrix B. As it happens to be a product of butterfly structured factors, a matrix of this type conveniently possesses the property of fast multiplication. For d = 4 an example of butterfly orthogonal matrix is

B (4) =    c 1 -s 1 0 0 s 1 c 1 0 0 0 0 c 3 -s 3 0 0 s 3 c 3       c 2 0 -s 2 0 0 c 2 0 -s 2 s 2 0 c 2 0 0 s 2 0 c 2    =    c 1 c 2 -s 1 c 2 -c 1 s 2 s 1 s 2 s 1 c 2 c 1 c 2 -s 1 s 2 -c 1 s 2 c 3 s 2 -s 3 s 2 c 3 c 2 -s 3 c 2 s 3 s 2 c 3 s 2 s 3 c 2 c 3 c 2   .

Definition 3.1. Let c i = cos θ i , s i = sin θ i for i = 1, . . . , d -1 be given. Assume d = 2 k with k > 0. Then an orthogonal matrix B (d) ∈ R d×d is defined recursively as follows

B (2d) = B (d) c d -B (d) s d B(d) s d B(d) c d , B (1) = 1,

where B(d) is the same as B (d) with indexes i shifted by d, e.g. 

B (2) = c 1 -s 1 s 1 c 1 , B(2) = c 3 -s 3 s 3 c 3 . Matrix B (d)

4. KERNELS:

This section gives examples on how quadrature rules can be applied to a number of kernels.

4.1. GAUSSIAN KERNEL:

Radial basis function (RBF) kernels are popular kernels widely used in kernel methods. Gaussian kernel is a widely exploited RBF kernel and has the following form:

k(x, y) = exp - x -y 2 2σ 2 .

In this case the integral representation has φ(w x) = [cos(w x), sin(w x)] . Since f xy (0) = 1, SR 3,3 rule for Gaussian kernel has the form (σ appears due to scaling):

SR 3,3 Q,ρ (f xy ) = 1 - d ρ 2 + d d + 1 d+1 j=1 f xy ρQv j σ ρ 2 ,

4.2. ARC-COSINE KERNELS:

Arc-cosine kernels were originally introduced by Cho & Saul (2009) upon studying the connections between deep learning and kernel methods. The integral representation of the b th -order arc-cosine kernel is

k b (x, y) = 2 R d φ b (w x)φ b (w y)p(w)dw,

where φ b (w x) = Θ(w x)(w x) b , Θ(•) is the Heaviside function and p is the density of the standard Gaussian distribution. Such kernels can be seen as an inner product between the representation produced by infinitely wide single layer neural network with random Gaussian weights. They have closed form expression in terms of the angle θ = cos -1 x y x y between x and y.

0 th -order arc-cosine kernel is given by k 0 (x, y) = 1 -θ π , 1 st -order kernel is given by k 1 (x, y) =

x y π (sin θ + (πθ) cos θ). Let φ 0 (w x) = Θ(w x) and φ 1 (w x) = max(0, w x), then we can rewrite the integral representation as follows:

k b (x, y) = 2 R d φ b (w x)φ b (w y)p(w)dw ≈ 2 n n i=1 SR 3,3

Qi,ρ i . For arccosine kernel of order 0 the value of the function φ 0 (0) = Θ(0) = 0.5 results in

SR 3,3 Q,ρ (f xy ) = 0.25 1 - d ρ 2 + d d + 1 d+1 j=1 f xy (ρQv j ) + f xy (-ρQv j ) 2ρ 2 .

In the case of arc-cosine kernel of order 1, the value of φ 1 (0) is 0 and the SR3,3 rule reduces to

SR 3,3 Q,ρ (f xy ) = d d + 1 d+1 j=1 f xy (|ρQv j |) 2ρ 2 .

4.3. EXPLICIT MAPPING:

The explicit mapping can be written as follows:

ψ(x) = [ a 0 f xy (0) a 1 f xy (w 1 x) . . . a D f xy (w D x) ]

where a 0 = 1 -d ρ2 

1 , for i = 1, . . . , D a i = 1 . To get D features one simply stacks n = D 2(d+1) such matrices

W j = ρ j (Q j V) -(Q j V)

so that W ∈ R 2n(d+1)×d , where only Q j ∈ R d×d and ρ j (j = 1, . . . , n) are generated randomly.

For example, in case of the Gaussian kernel the mapping can be rewritten as 2

ψ G (x) = [a 0 a 1 cos(w 1 x) . . . a D cos(w D x) a 1 sin(w 1 x) . . . a D sin(w D x)].

5. EXPERIMENTS:

We extensively study the proposed method on several established benchmarking datasets: Powerplant, LETTER, USPS, MNIST, CIFAR100 (Krizhevsky & Hinton (2009)), LEUKEMIA (Golub et al. (1999)). In Section 5.2 we show kernel approximation error across different kernels and number of features. We also report the quality of SVM models with approximate kernels on the same data sets in Section 5.3. The compared methods are described below.

5.1. METHODS:

We present a comparison of our method with estimators based on simple Monte Carlo 3 . The Monte Carlo approach has a variety of ways to generate samples: unstructured Gaussian Quadrature rules. Our main method that uses stochastic spherical-radial rules with Q = B4 (butterfly matrix) is denoted by B. As mentioned earlier we also include a variant of our algorithm that uses an orthogonal matrix Q based on a sequence of random reflectors (we denote it as H).

5.2. KERNEL APPROXIMATION:

To measure kernel approximation quality we use relative error in Frobenius norm K-K F K F , where K and K denote exact kernel matrix and its approximation. We run experiments for the kernel approximation on a random subset of a dataset (see Appendix D for details). Approximation was constructed for different number of SR samples n = D 2(d+1) , where d is an original feature space dimensionality and D is the new one. For the Gaussian kernel we set hyperparameter γ = 1 2σ 2 to the same value for all the approximants, while arc-cosine kernels have no hyperparameters.

We run experiments for each [kernel, dataset, n] tuple and plot 95% confidence interval around the mean value line. Figure 2 show results for kernel approximation error on LETTER, MNIST, CIFAR100 and LEUKEMIA datasets. Figure 3 : Accuracy/R 2 score using embeddings with three kernels (columns: arc-cosine 0, arccosine 1, Gaussian) on three datasets (rows: Powerplant, LETTER, USPS). Higher is better. The x-axis represents the factor to which we extend the original feature space, n = D 2(d+1) , where d is the dimensionality of the original feature space, D is the dimensionality of the new feature space. We drop one of our methods H here since its kernel approximation almost coincides with B.

We observe that for the most of the datasets and kernels the methods we propose in the paper (B, H) show better results than the baselines. They do coincide almost everywhere, which is expected, as the B method is only different from H in the choice of the matrix Q to facilitate speed up.

5.3. CLASSIFICATION/REGRESSION WITH NEW FEATURES:

We report accuracy and R 2 scores for the classification and regression tasks on the same data sets (see Figure 3 ). We examine the performance with the same setting (the number of runs for each [kernel, dataset, n] tuple) as in experiments for kernel approximation error, except now we map the whole dataset. We use Support Vector Machines to obtain predictions. We also drop one of our methods H here since its kernel approximation almost coincides with B.

Kernel approximation error does not fully define the final prediction accuracy -the best performing kernel matrix approximant not necessarily yields the best accuracy or R 2 score. However, the empirical results illustrate that our method delivers comparable and often the best quality on the final tasks. We also note that in many cases our method provides greater performance using less number of features n, e.g. LETTER and Powerplant datasets with arc-cosine kernel of the first order.

6. RELATED WORK:

The most popular methods for scaling up kernel methods are based on a low-rank approximation of the kernel using either data-dependent or independent basis functions. The first one includes Nyström method (Drineas & Mahoney (2005) ), greedy basis selection techniques (Smola & Schölkopf (2000) ), incomplete Cholesky decomposition (Fine & Scheinberg (2001) ).

The construction of basis functions in these techniques utilizes the given training set making them more attractive for some problems compared to Random Fourier Features approach. In general, data-dependent approaches perform better than data-independent approaches when there is a gap in the eigen-spectrum of the kernel matrix. The rigorous study of generalization performance of both approaches can be found in (Yang et al. (2012) ).

In data-independent techniques, the kernel function is approximated directly. Most of the methods (including the proposed approach) that follow this idea are based on Random Fourier Features (Rahimi & Recht (2008) ). They require so-called weight matrix that can be generated in a number of ways. Le et al. (2013) form the weight matrix as a product of structured matrices. It enables fast computation of matrix-vector products and speeds up generation of random features.

Another work (Felix et al. (2016) ) orthogonalizes the features by means of orthogonal weight matrix. This leads to less correlated and more informative features increasing the quality of approximation. They support this result both analytically and empirically. The authors also introduce matrices with some special structure for fast computations. 2015)) the authors make attempt to improve quality of the approximation of Random Fourier Features by optimizing sequences conditioning on a given dataset.

Among the recent papers there are works that, similar to our approach, use the numerical integration methods to approximate kernels. While Bach (2017) carefully inspects the connection between random features and quadratures, they did not provide any practically useful explicit mappings for kernels. Leveraging the connection Dao et al. (2017) propose several methods with Gaussian quadratures, among them three schemes are data-independent and one is data-dependent. The authors do not compare them with the approaches for random feature generation other than random Fourier features. The data-dependent scheme optimizes the weights for the quadrature points to yield better performance.

7. CONCLUSION:

In this work we proposed to apply advanced integration rule that allowed us to achieve higher quality of kernel approximation. Our derivation of the variance of the error implies the dependence of the error on the scale of data, which in case of Gaussian kernel can be interpreted as width of the kernel. However, as we have seen earlier, accuracy on the final task has no direct dependence on the approximation quality, so we can only speculate whether better approximated wide kernels deliver better accuracy compared to the poorer approximated narrow ones. It is interesting to explore this connection in the future work.

To speed up the computations we employed butterfly orthogonal matrices yielding the computational complexity O(d log d). Although the procedure we used to generate butterfly matrices claims to produce uniformly random orthogonal matrices, we found that it is not always so. However, the comparison of the method H (uses properly distributed orthogonal matrices) with method B (sometimes fails to do so) did not reveal any differences. We also leave it for the future investigation.

Our experimental study confirms that for many kernels on the most datasets the proposed approach delivers better kernel approximation. Additionally, the empirical results showed that the quality of the final task (classification/regression) is also higher than the state-of-the-art baselines. The connection between the final score and the kernel approximation error is to be explored as well.

Let us calculate the variance of the estimate.

EQUATION

Distribution of random variable ρ is χ(d + 2), therefore

E 1 ρ 2 = 1 d , E 1 ρ 4 = 1 d(d -2) , d > 2.

1. Now, let us calculate the variance VS i from the first term of equation ( 9)

EQUATION

where

M 1 = max sup z ∂ 4 g ∂z 4 1 (z 1 , z 2 ) , sup z ∂ 4 g ∂z 4 2 (z 1 , z 2 ) , sup z ∂ 4 g ∂z 2 1 ∂z 2 2 (z 1 , z 2 ) .

To calculate expectations of the form E(w x) k we will use expression for integral of monomial over unit sphere Baker (1997)

EQUATION

where k i = 2s i , s i ∈ Z + , |k| = i k i , σ d is a volume of an d-dimensional unit sphere. For example, let us show how to calculate E(w x) 4 :

E(w x) 4 = E i,j,k,m w i w j w k w m x i x j x k x m =

= thanks to symmetry all terms, for which at least one index doesn't coincide with other indices, are equal to 0. = E

  i w 4 i x 4 i + 3 i =j w 2 i w 2 j x 2 i x 2 j   = = 3 d(d + 2) i x 4 i + 3 d(d + 2) i =j x 2 i x 2 j = 3 d(d + 2) x 4 .

Using the same technique for (10) we obtain 

V(B i ) ≤ M 2 1 x 8

For the variance of A i we have

EQUATION

where

M 2 = max j=0,1,2 ∂ 2 g ∂z j 1 ∂z 2-j 2 (0, 0) .

Let's estimate covariance Cov(A i , B i ):

Cov(A i , B i ) ≤ E(A i B i ) + |EA i EB i |.

The first term of the right hand side:

E   2 ρ 2 2 j=0 c j (w T i x) j (w T i y) 2-j j!(2 -j)! 4 j=0 d i j (w T i x) j (w T i y) 4-j j!(4 -j)!   ≤ we assume that x ≥ y ≤ 210M 1 M 2 x 6 d 2 (d + 2)(d + 4

) .

The second term EA i EB i (again for x ≥ y )

|EA i EB i | ≤ 2M 1 M 2 x 6 d 3 (d + 2) .

Combining the derived inequalities we obtain

EQUATION

)

2. Now, let's examine the expectation of the second term in ( 9): Substituting ( 14) and ( 15) into ( 9) we obtain

EQUATION

V SR 3,3 Q,ρ (f xy ) ≤ (d + 1) 2.66M 2 1 x 8 d(d + 2) + 212M 1 M 2 x 6 d 3 (d + 2) + 24M 2 2 x 4 d 2 (d 2 -4) + + 1 4d(d -2) d(d -1) M 2 x 2 d 2 ≤ ≤ 2.66M 2 1 x 8 d + 212M 1 M 2 x 6 d 3 + (d + 95)M 2 2 x 4 4d 2 (d -2)

.

And finally

V 1 n n i=1 SR 3,3 Qi,ρi (f xy ) ≤ 2.66M 2 1 x 8 nd + 212M 1 M 2 x 6 nd 3 + (d + 95)M 2 2 x 4 4nd 2 (d -2)

.

C BUTTERFLY MATRICES GENERATION DETAILS C.1 NOT A POWER OF TWO:

We discuss here the procedure to generate butterfly matrices of size d × d when d is not a power of 2.

Let the number of butterfly factors k = log d . Then B (d) is constructed as a product of k factor matrices of size d × d obtained from k matrices used for generating B (2 k ) . For each matrix in the product for B (2 k ) , we delete the last 2 kd rows and columns. We then replace with 1 every c i in the remaining d × d matrix that is in the same column as deleted s i .

For the cases when d is not a power of two, the resulting B has deficient columns with zeros (Figure 4b , right), which introduces a bias to the integral estimate. To correct for this bias one may apply additional randomization by using a product BP, where P ∈ {0, 1} d×d is a permutation matrix. Even better, use a product of several BP's: B = (BP) 1 (BP) 2 . . . (BP) t . We set t = 3 in the experiments.

C.2 BUTTERFLY RANDOMIZATION:

The key to uniformly random orthogonal butterfly matrix B is the sequence of d-1 angles θ i . To get B (d) Haar distributed, we follow Fang & Li (1997) algorithm that first computes a uniform random point u from U d . It then calculates the angles by taking the ratios of the appropriate u coordinates θ i = ui ui+1 , followed by computing cosines and sines of the θ's.

D KERNEL APPROXIMATION ERRORS ON DIFFERENT DATA SETS:

Here we discuss the datasets that did not appear in the main body of the paper. Table 2 displays the settings for the experiments across the datasets. Figure 5 shows the results for the kernel approxi- Quasi-Monte Carlo integration boasts improved rate of convergence 1 /D compared to 1 / √ D of Monte Carlo, however, empirical results illustrate its performance is poorer than that of orthogonal random features (Felix et al. (2016) ). It also has larger constant factor hidden under O notation and higher complexity. For QMC the weight matrix M is generated as a transformation of quasi-random sequences. We run our experiments with Halton sequences in compliance with the previous work (see Figure 5 with QMC method included).

Although, for the arc-cosine kernels, our methods are the best performing estimators, for the Gaussian kernel the error is not always the lowest one and depends on the dataset, e.g. on the USPS dataset the lowest is Monte Carlo with ROM. However, for the most of the datasets we demonstrate superiority of our approach with this kernel.

We also notice that the dataset with a small amount of features, Powerplant, enjoys Halton and Orthogonal Random Features best, while ROM's convergence stagnates at some point. This could be due the small input feature space with d = 4 and we leave it for the future investigation.

E COMPARISON WITH GAUSSIAN QUADRATURES:

We also included subsampled dense grid method from Dao et al. (2017) into our comparison as it is the only data-independent approach from the paper that is shown to work well. We reimplemented code for the paper to the best of our knowledge since it is not open sourced. We run comparison on all datasets with Gaussian kernel. The Figure 6 illustrates the performance on the LETTER dataset across different expansions. We can see almost coinciding performance of the method (denoted GQ) with the baseline RFF (denoted G). For other datasets the figures are very similar to the case of LETTER, with RFF and GQ methods showing nearly matching relative error of kernel approximation. It should also be noted that explicit mappings produced by Gaussian quadratures do not possess any convenient structure and, thus, cannot boast any better computational complexity.

F WALLTIME EXPERIMENT:

We have also run the time measurement on our somewhat unoptimized implementation of the proposed method. Indeed, Figure 7 demonstrates that the method scales as theoretically predicted with larger dimensions thanks to the structured nature of the mapping. ). Thanks to structured explicit mapping of the proposed method B, it is favorable for higher dimensional data.

:

It may be the case when sampling ρ that 1 -d ρ

< 0, simple solution is just to resample ρ to satisfy the non-negativity of the expression.2 We do not use reflected points for Gaussian kernel as noted earlier, so W j = ρ j σ (Q j V) .

We also study quasi-Monte Carlo (Yang et al. (2014)) performance. See Appendix D for details.

B = (BP)1(BP)2 . . . (BP)3, where P is a permutation matrix, for explanation see Appendix C.

Figures and Tables:

g. D-dimensional approximation admits O(N D 2 ) training time, O(N D) memory and O(N ) prediction time. It is well known that as D → ∞, the inner product in (2) converges to exact kernel k(x, y). Recent research (Yang et al. (2014); Felix et al. (2016); Choromanski & Sindhwani (

d+1) , w i is the row in matrix W. The matrix W = ρ (QV) -(QV)

(Rahimi & Recht (2008)), structured Gaussian (Felix et al. (2016)), random orthogonal matrices (ROM) (Choromanski et al. (2017)).

Figure 2: Kernel approximation error across three kernels (columns: arc-cosine 0, arc-cosine 1, Gaussian) on three datasets: LETTER (d = 16), MNIST (d = 784), CIFAR100 (d = 3072) and LEUKEMIA(d = 7129). Lower is better. The x-axis represents the factor to which we extend the original feature space, n = D 2(d+1) , where d is the dimensionality of the original feature space, D is the dimensionality of the new feature space.

Choromanski et al. (2017) propose a generalization of the ideas from (Le et al. (2013)) and (Felix et al. (2016)), delivering an analytical estimate for the mean squared error (MSE) of approximation. All these works use simple Monte Carlo sampling. However, the convergence can be improved by changing Monte Carlo sampling to Quasi-Monte Carlo sampling. Following this idea Yang et al. (2014) apply quasi-Monte Carlo to Random Fourier Features. In (Yu et al. (

Figure 4: (a) Butterfly orthogonal matrix factors for d = 16. (b) Sparsity pattern for BPBPBP (left) and B (right), d = 15.

Figure 5: Kernel approximation error across three kernels (columns: arc-cosine 0, arc-cosine 1, Gaussian) on three datasets: Powerplant (d = 4), LETTER (d = 16), USPS (d = 256), MNIST (d = 784). Lower is better. The x-axis represents the factor to which we extend the original feature space, n = D 2(d+1) , where d is the dimensionality of the original feature space, D is the dimensionality of the new feature space.

Figure 6: Kernel approximation error for Gaussian kernel on LETTER dataset. Subsampled dense grid method (denoted GQ) from Dao et al. (2017) show very similar performance to the baseline random Fourier features (denoted G on the picture).

Figure 7: Walltime experiment measures the time spent on applying the explicit mapping across different The x-axis represents the 5 datasets with increasing input number of features: LETTER (d = 16), USPS (d = 256), MNIST (d = 784), CIFAR100 (d = 3072) and LEUKEMIA (d = 7129). Thanks to structured explicit mapping of the proposed method B, it is favorable for higher dimensional data.

Space and time complexity for different kernel approximation algorithms.

Table content: <table><tr><td>Method</td><td>Space</td><td>Time</td></tr><tr><td>ORF QMC ROM Quadrature based</td><td>O(Dd) O(Dd) O(d) O(d)</td><td>O(Dd) O(Dd) O(d log d) O(d log d)</td></tr></table>

Experimental settings for the datasets. N is the total number of objects, d is dimensionality of the original feature space.

Table content: <table><tr><td>Dataset</td><td>N</td><td>d</td><td colspan="2">Number of samples Number of runs</td></tr><tr><td>Powerplant</td><td>9568</td><td>4</td><td>550</td><td>500</td></tr><tr><td>LETTER</td><td>20000</td><td>16</td><td>550</td><td>500</td></tr><tr><td>USPS</td><td>9298</td><td>256</td><td>550</td><td>500</td></tr><tr><td>MNIST</td><td colspan="2">70000 784</td><td>50</td><td>50</td></tr><tr><td>CIFAR100</td><td colspan="2">60000 3072</td><td>50</td><td>50</td></tr><tr><td>LEUKEMIA</td><td>72</td><td>7129</td><td>10</td><td>10</td></tr></table>

Peer review (iclr2018):
 The authors offer a novel version of the random feature map approach to approximately solving large-scale kernel problems: each feature map evaluates the "fourier feature" corresponding to the kernel at a set of randomly sampled quadrature points. This gives an unbiased kernel estimator; they prove a bound its variance and provide experiment evidence that for Gaussian and arc-cos kernels, their suggested qaudrature rule outperforms previous random feature maps in terms of kernel approximation error and in terms of downstream classification and regression tasks. The idea is straightforward, the analysis seems correct, and the experiments suggest the method has superior accuracy compared to prior RFMs for shift-invariant kernels. The work is original, but I would say incremental, and the relevant literature is cited. The method seems to give significantly lower kernel approximation errors, but the significance of the performance difference in downstream ML tasks is unclear --- the confidence intervals of the different methods overlap sufficiently to make it questionable whether the relative complexity of this method is worth the effort. Since good performance on downstream tasks is the crucial feature that we want RFMs to have, it is not clear that this method represents a true improvement over the state-of-the-art. The exposition of the quadrature method is difficult to follow, and the connection between the quadrature rules and the random feature map is never explicitly stated: e.g. equation 6 says how the kernel function is approximated as an integral, but does not give the feature map that an ML practitioner should use to get that approximate integral. It would have been a good idea to include figures showing the time-accuracy tradeoff of the various methods, which is more important in large-scale ML applications than the kernel approximation error. It is not clear that the method is *not* more expensive in practice than previous methods (Table 1 gives superior asymptotic runtimes, but I would like to see actual run times, as evaluating the feature maps sound relatively complicated compared to other RFMs). On a related note, I would also like to have seen this method applied to kernels where the probability density in the Bochner integral was not the Gaussian density (e.g., the Laplacian kernel): the authors suggested that their method works there as well when one uses a Gaussian approximation of the density (which is not clear to me), --- and it may be the case that sampling from their quadrature distribution is faster than sampling from the original non-Gaussian density.