Title: NOT-SO-RANDOM FEATURES

Abstract: We propose a principled method for kernel learning, which relies on a Fourieranalytic characterization of translation-invariant or rotation-invariant kernels. Our method produces a sequence of feature maps, iteratively refining the SVM margin. We provide rigorous guarantees for optimality and generalization, interpreting our algorithm as online equilibrium-finding dynamics in a certain two-player min-max game. Evaluations on synthetic and real-world datasets demonstrate scalability and consistent improvements over related random features-based methods.

1. INTRODUCTION:

Choosing the right kernel is a classic question that has riddled machine learning practitioners and theorists alike. Conventional wisdom instructs the user to select a kernel which captures the structure and geometric invariances in the data. Efforts to formulate this principle have inspired vibrant areas of study, going by names from feature selection to multiple kernel learning (MKL).

We present a new, principled approach for selecting a translation-invariant or rotation-invariant kernel to maximize the SVM classification margin. We first describe a kernel-alignment subroutine, which finds a peak in the Fourier transform of an adversarially chosen data-dependent measure. Then, we define an iterative procedure that produces a sequence of feature maps, progressively improving the margin. The resulting algorithm is strikingly simple and scalable.

Intriguingly, our analysis interprets the main algorithm as no-regret learning dynamics in a zero-sum min-max game, whose value is the classification margin. Thus, we are able to quantify convergence guarantees towards the largest margin realizable by a kernel with the assumed invariance. Finally, we exhibit experiments on synthetic and benchmark datasets, demonstrating consistent improvements over related random features-based kernel methods.

1.1. RELATED WORK:

There is a vast literature on MKL, from which we use the key concept of kernel alignment (Cristianini et al., 2002) . Otherwise, our work bears few similarities to traditional MKL; this and much related work (e.g. Cortes et al. (2012) ; Gönen & Alpaydın (2011) ; Lanckriet et al. (2004) ) are concerned with selecting a kernel by combining a collection of base kernels, chosen beforehand. Our method allows for greater expressivity and even better generalization guarantees.

Instead, we take inspiration from the method of random features (Rahimi & Recht, 2007) . In this pioneering work, originally motivated by scalability, feature maps are sampled according to the Fourier transform of a chosen kernel. The idea of optimizing a kernel in random feature space was studied by Sinha & Duchi (2016) . In this work, which is most similar to ours, kernel alignment is optimized via importance sampling on a fixed, finitely supported proposal measure. However, the proposal can fail to contain informative features, especially in high dimension; indeed, they highlight efficiency, rather than showing performance improvements over RBF features.

Learning a kernel in the Fourier domain (without the primal feature maps) has also been considered previously: Oliva et al. (2016) and Yang et al. (2015b) model the Fourier spectrum parametrically, which limits expressivity; the former also require complicated posterior inference procedures. Bȃzȃvan et al. (2012) study learning a kernel in the Fourier domain jointly with regression parameters. They show experimentally that this locates informative frequencies in the data, without theoretical guarantees. Our visualizations suggest this approach can get stuck in poor local minima, even in 2 dimensions. Crammer et al. (2003) also use boosting to build a kernel sequentially; however, they only consider a basis of linear feature maps, and require costly generalized eigenvector computations. From a statistical view, Fukumizu et al. (2009) bound the SVM margin in terms of maximum mean discrepancy, which is equivalent to (unweighted) kernel alignment. Notably, their bound can be loose if the number of support vectors is small; in such situations, our theory provides a tighter characterization. Moreover, our attention to the margin goes beyond the usual objective of kernel alignment.

1.2. OUR CONTRIBUTIONS:

We present an algorithm that outputs a sequence of Fourier features, converging to the maximum realizable SVM classification margin on a labeled dataset. At each iteration, a pair of features is produced, which maximizes kernel alignment with a changing, adversarially chosen measure. As this measure changes slowly, the algorithm builds a diverse and informative feature representation.

Our main theorem can be seen as a case of von Neumann's min-max theorem for a zero-sum concave-linear game; indeed, our method bears a deep connection to boosting (Freund & Schapire, 1996; Schapire, 1999) . In particular, both the theory and empirical evidence suggest that the generalization error of our method decreases as the number of random features increases. In traditional MKL methods, generalization bounds worsen as base kernels are added.

Other methods in the framework of Fourier random features take the approach of approximating a kernel by sampling feature maps from a continuous distribution. In contrast, our method constructs a measure with small finite support, and realizes the kernel exactly by enumerating the associated finite-dimensional feature map; there is no randomness in the features.1 

2. PRELIMINARIES:



2.1. THE FOURIER DUAL MEASURE OF A KERNEL:

We focus on two natural families of kernels k(x, x ): translation-invariant kernels on X = R d , which depend only on xx , and rotation-invariant kernels on the hypersphere X = S d-1 , which depend only on x, x . These invariance assumptions subsume most widely-used classes of kernels; notably, the Gaussian (RBF) kernel satisfies both. For the former invariance, Bochner's theorem provides a Fourier-analytic characterization: Theorem 2.1 (e.g. Eq. (1), Sec. 1.4.3 in Rudin (2011)) . A translation-invariant continuous function k : R d × R d → R is positive semidefinite if and only if k(xx ) is the Fourier transform of a symmetric non-negative measure (where the Fourier domain is

EQUATION

for some λ ∈ L 1 (Ω), satisfying λ(ω) ≥ 0 and λ(ω) = λ(-ω) for all ω ∈ Ω.

A similar characterization is available for rotation-invariant kernels, where the Fourier basis functions are the spherical harmonics, a countably infinite family of complex polynomials which form an orthonormal basis for square-integrable functions X → C. 

EQUATION

for some λ ∈ L 1 (Ω), with λ(ω) ≥ 0 and λ(ω) = λ(-ω) for all valid index pairs ω = ( , m) ∈ Ω.

In each of these cases, we call this Fourier transform λ k (ω) the dual measure of k. This measure decomposes k into a non-negative combination of Fourier basis kernels. Furthermore, this decomposition gives us a feature map φ : X → L 2 (Ω, λ k ) whose image realizes the kernel under the codomain's inner product;2 that is, for all x, x ∈ X ,

EQUATION

Respectively, these feature maps are φ x (ω) = e i ω,x and φ x ( , m) = Y d ,m (x). Although they are complex-valued, symmetry of λ k allows us to apply the transformation {φ x (ω), φ x (-ω)} → {Re φ x (ω), Im φ x (ω)} to yield real features, preserving the inner product. The analogous result holds for spherical harmonics.

2.2. KERNEL ALIGNMENT:

In a binary classification task with n training samples (x i ∈ X , y i ∈ {±1}), a widely-used quantity for measuring the quality of a kernel k : X × X → R is its alignment (Cristianini et al., 2002; Cortes et al., 2012) ,3 defined by

γ k (P, Q) def = i,j∈[n] k(x i , x j )y i y j = y T G k y,

where y is the vector of labels and G k is the Gram matrix. Here, we let P = i:yi=1 δ xi and Q = i:yi=-1 δ xi denote the (unnormalized) empirical measures of each class, where δ x is the Dirac measure at x. When P, Q are arbitrary measures on X , this definition generalizes to

γ k (P, Q) def = X 2 k(x, x ) dP 2 (x, x ) + X 2 k(x, x ) dQ 2 (x, x ) - X 2 k(x, x ) dP(x) dQ(x ).

In terms of the dual measure λ k (ω), kernel alignment takes a useful alternate form, noted by Sriperumbudur et al. (2010) . Let µ denote the signed measure P -Q. Then, when k is translationinvariant, we have

EQUATION

Analogously, when k is rotation-invariant, we have

EQUATION

It can also be verified that λ k L1(Ω) = k(x, x), which is of course the same for all x ∈ X . In each case, the alignment is linear in λ k . We call v(ω) the Fourier potential, which is the squared magnitude of the Fourier transform of the signed measure µ = P -Q. This function is clearly bounded pointwise by (P(X ) + Q(X )) 2 .

3. ALGORITHMS:



3.1. MAXIMIZING KERNEL ALIGNMENT IN THE FOURIER DOMAIN:

First, we consider the problem of finding a kernel k (subject to either invariance) that maximizes alignment γ k (P, Q); we optimize the dual measure λ k (ω). Aside from the non-negativity and symmetry constraints from Theorems 2.1 and 2.2, we constrain λ k 1 = 1, as this quantity appears as a normalization constant in our generalization bounds (see Theorem 4.3). Maximizing γ k (P, Q) in this constraint set, which we call L ⊂ L 1 (Ω), takes the form of a linear program on an infinitedimensional simplex. Noting that v(ω) = v(-ω) ≥ 0, γ k is maximized by placing a Dirac mass at any pair of opposite modes ±ω * ∈ argmax ω v(ω).

At first, P and Q will be the empirical distributions of the classes, specified in Section 2.2. However, as Algorithm 2 proceeds, it will reweight each data point x i in the measures by α(i). Explicitly, the reweighted Fourier potential takes the form4 

v α (ω) def = n i=1 y i α i e ι ω,xi 2 or v α ( , m) def = n i=1 y i α i Y d ,m (x i ) 2 .

Due to its non-convexity, maximizing v α (ω), which can be interpreted as finding a global Fourier peak in the data, is theoretically challenging. However, we find that it is easy to find such peaks in our experiments, even in hundreds of dimensions. This arises from the empirical phenomenon that realistic data tend to be band-limited, a cornerstone hypothesis in data compression. An 2 constraint (or equivalently, 2 regularization; see Kakade et al. ( 2009)) can be explicitly enforced to promote band-limitedness; we find that this is not necessary in practice.

When a gradient is available (in the translation-invariant case), we use Langevin dynamics (Algorithm 1) to find the peaks of v(ω), which enjoys mild theoretical hitting-time guarantees (see Theorem 4.5). See Appendix A.3 for a discussion of the (discrete) rotation-invariant case.

Algorithm 1 Langevin dynamics for kernel alignment

1: Input: training samples S = {(x i , y i )} n i=1 , weights α ∈ R n . 2: Parameters: time horizon τ , diffusion rate ζ, temperature ξ. 3: Initialize ω 0 arbitrarily. 4: for t = 0, . . . , τ -1 do 5: Update ω t+1 ← ω t + ζ∇v α (ω t ) + 2ξ

ζ z, where z ∼ N (0, I). 6: end for 7: return ω * := argmax t v(ω t )

It is useful in practice to use parallel initialization, running m concurrent copies of the diffusion process and returning the best single ω encountered. This admits a very efficient GPU implementation: the multi-point evaluations v α (ω 1..m ) and ∇v α (ω 1..m ) can be computed from an (m, d) by (d, n) matrix product and pointwise trigonometric functions. We find that Algorithm 1 typically finds a reasonable peak within ∼100 steps.

3.2. LEARNING THE MARGIN-MAXIMIZING KERNEL FOR SVM:

Support vector machines (SVMs) are perhaps the most ubiquitous use case of kernels in practice. To this end, we propose a method that boosts Algorithm 1, building a kernel that maximizes the classification margin. Let k be a kernel with dual λ k , and Y def = diag(y). Write the dual l 1 -SVM objective,5 parameterizing the kernel by λ k :

F (α, λ k ) def = 1 T α - 1 2 α T YG k Yα = 1 T α - 1 2 γ k (αP, αQ) = 1 T α - 1 2 Ω λ k (ω) v α (ω) dω.

Thus, for a fixed α, F is equivalent to kernel alignment, and can be minimized by Algorithm 1. However, the support vector weights α are of course not fixed; given a kernel k, α is chosen to maximize F , giving the (reciprocal) SVM margin. In all, to find a kernel k which maximizes the margin under an adversarial choice of α, one must consider a two-player zero-sum game:

EQUATION

)

where L is the same constraint set as in Section 3.1, and K is the usual dual feasible set {0 α C, y T α = 0} with box parameter C.

In this view, we make some key observations. First, Algorithm 1 allows the min-player to play a pure-strategy best response to the max-player. Furthermore, a mixed strategy λ for the min-player is simply a translation-or rotation-invariant kernel, realized by the feature map corresponding to its support. Finally, since the objective is linear in λ and concave in α, there exists a Nash equilibrium (λ * , α * ) for this game, from which λ * gives the margin-maximizing kernel.

We can use no-regret learning dynamics to approximate this equilibrium. Algorithm 2 runs Algorithm 1 for the min-player, and online gradient ascent (Zinkevich, 2003) for the max-player. Intuitively (and as is visualized in our synthetic experiments), this process slowly morphs the landscape of v α to emphasize the margin, causing Algorithm 1 to find progressively more informative features.

At the end, we simply concatenate these features; contingent on the success of the kernel alignment steps, we have approximated the Nash equilibrium.

Algorithm 2 No-regret learning dynamics for SVM margin maximization

1: Input: training samples S = {(x i , y i )} n i=1 . 2: Parameters: box constraint C, # steps T , step sizes {η t }; parameters for Algorithm 1. 3: Set φ x (ω) = e i ω,x (translation-invariant), or φ x ( , m) = Y d ,m (x) (rotation-invariant). 4: Initialize α = Proj K [ C 2 • 1]. 5: for t = 1, . . . , T do 6:

Use Algorithm 1 (or other v α maximizer) on S with weights α t , returning ω t .

7::

Append two features {Re φ

x (ω t ), Im φ x (-ω t )} to each x i 's representation Φ(x i ). 8: Compute gradient g t := ∇ α F (α t , λ t ), where λ t = δ ωt + δ -ωt . 9: Update α t+1 ← Proj K [α t + η t g t ]. 10: end for 11: return features {Φ(x i ) ∈ R 2T } n i=1 , or dual measure λ := 1 2T T t=1 δ ωt + δ -ωt .

We provide a theoretical analysis in Section 4.1, and detailed discussion on heuristics, hyperparameters, and implementation details in depth in Appendix A. One important note is that the online gradient g t = 1-2Y Re( Φ t , Yα t Φ t ) is computed very efficiently, where Φ t (i) = φ xi (ω t ) ∈ C m is the vector of the most recently appended features. Langevin dynamics, easily implemented on a GPU, comprise the primary time bottleneck.

4. THEORY:



4.1. CONVERGENCE OF NO-REGRET LEARNING DYNAMICS:

We first state the main theoretical result, which quantifies the convergence properties of Algorithm 2. Theorem 4.1 (Main). Assume that at each step t, Algorithm 1 returns an ε t -approximate global maximizer ω t (i.e., v αt (ω t ) ≥ sup ω∈Ω v αt (ω)ε t ). Then, with a certain choice of step sizes η t , Algorithm 2 produces a dual measure λ ∈ L which satisfies

max α∈K F (α, λ) ≤ min λ∈L max α∈K F (α, λ) + T t=1 ε t 2T + 3(C + 2C 2 )n 4 √ T .

(Alternate form.) Suppose instead that v αt (ω t ) ≥ ρ at each time t. Then, λ satisfies

max α∈K F (α, λ) ≤ ρ + 3(C + 2C 2 )n 4 √ T .

We prove Theorem 4.1 in Appendix C. For convenience, we state the first version as an explicit margin bound (in terms of competitive ratio M/M * ): Corollary 4.2. Let M be the margin obtained by training an 1 linear SVM with the same C as in Algorithm 2, on the transformed samples {(Φ(x i ), y i )}. Then, M is (1δ)-competitive with M * , the maximally achievable margin by a kernel with the assumed invariance, with

δ ≤ t ε t T + 3(C + 2C 2 )n 4 √ T .

This bound arises from the regret analysis of online gradient ascent (Zinkevich, 2003) ; our analysis is similar to the approach of Freund & Schapire (1996) , where they present a boosting perspective.

When using an 2 -SVM, the final term can be improved to O( log T T ) (Hazan et al., 2007) . For a general overview of results in the field, refer to Hazan (2016) .

4.2. GENERALIZATION BOUNDS:

Finally, we state two (rather distinct) generalization guarantees. Both depend mildly on a bandwidth assumption ω 2 ≤ R ω and the norm of the data R x def = max i x i . First, we state a margindependent SVM generalization bound, due to Koltchinskii & Panchenko (2002) . Notice the appearance of λ k 1 = R λ , justifying our choice of normalization constant for Algorithm 1. Intriguingly, the end-to-end generalization error of our method decreases with an increasing number of random features, since the margin bound is being refined during Algorithm 2.

Theorem 4.3 (Generalization via margin). For any SVM decision function f : X → R with a kernel k λ constrained by λ 1 ≤ R λ trained on samples S drawn i.i.d. from distribution D, the generalization error is bounded by

Pr (x,y)∼D [yf (x) ≤ 0] ≤ min θ 1 n n i=1 1 yif (xi)≤θ + 6R ω R x θ R λ n + 3 log 2 δ n .

The proof can be found in Appendix D. Note that this improves on the generic result for MKL, from Theorem 2 in (Cortes et al., 2010) , which has a √ log T dependence on the number of base kernels T . This improvement stems from the rank-one property of each component kernel.

Next, we address another concern entirely: the sample size required for v(ω) to approximate the ideal Fourier potential v ideal (ω), the squared magnitude of the Fourier transform of the signed measure P -Q arising from the true distribution. For the shift-invariant case: . We have that for all ω : ω ≤ R ω , with probability at least

1 -δ, v(ω) n 2 -v ideal (ω) ≤ ε.

The O(•) suppresses factors polynomial in R x and R ω .

The full statement and proof are standard, and deferred to Appendix E. In particular, this result allows for a mild guarantee of polynomial hitting-time on the locus of approximate local maxima of v ideal (as opposed to the empirical v). Adapting the main result from Zhang et al. ( 2017):

Theorem 4.5 (Langevin hitting time). Let ω τ be the output of Algorithm 1 on v α (ω), after τ steps. Algorithm 1 finds an approximate local maximum of v ideal in polynomial time. That is, with U being the set of ε-approximate local maxima of v ideal (ω), some ω t satisfies

v(ω t ) ≥ inf d(U,ω)≤∆ v(ω) for some τ ≤ O(poly (R x , R ω , d, ζ, ξ, ε -1 , ∆ -1 , log(1/δ))

, with probability at least 1δ.

Of course, one should not expect polynomial hitting time on approximate global maxima; Roberts & Tweedie (1996) give asymptotic mixing guarantees.

5. EXPERIMENTS:

In this section, we highlight the most important and illustrative parts of our experimental results.

For further details, we provide an extended addendum to the experimental section in Appendix A. The code can be found at github.com/yz-ignescent/Not-So-Random-Features. First, we exhibit two simple binary classification tasks, one in R 2 and the other on S 2 , to demonstrate the power of our kernel selection method. As depicted in Figure 1 , we create datasets with sharp boundaries, which are difficult for the standard RBF and arccosine (Cho & Saul, 2009) kernel. In both cases, n train = 2000 and n test = 50000.6 Used on a 1 -SVM classifier, these baseline kernels saturate at 92.1% and 95.1% test accuracy, respectively; they are not expressive enough.

On the R 2 "windmill" task, Algorithm 2 chooses random features that progressively refine the decision boundary at the margin. By T = 1000, it exhibits almost perfect classification (99.7% training, 99.3% test). Similarly, on the S 2 "checkerboard" task, Algorithm 2 (with some adaptations described Published as a conference paper at ICLR 2018 in Appendix A.3) reaches almost perfect classification (99.7% training, 99.1% test) at T = 100, supported on only 29 spherical harmonics as features.

We provide some illuminating visualizations. Figures 1a and 1b show the evolution of the dual weights, random features, and classifier. As the theory suggests, the objective evolves to assign higher weight to points near the margin, and successive features improve the classifier's decisiveness in challenging regions (1a, bottom). Figure 1c visualizes some features from the S 2 experiment.

Next, we evaluate our kernel on standard benchmark binary classification tasks. Challenging label pairs are chosen from the MNIST (LeCun et al., 1998) and CIFAR-10 (Krizhevsky, 2009) datasets; each task consists of ∼10000 training and ∼2000 test examples; this is considered to be large-scale for kernel methods. Following the standard protocol from Yu et al. ( 2016), 512-dimensional HoG features (Dalal & Triggs, 2005) are used for the CIFAR-10 tasks instead of raw images. Of course, our intent is not to show state-of-the-art results on these tasks, on which deep neural networks easily dominate. Instead, the aim is to demonstrate viability as a scalable, principled kernel method.

We compare our results to baseline random features-based kernel machines: the standard RBF random features (RBF-RF for short), and the method of Sinha & Duchi (2016) (LKRF),7 using the same 1 -SVM throughout. As shown by Table 1 , our method reliably outperforms these baselines, most significantly in the regime of few features. Intuitively, this lines up with the expectation that isotropic random sampling becomes exponentially unlikely to hit a good peak in high dimension; our method searches for these peaks. Furthermore, our theory predicts that the margin keeps improving, regardless of the dimensionality of the feature maps. Indeed, as our classifier saturates on the training data, test accuracy continues increasing, without overfitting. This decoupling of generalization from model complexity is characteristic of boosting methods.

In practice, our method is robust with respect to hyperparameter settings. As well, to outperform both RBF-RF and LKRF with 5000 features, our method only needs ∼ 100 features. Our GPU implementation reaches this point in ∼30 seconds. See Appendix A.1 for more tuning guidelines.

6. CONCLUSION:

We have presented an efficient kernel learning method that uses tools from Fourier analysis and online learning to optimize over two natural infinite families of kernels. With this method, we show meaningful improvements on benchmark tasks, compared to related random features-based methods. Many theoretical questions remain, such as accelerating the search for Fourier peaks (e.g. Hassanieh et al. (2012) ; Kapralov (2016) ). These, in addition to applying our learned kernels to state-of-the-art methods (e.g. convolutional kernel networks (Mairal et al., 2014; Yang et al., 2015a; Mairal, 2016) ), prove to be exciting directions for future work.

A.3 SAMPLING SPHERICAL HARMONICS:

As we note in Section 3.1, it is unclear how to define gradient Langevin dynamics on v(l, m) in the inner-product case, since no topology is available on the indices (l, m) of the spherical harmonics.

One option is to emulate Langevin dynamics, by constructing a discrete Markov chain which mixes to λ(l, m) ∝ e βv (l,m) .

However, we find in our experiments that it suffices to compute λ(l, m) by examining all values of v(l, m) with j no more than some threshold J. One should view this as approximating the kerneltarget alignment objective function via Fourier truncation. This is highly parallelizable: it involves approximately N (m, d) degree-J polynomial evaluations on the same sample data, which can be expressed using matrix multiplication. In our experiments, it sufficed to examine the first 1000 coefficients; we remark that it is unnatural in any real-world datasets to expect that v(ω) only has large values outside the threshold J.

Under Fourier truncation, the domain of λ becomes a finite-dimensional simplex. In the gametheoretic view of 4.1, an approximate Nash equilibrium becomes concretely achievable via Nesterov's excessive gap technique (Nesterov, 2005; Daskalakis et al., 2015) , given that the kernel player's actions are restricted to a mixed strategy over a finite set of basis kernels.

Finally, we note a significant advantage to this setting, where we have a discrete set of Fourier coefficients: the same feature might be found multiple times. When a duplicate feature is found, it need not be concatenated to the representation; instead, the existing feature is scaled appropriately. This accounts for the drastically smaller support of features required to achieve near-perfect classification accuracy.

B MORE ON SPHERICAL HARMONICS:

In this section, we go into more detail about the spherical harmonics in d dimensions. Although all of this material is standard in harmonic analysis, we provide this section for convenience, isolating only the relevant facts.

B.1 SPHERICAL HARMONIC EXPANSION OF A ROTATION-INVARIANT KERNEL:

First, we provide a proof sketch for Theorem 2.2. We rely on the following theorem, an analogue of Bochner's theorem on the sphere, which characterizes rotation-invariant kernels: Theorem B.1 (Schoenberg (1942)). A continuous rotation-invariant function k(x, x ) = k( x, x ) on S d-1 × S d-1 is positive semi-definite if and only if its expansion into Gegenbauer polynomials P d i has only non-negative coefficients, i.e.

EQUATION

with λ m ≥ 0, ∀m ∈ N + .

The Gegenbauer polynomials P d m : [-1, 1] → R are a generalization of the Legendre and Chebyshev polynomials; they form an orthogonal basis on [-1, 1] under a certain weighting function (see Stein & Weiss (2016) ; Gallier (2009) ; Müller (2012) for details). Note that we adopt a different notation from Schoenberg (1942) , which denotes our P d m by P (d-1)/2 m . Unlike Theorem 2.1, Theorem B.1 alone does not provide a clear path for constructing a feature map for inner-product kernel, because the inputs x, x of the kernel are still coupled after the expansion into P d m , which acts on the inner product x, x instead of on x and x separately. Fortunately, (see, e.g., Proposition 1.18 in Gallier (2009) ), the Gegenbauer polynomials evaluated on an inner product admits a decomposition into spherical harmonics Y d ,m :

EQUATION

where |S d-1 | denotes the surface area of S d-1 , and N (d, ) = d-1+ -d-1+ -2 . Here, we use the normalization convention that S d-1 |Y d ,m (x)| 2 dx = 1. From the existence of this expansion follows the claim from Theorem 2.2. For a detailed reference, see Müller (2012) .

B.2 PAIRING THE SPHERICAL HARMONICS:

We now specify the involution ω → -ω on indices of spherical harmonics, which gives a pairing that takes the role of opposite Fourier coefficients in the X = R d case. In particular, the Fourier transform λ of a real-valued function k : R n → R satisfies λ(ω) = λ(-ω), so that the dual measure can be constrained to be symmetric. Now, consider the X = S d-1 case, where the Fourier coefficients are on the set Ω of valid indices of spherical harmonics. We would like a permutation on the indices σ so that σ(σ(ω)) = ω, and λ(ω) = λ(σ(ω)) whenever λ is the spherical harmonic expansion of a real kernel.

For a fixed dimension d and degree , the spherical harmonics form an N := N (d, )-dimensional orthogonal basis for a vector space V of complex polynomials from S d-1 to C, namely thehomogeneous polynomials (with domain extended to R n ) which are harmonic (eigenfunctions of the Laplace-Beltrami operator ∆ S d-1 ). This basis is not unique; an arbitrary choice of basis may not contain pairs of conjugate functions.

Fortunately, such a basis exists constructively, by separating the Laplace-Beltrami operator over spherical coordinates (θ 1 , . . . , θ d-1 ). Concretely, for a fixed d and , the d-dimensional spherical harmonic, indexed by

| 1 | ≤ 2 ≤ 3 ≤ . . . ≤ d-1 , is defined as: Y 1,..., N def = 1 √ 2π e i 1θ1 d-1 k=2 k P d-1 k (θ k ),

where the functions in the product come from a certain family of associated Legendre functions in sin θ. For a detailed treatment, we adopt the construction and conventions from Higuchi (1987) .

In the scope of this paper, the relevant fact is that the desired involution is well-defined in all dimensions. Namely, consider the permutation σ that sends a spherical harmonic indexed by ω =

(d, , 1 , 2 , . . . , d-1 ) to σ(ω) = (d, , -1 , 2 , . . . , d-1 ). Then, Y -ω (x) def = Y σ(ω) (x) = Y ω (x) for all ω, x.

The symmetry condition in Theorem 2.2 follows straightforwardly. By orthonormality, we know that every square-integrable function f : R n → C has a unique decomposition into spherical harmonics, with coefficients λ ω = f, Y ω , so that λ -ω = f, Y ω = λ ω . When f is real-valued, we conclude that λ -ω = -λ ω , as claimed.

C PROOF OF THE MAIN THEOREM:

In this section, we prove the main theorem, which quantifies convergence of Algorithm 2 to the Nash equilibrium. We restate it here: Theorem 4.1. Assume that during each timestep t, the call to Algorithm 1 returns an ε t -approximate global maximizer ω t (i.e. vαt (ω t ) ≥ max ω∈K vαt (ω)ε t ). Then, Algorithm 2 returns a dual measure λ, which satisfies

max α∈K F (α, λ) ≤ min λ∈L max α∈K F (α, λ) + T t=1 ε t 2T + O 1 √ T .

Alternatively with the assumption that at each timestep t, vαt (ω t ) ≥ ρ, λ, satisfies

max α∈K F (α, λ) ≤ ρ + O 1 √ T .

If Algorithm 2 is used on a l 2 -SVM, the regret bound can be improved to be O( log T T ).

Proof. We will make use of the regret bound of online gradient ascent (see, e.g., (Hazan, 2016) ).

Here we only prove the theorem in the case for l 1 -SVM with box constraint C, under the assumption of ε t -approximate optimality of Algorithm 1. Extending the proof to other cases is straightfoward.

Lemma C.1 (Regret bound for online gradient ascent). Let D be the diameter of the constraint set K, and G a Lipschitz constant for an arbitrary sequence of concave functions f t (α) def = F (α, λ t ) on K. Online gradient ascent on f t (α), with step size schedule η t = G D √ t , guarantees the following for all T ≥ 1:

1

T T t=1 f t (α t ) ≥ max α∈K 1 T T t=1 f t (α) - 3GD 2 √ T .

Here, D ≤ C √ n by the box constraint, and we have

G ≤ sup α∈K, ω∈R d 1 -2Y Re( Φ t , Yα t Φ t ) 2 ≤ (1 + 2C) √ n. Thus, our regret bound is 3n(C+2C 2 ) 2 √ T ; that is, for all T ≥ 1, 1 T T t=1 F (α t , δ ωt + δ -ωt ) ≥ max α∈K 1 T T t=1 F (α, δ ωt + δ -ωt )) - 3n(C + 2C 2 ) 2 √ T .

Since at each timestep t, vαt (ω t ) ≥ max ω∈K vαt (ω)ε t , we have by assumption

F (α t , δ ωt + δ -ωt ) ≤ min ω∈K F (α t , δ ω + δ -ω ) + ε t ,

and by concavity of F in α, we know that, for any fixed sequence α 1 , α 2 , . . . , α T , min Proof. For a fixed sample S = {x i } n i=1 , we compute the empirical Rademacher complexity of the composite hypothesis class. Below, σ ∈ R n is a vector of i.i.d. Rademacher random variables. 

RS:

= 1 n E σ     sup λ 1 ≤R λ α G k λ α≤1 n i=1 n j=1 σ i α j k λ (x i , x j )     = 1 n E σ sup λ 1≤Rλ sup α G k λ α≤1 σ G k λ α = 1 n E σ sup λ 1≤Rλ σ G k λ σ = 1 n E σ sup λ 1 ≤R λ σ G k λ σ = 1 n E σ    sup

≤ 2R ω R x √ n ,

:

We note a superficial resemblance to quasi-Monte Carlo methods(Avron et al.,

2016); however, these are concerned with accelerating convergence rates rather than learning a kernel.

In fact, such a pair (φ, λ) exists for any psd kernel k; seeDai et al. (2014) orDevinatz (1953).

Definitions vary up to constants and normalization, depending on the use case.

Whenever i is an index, we will denote the imaginary unit by ι.

Of course, our method applies to l2 SVMs, and has even stronger theoretical guarantees; see Section 4.1.

ntest is chosen so large to measure the true generalization error.

Traditional MKL methods are not tested here, as they are noticeably (> 100 times) slower.

Figures and Tables:

Theorem 4.4 (Generalization of the potential). Let S be a set of i.i.d. training samples on D, with n ≥ O

Evolution of Algorithm 2. Top: The windmill dataset, weighted by dual variables αt. Middle: 2t random features (magenta; last added ±ωt in cyan), overlaying the Fourier potential vα t (ω) in the background. Bottom: Decision boundary of a hinge-loss classifier trained on the kernel, showing refinement at the margin. Contours indicate the value of the margin. (b) Detailed evolution of the Fourier potential vα(ω). Left: Initial potential v(ω), with uniform weights. Right: Reweighted vα(ω) at t = 300. Note the larger support of peaks.(c) Toy dataset on the sphere, with top-3 selected features: spherical harmonics Y0,11, Y4,9, Y8,27.

Figure 1: Visualizations and experiments on synthetic data.

(α t , δ ωt + δ -ωt ) -(α, δ ωt + δ -ωt ) -(α, δ ωt + δ -ωt ) -(α, δ ωt + δ -ωt ) = F (α, λ) by linearity of F in the dual measure; here, λ = 1 2T T t=1 δ ωt +δ -ωt is the approximate Nash equilibrium found by the no-regret learning dynamics.D PROOF OF THEOREM 4.3In this section, we compute the Rademacher complexity of the composition of the learned kernel and the classifier, proving Theorem 4.3.

For any SVM decision function f : X → R with a kernel k λ constrained by λ 1 ≤ R λ trained on samples S drawn i.i.d. from distribution D, the generalization error is bounded by Pr

( ω, x i ) .Note that each term in the formula is the empirical Rademacher complexity of the composition of a 1-Lipschitz function (cosine or sine) with a linear function with 2 norm bounded by R ω . By the composition property of Rademacher complexity, we have the following: i cos( ω, x i )

Comparison on binary classification tasks. Our method is compared against standard RBF random features (RBF-RF), as well as the method from Sinha & Duchi (2016) (LKRF). Left: Performance is measured on the CIFAR-10 automobile vs. truck task, varying the number of features m. Standard deviations over 10 trials are shown, demonstrating high stability.

Table content: <table><tr><td/><td>0.95 1.00</td><td/><td/><td>Dataset</td><td>Method</td><td>100</td><td>Number of features m 500 1000 2000</td><td>5000</td></tr><tr><td>Training accuracy</td><td>0.80 0.85 0.90</td><td/><td/><td>MNIST (1-7) MNIST (4-9)</td><td colspan="2">RBF-RF 97.42 99.02 99.10 99.42 99.31 LKRF 97.60 98.95 99.05 99.21 99.37 ours 99.05 99.28 99.45 99.63 99.65 RBF-RF 87.72 94.11 96.95 97.36 98.08 LKRF 88.42 93.99 95.88 97.18 97.76 ours 93.02 96.31 97.68 98.72 99.01</td></tr><tr><td/><td>0.75 0.70 0.90</td><td/><td>RBF-RF LKRF ours</td><td>MNIST (5-6) CIFAR-10 (plane-bird)</td><td colspan="2">RBF-RF 91.96 97.63 97.85 98.51 98.82 LKRF 92.78 97.60 97.79 98.36 98.54 ours 96.62 98.67 99.18 99.32 99.64 RBF-RF 71.75 78.30 79.85 80.15 81.25 LKRF 73.30 78.80 79.85 80.80 81.55 ours 81.80 82.15 83.00 84.25 85.15</td></tr><tr><td>Test accuracy</td><td>0.75 0.80 0.85</td><td/><td/><td>CIFAR-10 (auto-truck) CIFAR-10 (dog-frog)</td><td colspan="2">RBF-RF 69.70 75.75 82.55 83.35 85.10 LKRF 71.89 77.33 83.62 84.60 84.95 ours 84.80 85.90 86.00 88.25 88.80 RBF-RF 70.80 79.95 81.85 82.50 82.95 LKRF 70.75 78.95 81.25 81.35 83.00 ours 81.95 84.05 85.40 85.95 86.25</td></tr><tr><td/><td>0.65 0.70</td><td>0</td><td>1000 Number of features m 2000 3000 4000 RBF-RF 5000 LKRF ours</td><td>CIFAR-10 (auto-ship) CIFAR-10 (auto-deer)</td><td colspan="2">RBF-RF 70.85 78.00 79.15 81.80 82.55 LKRF 70.55 77.90 79.55 80.95 82.85 ours 81.00 83.70 84.40 86.70 87.35 RBF-RF 84.65 90.00 91.85 92.45 92.65 LKRF 84.85 89.75 91.05 92.65 92.90 ours 92.40 93.60 93.55 94.05 94.90</td></tr></table>

Peer review (iclr2018):
 In this paper the authors consider learning directly Fourier representations of shift/translation invariant kernels for machine learning applications. They choose the alignment of the kernel to data as the objective function to optimize. They empirically verify that the features they learned lead to good quality SVM classifiers. My problem with that paper is that even though at first glance learning adaptive feature maps seems to be an attractive approach, authors* contribution is actually very little. Below I list some of the key problems. First of all the authors claim in the introduction that their algorithm is very fast and with provable theoretical guarantees. But in fact later they admit that the problem of optimizing the alignment is a non-convex problem and the authors end up with a couple of heuristics to deal with it. They do not really provide any substantial theoretical justification why these heuristics work in practice even though they observe it empirically. The assumptions that large Fourier peaks happen close to origin is probably well-justified from the empirical point of view, but it is a hack, not a well established well-grounded theoretical method (the authors claim that in their experiments they found it easy to find informative peaks, even in hundreds of dimensions, but these experiments are limited to the SVM setting, I have no idea how these empirical findings would translate to other kernelized algorithms using these adaptive features). The Langevin dynamics algorithm used by the authors to find the peaks (where the gradient is available) gives only weak theoretical guarantees (as the authors actually admit) and this is a well known method, certainly not a novelty of that paper. Finally, the authors notice that "In the rotation-invariant case, where Ω is a discrete set, heuristics are available". That is really not very informative (the authors refer to the Appendix so I carefully read that part of the Appendix, but it is extremely vague, it is not clear at all how the Langevin dynamics can be "emulated" by a discrete Markov chain that mixes fast; the authors do not provide any justification of that approach, what is the mixing time ?; how the "good emulation property" is exactly measured ?). In the conclusions the authors admit that: "Many theoretical questions remain, such as accelerating the search for Fourier peaks". I think that the problem of accelerating this approach is a critical point that this publication is missing. Without this, it is actually really hard to talk about general mechanism of learning adaptive Fourier features for kernel algorithms (which is how the authors present their contribution); instead we have a method heavily customized and well-tailored to the (not particularly exciting) SVM scenario (with optimization performed by the standard annealing method; it is not clear at all whether for other downstream kernel applications this approach for optimizing the alignment would provide good quality models) that uses lots of task specific hacks and heuristics to efficiently optimize the alignment. Another problem is that it is not clear at all to me how authors* approach can be extended to non shift-invariant kernels that do not benefit from Bochner*s Theorem. Such kernels are very related to neural networks (for instance PNG kernels with linear rectifier nonlinearities correspond to random layers in NNs with ReLU) and in the NN context are much more interesting that radial basis function or in general shift-invariant kernels. A general kernel method should address this issue (the authors just claim in the conclusions that it would be interesting to explore the NN context in more detail). To sum it up, it is a solid submission, but in my opinion without a substantial contribution and working only in a very limited setting when it is heavily relying on many unproven hacks and heuristics.