Fastfood: Approximate Kernel Expansions in Loglinear Time

Introduction

Kernel methods have proven to be a highly successful technique for solving many problems in machine learning, ranging from classification and regression to sequence annotation and feature extraction (Boser et al., 1992; Cortes and Vapnik, 1995; Vapnik et al., 1997; Taskar et al., 2004; Schölkopf et al., 1998). At their heart lies the idea that inner products in high-dimensional feature spaces can be computed in implicit form via kernel function $k$:

Here $\phi:\mathcal{X}\to\mathcal{F}$ is a feature map transporting elements of the observation space $\mathcal{X}$ into a possibly infinite-dimensional feature space $\mathcal{F}$. This idea was first used by Aizerman et al. (1964) to show nonlinear separation. There exists a rich body of literature on Reproducing Kernel Hilbert Spaces (RKHS) (Aronszajn, 1944; Wahba, 1990; Micchelli, 1986) and one may show that estimators using norms in feature space as penalty are equivalent to estimators using smoothness in an RKHS (Girosi, 1998; Smola et al., 1998a). Furthermore, one may provide a Bayesian interpretation via Gaussian Processes. See e.g. (Williams, 1998; Neal, 1994; MacKay, 2003) for details.

More concretely, to evaluate the decision function $f(x)$ on an example $x$, one typically employs the kernel trick as follows

This has been viewed as a strength of kernel methods, especially in the days that datasets consisted of ten thousands of examples. This is because the Representer Theorem of Kimeldorf and Wahba (1970) states that such a function expansion in terms of finitely many coefficients must exist under fairly benign conditions even whenever the space is infinite dimensional. Hence we can effectively perform optimization in infinite dimensional spaces. This trick that was also exploited by Schölkopf et al. (1998) for evaluating PCA. Frequently the coefficient space is referred to as dual space. This arises from the fact that the coefficients are obtained by solving a dual optimization problem.

Unfortunately, on large amounts of data, this expansion becomes a significant liability for computational efficiency. For instance, Steinwart and Christmann (2008) show that the number of nonzero $\alpha_{i}$ (i.e., $N$, also known as the number of “support vectors”) in many estimation problems can grow linearly in the size of the training set. As a consequence, as the dataset grows, the expense of evaluating $f$ also grows. This property makes kernel methods expensive in many large scale problems: there the sample size $m$ may well exceed billions of instances. The large scale solvers of Fan et al. (2008) and Matsushima et al. (2012) work in primal space to sidestep these problems, albeit at the cost of limiting themselves to linear kernels, a significantly less powerful function class.

Related Work

Numerous methods have been proposed to mitigate this issue. To compare computational cost of these methods we make the following assumptions:

We have $m$ observations and access to an $O(m^{\beta})$ with $\beta\geq 1$ algorithm for solving the optimization problem at hand. In other words, the algorithm is linear or worse. This is a reasonable assumption — almost all data analysis algorithm need to inspect the data at least once to draw inference.

Data has $d$ dimensions. For simplicity we assume that it is dense with density rate $\rho$, i.e. on average $O(\rho d)$ coordinates are nonzero.

The number of nontrivial basis functions is $O(\gamma m)$. This is well motivated by Steinwart and Christmann (2008) and it also follows from the fact that e.g. in regularized risk minimization the subgradient of the loss function determines the value of the associated dual variable.

We denote the number of (nonlinear) basis functions by $n$.

Burges (1996) focused on compressing function expansions after the problem was solved by means of reduced-set expansions. That is, one first solves the full optimization problem at $O(m^{\beta+1}\rho d)$ cost and subsequently one minimizes the discrepancy between the full expansion and an expansion on a subset of basis functions. The exponent of $m^{\beta+1}$ arises from the fact that we need to compute $O(m)$ kernels $O(m^{\beta})$ times. Evaluation of the reduced function set costs at least $O(n\rho d)$ operations per instance and $O(n\rho d)$ storage, since each kernel function $k(x_{i},\cdot)$ requires storage of $x_{i}$.

Low Rank Expansions

Subsequent work by Smola and Schölkopf (2000); Fine and Scheinberg (2001) and Williams and Seeger (2001) aimed to reduce memory footprint and complexity by finding subspaces to expand functions. The key difference is that these algorithms reduce the function space before seeing labels. While this is suboptimal, experimental evidence shows that for well designed kernels the basis functions extracted in this fashion are essentially as good as reduced set expansions. This is to be expected. After all, the kernel encodes our prior belief in which function space is most likely to capture the relevant dependencies between covariates and labels. These projection-based algorithms generate an $n$-dimensional subspace:

Compute the kernel matrix $K_{nn}$ on an $n$-dimensional subspace at $O(n^{2}\rho d)$ cost.

The matrix $K_{nn}$ is inverted at $O(n^{3})$ cost.

For all observations one computes an explicit feature map by projecting data in RKHS onto the set of $n$ basis vectors via $\phi(x)=K_{n}n^{-\frac{1}{2}}\left[k(x_{1},x),\ldots,k(x_{n},x)\right]$. That is, training proceeds at $O(n\rho m^{\beta}+n^{2}m)$ cost.

Prediction costs $O(n\rho d)$ computation and $O(n\rho d)$ memory, as in reduced set methods, albeit with a different set of basis functions.

Note that these methods temporarily require $O(n^{2})$ storage during training, since we need to be able to multiply with the inverse covariance matrix efficiently. This allows for solutions to problems where $m$ is in the order of millions and $n$ is in the order of thousands: for $n=10^{4}$ we need approximately $1$GB of memory to store and invert the covariance matrix. Preprocessing can be parallelized efficiently. Obtaining a minimal set of observations to project on is even more difficult and only the recent work of Das and Kempe (2011) provides usable performance guarantees for it.

Multipole Methods

Fast multipole expansions (Lee and Gray, 2009; Gray and Moore, 2003) offer one avenue for efficient function expansions whenever the dimensionality of the underlying space is relatively modest. However, for high dimensions they become computationally intractable in terms of space partitioning, due to the curse of dimensionality. Moreover, they are typically tuned for localized basis functions, specifically the Gaussian RBF kernel.

Random Subset Kernels

While the paper suggests that the algorithm is scalable to large amounts of data, it suffers from essentially the same problem as other feature generation methods insofar as it needs to evaluate set membership for each of the partitions for all data, hence we have an $O(knm)$ computational cost for $n$ partitions into $k$ sets on $m$ observations. Even this estimate is slightly optimistic since we assume that computing the partitions is independent of the dimensionality of the data. In summary, while the function class is potentially promising, its computational cost considerably exceeds that of the other algorithms discussed below, hence we do not investigate it further.

Random Kitchen Sinks

A promising alternative was proposed by Rahimi and Recht (2009) under the moniker of Random Kitchen Sinks. In contrast to previous work the authors attempt to obtain an explicit function space expansion directly. This works for translation invariant kernel functions by performing the following operations:

Generate a (Gaussian) random matrix $M$ of size $n\times d$.

For each observation $x$ compute $Mx$ and apply a nonlinearity $\psi$ to each coordinate separately, i.e. $\phi_{i}(x)=\psi([Mx]_{i})$.

The approach requires $O(n\times d)$ storage both at training and test time. Training costs $O(m^{\beta}n\rho d)$ operations and prediction on a new observation costs $O(n\rho d)$. This is potentially much cheaper than reduced set kernel expansions. The experiments in (Rahimi and Recht, 2009) showed that performance was very competitive with conventional RBF kernel approaches while providing dramatically simplified code.

Note that explicit spectral finite-rank expansions offer potentially much faster rates of convergence, since the spectrum decays as fast as the eigenvalues of the associated regularization operator Williamson et al. (2001). Nonetheless Random Kitchen Sinks are a very attractive alternative due to their simple construction and the flexility in synthesizing kernels with predefined smoothness properties.

Fastfood

Our approach hews closely to random kitchen sinks. However, it succeeds at overcoming their key obstacle — the need to store and to multiply by a random matrix. This way, fastfood, accelerates Random Kitchen Sinks from $O(nd)$ to $O(n\log d)$ time while only requiring $O(n)$ rather than $O(nd)$ storage. The speedup is most significant for large input dimensions, a common case in many large-scale applications. For instance, a tiny 32x32x3 image in the CIFAR-10 (Krizhevsky, 2009) already has 3072 dimensions, and non-linear function classes have shown to work well for MNIST (Schölkopf and Smola, 2002) and CIFAR-10. Our approach relies on the fact that Hadamard matrices, when combined with Gaussian scaling matrices, behave very much like Gaussian random matrices. That means these two matrices can be used in place of Gaussian matrices in Random Kitchen Sinks and thereby speeding up the computation for a large range of kernel functions. The computational gain is achieved because unlike Gaussian random matrices, Hadamard matrices admit FFT-like multiplication and require no storage.

We prove that the Fastfood approximation is unbiased, has low variance, and concentrates almost at the same rate as Random Kitchen Sinks. Moreover, we extend the range of applications from radial basis functions $k(\left\|x-x^{\prime}\right\|)$ to any kernel that can be written as dot product $k(\left\langle x,x^{\prime}\right\rangle)$. Extensive experiments with a wide range of datasets show that Fastfood achieves similar accuracy to full kernel expansions and Random Kitchen Sinks while being 100x faster with 1000x less memory. These improvements, especially in terms of memory usage, make it possible to use kernel methods even for embedded applications.

Our experiments also demonstrate that Fastfood, thanks to its speedup in training, achieves state-of-the-art accuracy on the CIFAR-10 dataset (Krizhevsky, 2009) among permutation-invariant methods. Table 1 summarizes the computational cost of the above algorithms.

Having an explicit function expansion is extremely beneficial from an optimization point of view. Recent advances in both online (Ratliff et al., 2007) and batch (Teo et al., 2010; Boyd et al., 2010) subgradient algorithms summarily rely on the ability to compute gradients in the feature space $\mathcal{F}$ explicitly.

Kernels and Regularization

For concreteness and to allow for functional-analytic tools we need to introduce some machinery from regularization theory and functional analysis. The derivation is kept brief but we aim to be self-contained. A detailed overview can be found e.g. in the books of Schölkopf and Smola (2002) and Wahba (1990).

When solving a regularized risk minimization problem one needs to choose a penalty on the functions employed. This can be achieved e.g. via a simple norm penalty on the coefficients

Alternatively we could impose a smoothness requirement which emphasizes simple functions over more complex ones via

One may show that the choice of feature map $\phi(x)$ and RKHS norm $\left\|\cdot\right\|_{\mathcal{H}}$ are connected. This is formalized in the reproducing property

In other words, inner products in feature space $\mathcal{F}$ can be viewed as inner products in the RKHS. An immediate consequence of the above is that $k(x,x^{\prime})=\left\langle k(x,\cdot),k(x^{\prime},\cdot)\right\rangle_{\mathcal{H}}$. It also means that whenever norms can be written via regularization operator $P$, we may find $k$ as the Greens function of the operator. That is, whenever $\left\|f\right\|_{\mathcal{H}}^{2}=\left\|Pf\right\|^{2}_{L_{2}}$ we have

That is, $P^{\dagger}Pk(x,\cdot)$ as like a delta distribution on $f\in\mathcal{H}$. This allows us to identify $P^{\dagger}P$ from $k$ and vice versa (Smola et al., 1998a; Girosi, 1998; Girosi et al., 1995; Girosi and Anzellotti, 1993; Wahba, 1990). Note, though, that this need not uniquely identify $P$, a property that we will be taking advantage of when expressing a given kernel in terms of global and local basis functions. For instance, any isometry $U$ with $U^{\top}U=\mathbf{1}$ generates an equivalent $P^{\prime}=UP$. In other words, there need not be a unique feature space representation that generates a given kernel (that said, all such representations are equivalent).

2 Mercer’s Theorem and Feature Spaces

A key tool is the theorem of Mercer (1909) which guarantees that kernels can be expressed as an inner product in some Hilbert space.

The key idea of Rahimi and Recht (2008) is to use sampling to approximate the sum in (6). Note that for trace-class kernels, i.e. for kernels with finite $\sum_{j}\lambda_{j}$ we can normalize the sum to mimic a probability distribution, i.e. we have

Consequently the following approximation converges for $n\to\infty$ to the true kernel

Note that the basic connection between random basis functions was well established, e.g., by Neal (1994) in proving that the Gaussian Process is a limit of an infinite number of basis functions. A related strategy can be found in the so-called ‘empirical’ kernel map (Tsuda et al., 2002; Schölkopf and Smola, 2002) where kernels are computed via

for $x_{i}$ often drawn from the same distribution as the training data. An explicit expression for this map is given e.g. in (Smola et al., 1998b). The expansion (8) is possible whenever the following conditions hold:

An inner product expansion of the form (6) is known for a given kernel $k$.

The basis functions $\phi_{j}$ are sufficiently inexpensive to compute.

The norm $\left\|\lambda\right\|_{1}$ exists, i.e., $k$ corresponds to a trace class operator Kreyszig (1989).

Although condition 2 is typically difficult to achieve, there exist special classes of expansions that are computationally attractive. Specifically, whenever the kernels are invariant under the action of a symmetry group, we can use the eigenfunctions of its representation to diagonalize the kernel.

3 Kernels via Symmetry Groups

Of particular interest in our case are kernels with some form of group invariance since in these cases it is fairly straightforward to identify the basis functions $\phi_{i}(x)$. The reason is that whenever $k(x,x^{\prime})$ is invariant under a symmetry group transformation of its arguments, it means that we can find a matching eigensystem efficiently, simply by appealing to the functions that decompose according to the irreducible representation of the group.

For details see e.g. Berg et al. (1984). This means that knowledge of a group invariance dramatically simplifies the task of finding an eigensystem that satisfies the Mercer decomposition. Moreover, by construction unitary representations are orthonormal.

To make matters more concrete, consider translation invariant kernels

The matching symmetry group is translation group with the Fourier basis admitting a unitary irreducible representation. Corresponding kernels can be expanded

This expansion is particularly simple since the translation group is Abelian. By construction the function $\lambda(z)$ is obtained by applying the Fourier transform to $k(x,0)$ — in this case the above expansion is simply the inverse Fourier transform. We have

This is a well studied problem and for many kernels we may obtain explicit Fourier expansions. For instance, for Gaussians it is a Gaussian with the inverse covariance structure. For the Laplace kernel it yields the damped harmonic oscillator spectrum. That is, good choices of $\lambda$ are

Here the first follows from the fact that Fourier transforms of Gaussians are Gaussians and the second equality follows from the fact that the Fourier spectrum of Bessel functions can be expressed as multiple convolution of the unit sphere. For instance, this includes the Bernstein polynomials as special case for one-dimensional problems. For a detailed discussion of spectral properties for a broad range of kernels see e.g. (Smola, 1998).

Spherical Harmonics

Kernels that are rotation invariant can be written as an expansion of spherical harmonics. (Smola et al., 2001, Theorem 5) shows that dot-product kernels of the form $k(x,x^{\prime})=\kappa(\left\langle x,x^{\prime}\right\rangle)$ can be expanded in terms of spherical harmonics. This provides necessary and sufficient conditions for certain families of kernels. Since Smola et al. (2001) derive an incomplete characterization involving an unspecified radial contribution we give a detailed derivation below.

Here $Y_{n,j}^{d}$ are orthogonal polynomials of degree $n$ on the $d$-dimensional sphere. Moreover, $N(d,n)$ denotes the number of linearly independent homogeneous polynomials of degree $n$ in $d$ dimensions, and $\Omega_{d-1}$ denotes the volume of the $d-1$ dimensional unit ball. $L_{n,d}$ denotes the Legendre polynomial of degree $n$ in $d$ dimensions. Finally, $\lambda_{n}$ denotes the expansion coefficients of $\kappa$ in terms of $L_{n,d}$.

Equality between the two expansions follows from the addition theorem of spherical harmonics of order $n$ in $d$ dimensions. Hence, we only need to show that for $\left\|x\right\|=\left\|x^{\prime}\right\|=1$ the expansion $\kappa(\left\langle x,x^{\prime}\right\rangle)=\sum_{n}\lambda_{n}L_{n,d}(\left\langle x,x^{\prime}\right\rangle)$ holds.

First, observe, that such an expansion is always possible since the Legendre polynomials are orthonormal with respect to the measure induced by the $d-1$ dimensional unit sphere, i.e. with respect to $(1-t^{2})^{\frac{d-3}{2}}$. See e.g. (Hochstadt, 1961, Chapter 3) for details. Hence they form a complete basis for one-dimensional expansions of $\kappa(\xi)$ in terms of $L_{n,d}(\xi)$. Since $\kappa$ is analytic, we can extend the homogeneous polynomials radially by expanding according to (16). This proves the correctness.

To show that this expansion provides necessary and sufficient conditions for positive semidefiniteness, note that $Y_{l,n}^{d}$ are orthogonal polynomials. Hence, if we had $\lambda_{n}<0$ we could use any matching $Y_{l,n}^{d}$ to falsify the conditions of Mercer’s theorem.

Finally, the last equality follows from the fact that $\int_{S_{d-1}}Y_{l,n}^{d}Y_{l^{\prime},n}^{d}=\delta_{l,l^{\prime}}$, i.e. the functions $Y_{l,n}^{d}$ are orthogonal polynomials. Moreover, we use the series expansion of $L_{n,d}$ that also established equality between the first and second line. ∎

The integral representation of (17) may appear to be rather cumbersome. Quite counterintuitively, it holds the key to a computationally efficient expansion for kernels depending on $\left\langle x,x^{\prime}\right\rangle$ only. This is the case since we may sample from a spherically isotropic distribution of unit vectors $z$ and compute Legendre polynomials accordingly. As we will see, computing inner products with spherically isotropic vectors can be accomplished very efficiently using a construction described in Section 4.

Denote by $\lambda_{n}$ the coefficients obtained by a Legendre polynomial series expansion of $\kappa(\left\langle x,x^{\prime}\right\rangle)$ and let $N(d,n)=\frac{(d+n-1)!}{n!(d-1)!}$ be the number of linearly independent homogeneous polynomials of degree $n$ in $d$ variables. Draw $z_{i}\sim S_{d-1}$ uniformly from the unit sphere and draw $n_{i}$ from a spectral distribution with $p(n)\propto\lambda_{n}N(d,n)$. Then

In other words, provided that we are able to compute the Legendre polynomials $L_{n,d}$ efficiently, and provided that it is possible to draw from the spectral distribution of $\lambda_{n}N(d,n)$, we have an efficient means of computing dot-product kernels.

For kernels on the symmetric group that are invariant under group action, i.e. kernels satisfying $k(x,x^{\prime})=k(g\circ x,g\circ x^{\prime})$ for permutations, expansions using Young Tableaux can be found in (Huang et al., 2007). A very detailed discussion of kernels on symmetry groups is given in (Kondor, 2008, Section 4). However, efficient means of computing such kernels rapidly still remains an open problem.

4 Explicit Templates

In some cases expanding into eigenfunctions of a symmetry group may be undesirable. For instance, the Fourier basis is decidedly nonlocal and function expansions using it may exhibit undesirable local deviations, effectively empirical versions of the well-known Gibbs phenomenon. That is, local changes in terms of observations can have far-reaching global effects on observations quite distant from the observed covariates.

This makes it desirable to expand estimates in terms of localized basis functions, such as Gaussians, Epanechikov kernels, B-splines or Bessel functions. It turns out that the latter is just as easily achievable as the more commonplace nonlocal basis function expansions. Likewise, in some cases the eigenfunctions are expensive to compute and it would be desirable to replace them with possibly less statistically efficient alternatives that offer cheap computation.

Consequently we generalize the above derivation to general nonlinear function classes dependent on matrix multiplication or distance computation with respect to spherically symmetric sets of instances. The key is that the feature map depends on $x$ only via

That is, the feature map depends on $x$ and $z$ only in terms of their norms and an inner product between both terms. Here the dominant cost of evaluating $\phi_{z}(x)$ is the inner product $x^{\top}z$. All other operations are $O(1)$, provided that we computed $\left\|x\right\|$ and $\left\|z\right\|$ previously as a one-off operation. Eq. (19) includes the squared distance as a special case:

Here $\kappa$ is suitably normalized, such as $\int dz\kappa(z)=1$. In other words, we expand $x$ in terms of how close the observations are to a set of well-defined anchored basis functions. It is clear that in this case

is a kernel function since it can be expressed as an inner product. Moreover, provided that the basis functions $\phi_{z}(x)$ are well bounded, we can use sampling from the (normalized) measure $\mu(z)$ to obtain an approximate kernel expansion

Note that there is no need to obtain an explicit closed-form expansion in (22). Instead, it suffices to show that this expansion is well-enough approximated by draws from $\mu(z)$.

This is a locally weighted variant of the conventional Gaussian RBF kernel, e.g. as described by Haussler (1999). While this loses its translation invariance, one can easily verify that for $b\to 0$ it converges to the conventional kernel. Note that the key operation in generating an explicit kernel expansion is to evaluate $\left\|z_{i}-x\right\|$ for all $i$. We will explore settings where this can be achieved for $n$ locations $z_{i}$ that are approximately random at only $O(n\log d)$ cost, where $d$ is the dimensionality of the data. Any subsequent scaling operation is $O(n)$, hence negligible in terms of aggregate cost. Finally note that by dividing out the terms related only to $\left\|x\right\|$ and $\left\|x^{\prime}\right\|$ respectively we obtain a ’proper’ Gaussian RBF kernel. That is, we use the following features:

Weighting functions $\mu(z)$ that are more spread-out than a Gaussian will yield basis function expansions that are more adapted to heavy-tailed distributions. It is easy to see that such expansions can be obtained simply by specifying an algorithm to draw $z$ rather than having to express the kernel $k$ in closed form at all.

Polynomial Expansions

One of the main inconveniences in computational evaluation of Corollary 4 is that we need to evaluate the associated Legendre polynomials $L_{n,d}(\xi)$ directly. This is costly since currently there are no known $O(1)$ expansions for the associate Legendre polynomials, although approximate $O(1)$ variants for the regular Legendre polynomials exist (Bogaert et al., 2012). This problem can be alleviated by considering the following form of polynomial kernels:

In this case we only need the ability to draw from the uniform distribution over the unit sphere to compute a kernel. The price to be paid for this is that the effective basis function expansion is rather more complex. To compute it we use the following tools from the theory of special functions.

For odd $a$ the integral vanishes, which follows immediately from the dependence on $x^{a}$ and the symmetric domain of integration $$.

That is, we decompose $x$ into its first coordinate $x_{1}$ and the remainder $x_{2}$ that lies on $S_{d-2}$ with suitable rescaling by $(1-x_{1})^{\frac{d-3}{2}}$. Note the exponent of $\frac{d-3}{2}$ that arises from the curvature of the unit sphere. See e.g. (Hochstadt, 1961, Chapter 6) for details.

While (28) offers a simple expansion for sampling, it is not immediately useful in terms of describing the kernel as a function of $\left\langle x,x^{\prime}\right\rangle$. For this we need to solve the integral in (28). Without loss of generality we may assume that $x=(x_{1},x_{2},0,\ldots 0)$ and that $x^{\prime}=(1,0,\ldots 0)$ with $\left\|x\right\|=\left\|x^{\prime}\right\|=1$. In this case a single summand of (28) becomes

Using the fact that $x_{1}=\theta$ and $x_{2}=\sqrt{1-\theta^{2}}$ we have the full expansion of (28) via

The above form is quite different from commonly used inner-product kernels, such as an inhomogeneous polynomial $\left(\left\langle x,x^{\prime}\right\rangle+d\right)^{p}$. That said, the computational savings are considerable and the expansion bears sufficient resemblance to warrant its use due to significantly faster evaluation.

Sampling Basis Functions

We now discuss computationally efficient strategies for approximating the function expansions introduced in the previous section, beginning with Random Kitchen Sinks of Rahimi and Recht (2008), as described in Section 3.2. Direct use for Gaussian RBF kernels yields the following algorithm to approximate kernel functions by explicit feature construction:

As discussed previously, and as shown by Rahimi and Recht (2009), the associated feature map converges in expectation to the Gaussian RBF kernel. Moreover, they also show that this convergence occurs with high probability and at the rate of independent empirical averages. While this allows one to use primal space methods, the approach remains limited by the fact that we need to store $Z$ and, more importantly, we need to compute $Zx$ for each $x$. That is, each observation costs $O(n\cdot d)$ operations. This seems wasteful, given that we are really only multiplying $x$ with a ‘random’ matrix $Z$, hence it seems implausible to require a high degree of accuracy for $Zx$.

The above idea can be improved to extend matters beyond a Gaussian RBF kernel and to reduce the memory footprint in computationally expensive settings. We summarize this in the following two remarks:

To avoid storing the Gaussian random matrix $Z$ we recompute $Z_{ij}$ on the fly. Assume that we have access to a random number generator which takes samples from the uniform distribution $\xi\sim U$ as input and emits samples from a Gaussian, e.g. by using the inverse cumulative distribution function $z=F^{-1}(\xi)$. Then we may replace the random number generator by a hash function via $\xi_{ij}=N^{-1}h(i,j)$ where $N$ denotes the range of the hash, and subsequently $Z_{ij}=F^{-1}(\xi_{ij})$.

Unfortunately this variant is computationally even more costly than Random Kitchen Sinks, its only benefit being the $O(n)$ memory footprint relative to the $O(nd)$ footprint for random kitchen sinks. To make progress, a more effective approximation of the Gaussian random matrix $Z$ is needed.

2 Fastfood

The coefficients for $S,G,B$ are computed once and stored. On the other hand, the Walsh-Hadamard matrix is never computed explicitly. Instead we only multiply by it via the fast Hadamard transform, a variant of the FFT which allows us to compute $H_{d}x$ in $O(d\log d)$ time. The Hadamard matrices are defined as follows:

When $n>d,$ we replicate (33) for $n/d$ independent random matrices $V_{i}$ and stack them via $V^{T}=[V_{1},V_{2},\ldots V_{n/d}]^{T}$ until we have enough dimensions. The feature map for Fastfood is then defined as

Before proving that in expectation this transform yields a Gaussian random matrix, let us briefly verify the computational efficiency of the method.

The features of (34) can be computed at $O(n\log d)$ cost using $O(n)$ permanent storage for $n\geq d$.

Storing the matrices $S,G,B$ costs $3n$ entries and $3n$ operations for a multiplication. The permutation matrix $\Pi$ costs $n$ entries and $n$ operations. The Hadamard matrix itself requires no storage since it is only implicitly represented. Furthermore, the fast Hadamard transforms costs $O(n\log d)$ operations to carry out since we have $O(d\log d)$ per block and $n/d$ blocks. Computing the Fourier basis for $n$ numbers is an $O(n)$ operation. Hence the total CPU budget is $O(n\log d)$ and the storage is $O(n)$. ∎

Note that the construction of $V$ is analogous to that of Dasgupta et al. (2011). We will use these results in establishing a sufficiently high degree of decorrelation between rows of $V$. Also note that multiplying with a longer chain of Walsh-Hadamard matrices and permutations would yield a distribution closer to independent Gaussians. However, as we shall see, two matrices provide a sufficient amount of decorrelation.

3 Basic Properties

Now that we showed that the above operation is fast, let us give some initial indication why it is also useful and how the remaining matrices $S,G,B,\Pi$ are defined.

This is a diagonal matrix with $B_{ii}\in\left\{\pm 1\right\}$ drawn iid from the uniform distribution over $\left\{\pm 1\right\}$. The initial $HBd^{-\frac{1}{2}}$ acts as an isometry that densifies the input, as pioneered by Ailon and Chazelle (2009).

It ensures that the rows of the two Walsh-Hadamard matrices are incoherent relative to each other. $\Pi$ can be stored efficiently as a lookup table at $O(d)$ cost and it can be generated by sorting random numbers.

This is a diagonal matrix whose elements $G_{ii}\sim\mathcal{N}(0,1)$ are drawn iid from a Gaussian. The next Walsh-Hadamard matrices $H$ will allow us to ’recycle’ $n$ Gaussians to make the resulting matrix closer to an iid Gaussian. The goal of the preconditioning steps above is to guarantee that no single $G_{ii}$ can influence the output too much and hence provide near-independence.

Note that the length of all rows of $HG\Pi HB$ are constant as equation (36) shows below. In the Gaussian case $S$ ensures that the length distribution of the row of $V$ are independent of each other. In the more general case, one may also adjust the capacity of the function class via a suitably chosen scaling matrix $S$. That is, large values in $S_{ii}$ correspond to high complexity basis functions whereas small $S_{ii}$ relate to simple functions with low total variation. For the RBF kernel we choose

Thus $s_{i}$ matches the radial part of a normal distribution and we rescale it using the Frobenius norm of $G$.

We now analyze the distribution of entries in $V$.

Rescaling the length of a Gaussian vector using (35) retains Gaussianity. Hence the rows of $SHG\Pi HB$ are Gaussian, albeit not independent.

The expected feature map recovers the Gaussian RBF kernel, i.e.,

Moreover, the same holds for $V^{\prime}=\frac{1}{\sigma\sqrt{d}}HG\Pi HB$.

We already proved that any given row in $V$ is a random Gaussian vector with distribution $\mathcal{N}(0,\sigma^{-2}I_{d})$, hence we can directly appeal to the construction of Rahimi and Recht (2008). This also holds for $V^{\prime}$. The main difference being that the rows in $V^{\prime}$ are considerably more correlated. Note that by assembling several $d\times d$ blocks to obtain an $n\times d$ matrix this property is retained, since each block is drawn independently. ∎

4 Changing the Spectrum

Changing the kernel from a Gaussian RBF to any other radial basis function kernel is straightforward. After all, $HG\Pi HB$ provides a approximately spherically uniformly distributed random vectors of the same length. Rescaling each direction of projection separately costs only $O(n)$ space and computation. Consequently we are free to choose different coefficients $S_{ii}$ rather than (35). Instead, we may use

Here $c$ is a normalization constant and $\lambda(r)$ is the radial part of the spectral density function of the regularization operator associated with the kernel.

A key advantage over a conventional kernel approach is that we are not constrained by the requirement that the spectral distributions be analytically computable. Even better, we only need to be able to sample from the distribution rather than compute its Fourier integral in closed form.

While this may appear costly, it only needs to be carried out once at initialization time and it allows us to sidestep computing the convolution entirely. After that we can store the coefficients $S_{ii}$. Also note that this addresses a rather surprising problem with the Gaussian RBF kernel — in high dimensional spaces draws from a Gaussian are strongly concentrated on the surface of a sphere. That is, we only probe the data with a fixed characteristic length. The Matern kernel, on the other hand, spreads its capacity over a much larger range of frequencies.

5 Inner Product Kernels

We now put Theorem 3 and Corollary 4 to good use. Recall that the latter states that any dot-product kernel can be obtained by taking expectations over draws from the degree of corresponding Legendre polynomial and over a random direction of reference, as established by the integral representation of (17).

It is understood that the challenging part is to draw vectors uniformly from the unit sphere. Note, though, that it is this very operation that Fastfood addresses by generating pseudo-Gaussian vectors. Hence the modified algorithm works as follows:

Note that the equality $L_{n_{i},d}(t_{i})=r^{n_{i}}L_{n_{i},d}(t_{i}/r)$ follows from the fact that $L_{n_{i},d}$ is a homogeneous polynomial of degree $n_{i}$. The second representation may sometimes be more effective for reasons of numerical stability. As can be seen, this relies on access to efficient Legendre polynomial computation. Recent work of Bogaert et al. (2012) shows that (quite surprisingly) this is possible in $O(1)$ time for $L_{n}(t)$ regardless of the degree of the polynomial. Extending these guarantees to associated Legendre polynomials is unfortunately rather nontrivial. Hence, a direct expansion in terms of $\left\langle x,v\right\rangle^{d}$, as discussed previously, or brute force computation may well be more effective.

We conclude our reasoning by providing an extension of the above argument to the symmetric group. Clearly, by treating permutation matrices $\Pi\in C_{n}$ as $d\times d$ dimensional vectors, we can use them as inputs to a dot-product kernel. Subsequently, taking inner products with random reference vectors of unit length yields kernels which are dependent on the matching between permutations only.

Analysis

The next step is to show that the feature map is well behaved also in terms of decorrelation between rows of $V$. We focus on Gaussian RBF kernels in this context.

When using random kitchen sinks, the variance of the feature map is at least $O(1/n)$ since we draw $n$ samples iid from the space of parameters. In the following we show that the variance of fastfood is comparable, i.e. it is also $O(1/n)$, albeit with a dependence on the magnitude of the magnitude of the inputs of the feature map. This guarantee matches empirical evidence that both algorithms perform equally well as the exact kernel expansion.

For convenience, since the kernel values are real numbers, let us simplify terms and rewrite the inner product in terms of a sum of cosines. Trigonometric reformulation yields

Let $v=\frac{x-x^{\prime}}{\sigma}$ and let $\psi_{j}(v)=\cos[Vv]_{j}$ denote the estimate of the kernel value arising from the $j$th pair of random features for each $j\in\left\{1\ldots d\right\}$. Then for each $j$ we have

where $C(\alpha)=6\alpha^{4}\left[e^{-\alpha^{2}}+\frac{\alpha^{2}}{3}\right]$ depends on the scale of the argument of the kernel.

We decompose $Vv$ into a sequence of terms $w=d^{-\frac{1}{2}}HBv$ and $u=\Pi w$ and $z=HGu$. Hence we have $\psi_{j}(v)=\cos(z_{j})$. Note that $\left\|u\right\|=\left\|v\right\|$ since by construction $d^{-\frac{1}{2}}H$, $B$ and $\Pi$ are orthonormal matrices.

Observe that the marginal distribution of each $z_{j}$ is $\mathcal{N}(0,\left\|v\right\|^{2})$ since each element of $H$ is $\pm 1$. Thus the joint distribution of $z_{j}$ and $z_{t}$ is a Gaussian with mean and covariance

where $h\sim\mathcal{N}(0,1)$ and $a^{2}_{\pm}=\left\|L^{\top}[1,\pm 1]\right\|^{2}=2\left\|v\right\|^{2}(1\pm\rho)$. That is, after applying the addition theorem we explicitly computed the now one-dimensional Gaussian integrals.

We compute the first moment analogously. Since by construction $z_{j}$ and $z_{j}$ have zero mean and variance $\left\|v\right\|^{2}$ we have that

Combining both terms we obtain that the covariance can be written as

Taylor Expansion

To prove the first claim note that here $j=t$, since we are computing the variance of a single feature. Correspondingly $\rho(u)=1$. Plugging this into (41) and simplifying terms yields the first claim of (39).

To prove our second claim, observe that from the Taylor series of $\cosh$ with remainder in Lagrange form, it follows that there exists $\eta\in[-\left\|v\right\|^{2}|\rho|,\left\|v\right\|^{2}|\rho|]$ such

where $B(\left\|v\right\|)=\frac{1}{2}+\frac{\sinh(\left\|v\right\|^{2})\left\|v\right\|^{2}}{6}$. Plugging this into (41) yields

Note that the above is still conditioned on $u$. What remains is to bound $\mathbf{E}_{u}[\rho^{2}]$, which is small if $\mathbf{E}[\left\|u\right\|_{4}^{4}]$ is small. The latter is ensured by $HB$, which acts as a randomized preconditioner: Since $G$ is diagonal and $G_{ii}\sim\mathcal{N}(0,1)$ independently we have

Now recall that $u=\Pi w$ and that $\Pi$ is a random permutation matrix. Therefore $u_{i}=w_{\pi(i)}$ for a randomly chosen permutation $\pi$ and thus the distribution of $\rho\left\|v\right\|^{2}$ and $2\sum_{i\in R}w_{i}^{2}-\left\|v\right\|^{2}$ where $R$ is a randomly chosen subset of size $\frac{d}{2}$ in $\left\{1\ldots d\right\}$ are the same. Let us fix (condition on) $w$. Since $2\mathbf{E}_{R}\left[\sum_{i\in R}w_{i}^{2}\right]=\left\|v\right\|^{2}$ we have that

Now let $\delta_{i}=1$ if $i\in R$ and otherwise. Note that $\mathbf{E}_{\delta}(\delta_{i})=\frac{1}{2}$ and if $j\neq k$ then $\mathbf{E}_{\delta}(\delta_{i}\delta_{k})\leq\frac{1}{4}$ as $\delta_{i}$ are (mildly) negatively correlated. From $\left\|w\right\|=\left\|v\right\|$ it follows that

From the two equations above it follows that

Bounding the fourth moment of ‖w‖norm𝑤\left\|w\right\|

Let $b_{i}=B_{ii}$ be the independent $\pm 1$ random variables of $B$. Using the fact $w_{i}=\frac{1}{\sqrt{d}}\sum_{t=1}^{d}H_{it}b_{t}v_{t}$ and that $b_{i}$ are independent with similar calculations to the above it follows that

which shows that $\frac{1}{\sqrt{d}}HB$ acts as preconditioner that densifies the input. Putting it all together we have

Combining the latter with the already proven first claim establishes the second claim. ∎

Denote by $V,V^{\prime}$ Gauss-like matrices of the form

Moreover, let $C(\alpha)=6\alpha^{4}\left[e^{-\alpha^{2}}+\frac{\alpha^{2}}{3}\right]$ be a scaling function. Then for the feature maps obtained by stacking $n/d$ iid copies of either $V$ or $V^{\prime}$ we have

Since $\phi^{\prime}(x)^{\top}\phi^{\prime}(x^{\prime})$ is the average of $n/d$ independent estimates, each arising from $2d$ features. Hence we can appeal to Theorem 9 for a single block, i.e. when $n=d$. The near-identical argument for $V$ is omitted. ∎

2 Concentration

The following theorem shows that for a given error probability $\delta$, the approximation error of a $d\times d$ block of Fastfood is at most logarithmically larger than the error of Random Kitchen Sinks. That is, it is only logarithmically weaker. We believe that this bound is pessimistic and could be further improved with considerable analytic effort. That said, the $O(m^{-\frac{1}{2}})$ approximation guarantees to the kernel matrix are likely rather conservative when it comes to generalization performance, as we found in experiments. In other words, we found that the algorithm works much better in practice than in theory, as confirmed in Section 6. Nonetheless it is important to establish tail bounds, not to the least since this way improved guarantees for random kitchen sinks also immediately benefit fastfood.

Our key tool is concentration of Lipschitz continuous functions under the Gaussian measure Ledoux (1996). We ensure that Fastfood construct has a small Lipschitz constant using the following lemma, which is due to Ailon and Chazelle (2009).

In other words, with high probability, the largest elements of $d^{-\frac{1}{2}}HBx$ are with high probability no larger than what one could expect if all terms were of even size as per the $\left\|x\right\|_{2}$ norm.

[Theorem 11] Since both $k$ and $\hat{k}$ are shift invariant we set $v=\sigma(x-x^{\prime})$ and write $k(v)=k(x,x^{\prime})$ and $\hat{k}(v)=\hat{k}(x,x^{\prime})$ to simplify the notation. Set $u=\Pi d^{-\frac{1}{2}}HBv$, and $z=HGu$ and define

From $u=\Pi d^{-\frac{1}{2}}HBv$ and $\left\|\Pi w\right\|_{\infty}=\left\|w\right\|_{\infty}$ combined with Lemma 12 it follows that

holds with probability at least $1-\delta$, where the probability is over the choice of $B$.Note that in contrast to Theorem 9, the permutation matrix $\Pi$ does not play a role in the proof of Theorem 11. Now condition on (50). From inequality (49) we have that the function

is Lipschitz continuous with Lipschitz constant

Hence from Theorem 13 and from the independently chosen $G_{jj}\sim\mathcal{N}(0,1)$ it follows that

Combining inequalities (51) and (52) with the union bound concludes the proof. ∎

Experiments

In the following we assess the performance of Random Kitchen Sinks and Fastfood. The results show that Fastfood performs as well as Random Kitchen Sinks in terms of accuracy. Fastfood, however, is orders of magnitude faster and exhibits a significantly lower memory footprint. For simplicity, we focus on penalized least squares regression since in this case we are able to compute exact solutions and are independent of any other optimization algorithms. We also benchmark Fastfood on CIFAR-10 (Krizhevsky, 2009) and observe that it achieves state-of-the-art accuracy. This advocates for the use of non-linear expansions even when $d$ is large.

We begin by investigating how well our features can approximate the exact kernel computation as $n$ increases. For that purpose, we uniformly sample 4000 vectors from $^{10}$. We compare the exact kernel values to Random Kitchen Sinks and Fastfood.

The results are shown in Figure 1. We used the absolute difference between the exact kernel and the approximation to quantify the error (the relative difference also exhibits similar behavior and is thus not shown due to space constraints). The results are presented as averages, averaging over 4000 samples. As can be seen, as $n$ increases, both Random Kitchen Sinks and Fastfood converge quickly to the exact kernel values. Their performance is indistinguishable, as expected from the construction of the algorithm.

Note, though, that fidelity in approximating $k(x,x^{\prime})$ does not imply generalization performance (unless the bounds are very tight). To assess this, we carried out experiments on all regression datasets from the UCI repository (Frank and Asuncion, 2010) that are not too tiny, i.e., that contained at least $4,000$ instances.

We investigate estimation accuracy via Gaussian process regression (Rasmussen and Williams, 2006) using approximated kernel computation methods and we compare this to exact kernel computation whenever the latter is feasible. For completeness, we compare the following methods:

uses the exact Gaussian RBF kernel, that is $k(x,x^{\prime})=\exp\left(-\left\|x-x^{\prime}\right\|^{2}/2\sigma^{2}\right)$. This is possible, albeit not practically desirable due to its excessive cost, on all but the largest datasets where the kernel matrix does not fit into memory.

uses the Nystrom approximation of the kernel matrix (Williams and Seeger, 2001). These methods have received recent interest due to the improved approximation guarantees of Jin et al. (2011) which indicate that approximation rates faster than $O(n^{-\frac{1}{2}})$ are achievable. Hence, theoretically, the Nystrom method could have a significant accuracy advantage over Random Kitchen Sinks and Fastfood when using the same number of basis functions, albeit at exponentially higher cost of $O(d)$ vs. $O(\log d)$ per function. We set $n=2,048$ to retain a computationally feasible feature projection.

uses the the Gaussian random projection matrices of Rahimi and Recht (2008). As before, we use $n=2,048$ basis functions. Note that this is a rather difficult setting for Random Kitchen Sinks relative to the Nystrom decomposition, since the basis functions obtained in the latter are arguably better in terms of approximating the kernel. Hence, one would naively expect slightly inferior performance from Random Kitchen Sinks relative to direct Hilbert Space methods.

(Hadamard features) uses the random matrix given by $SHG\Pi HB$, again with $n=2,048$ dimensions. Based on the above reasoning one would expect that the performance of the Hadamard features is even weaker than that of Random Kitchen Sinks since now the basis functions are no longer even independently drawn from each other.

(Fourier features) uses a variant of the above construction. Instead of combining two Hadamard matrices, a permutation and Gaussian scaling, we use a permutation in conjunction with a Fourier Transform matrix $F$: the random matrix given by $V=\Pi FB$. The motivation is the Subsampled Random Fourier Transform, as described by Tropp (2010): by picking a random subset of columns from a (unitary) Fourier matrix, we end up with vectors that are almost spatially isotropic, albeit with slightly more dispersed lengths than in Fastfood. We use this heuristic for comparison purposes.

uses the exact polynomial kernel, that is $k(x,x^{\prime})=(\left\langle z,x\right\rangle+1)^{d}$, with $d=10$. Similar to the case of Exact RBF, this method is only practical on small datasets.

uses the Fastfood trick via Spherical Harmonics to approximate the polynomial kernels.

The results of the comparison are given in Table 3. As can be seen, and contrary to the intuition above, there is virtually no difference between the exact kernel, the Nystrom approximation, Random Kitchen Sinks and Fastfood. In other words, Fastfood performs just as well as the exact method, while being substantially cheaper to compute. Somewhat surprisingly, the Fourier features work very well. This indicates that the concentration of measure effects impacting Gaussian RBF kernels may actually be counterproductive at their extreme. This is corroborated by the good performance observed with the Matern kernel.

In Figure 2, we show regression performance as a function of the number of basis functions $n$ on the CPU dataset. As is evident, it is necessary to have a large $n$ in order to learn highly nonlinear functions. Interestingly, although the Fourier features do not seem to approximate the Gaussian RBF kernel, they perform well compared to other variants and improve as $n$ increases. This suggests that learning the kernel by direct spectral adjustment might be a useful application of our proposed method.

2 Speed of kernel computations

In the previous experiments, we observe that Fastfood is on par with exact kernel computation, the Nystrom method, and Random Kitchen Sinks. The key point, however, is to establish whether the algorithm offers computational savings.

For this purpose we compare Random Kitchen Sinks using Eigenhttp://eigen.tuxfamily.org/index.php?title=Main_Page and our method using Spiralhttp://spiral.net. Both are highly optimized numerical linear algebra libraries in C++. We are interested in the time it takes to go from raw features of a vector with dimension $d$ to the label prediction of that vector. On a small problem with $d=1,024$ and $n=16,384$, performing prediction with Random Kitchen Sinks takes 0.07 seconds. Our method is around 24x faster, taking only 0.003 seconds to compute the label for one input vector. The speed gain is even more significant for larger problems, as is evident in Table 2. This confirms experimentally the $O(n\log d)$ vs. $O(nd)$ runtime and the $O(n)$ vs. $O(nd)$ storage of Fastfood relative to Random Kitchen Sinks. In other words, the computational savings are substantial for large input dimensionality $d$.

3 Random features for CIFAR-10

To understand the importance of nonlinear feature expansions for a practical application, we benchmarked Fastfood, Random Kitchen Sinks on the CIFAR-10 dataset Krizhevsky (2009) which has 50,000 training images and 10,000 test images. Each image has 32x32 pixels and 3 channels ($d=3072$). In our experiments, linear SVMs achieve 42.3% accuracy on the test set. Non-linear expansions improve the classification accuracy significantly. In particular, Fastfood FFT (“Fourier features”) achieve 63.1% while Fastfood (“Hadamard features”) and Random Kitchen Sinks achieve 62.4% with an expansion of $n=16,384$. These are also best known classification accuracies using permutation-invariant representations on this dataset. In terms of speed, Random Kitchen Sinks is 5x slower (in total training time) and 20x slower (in predicting a label given an image) compared to both Fastfood and and Fastfood FFT. This demonstrates that non-linear expansions are needed even when the raw data is high-dimensional, and that Fastfood is more practical for such problems.

In particular, in many cases, linear function classes are used because they provide fast training time, and especially test time, but not because they offer better accuracy. The results on CIFAR-10 demonstrate that Fastfood can overcome this obstacle.

Summary

We demonstrated that it is possible to compute $n$ nonlinear basis functions in $O(n\log d)$ time, a significant speedup over the best competitive algorithms. This means that kernel methods become more practical for problems that have large datasets and/or require real-time prediction. In fact, Fastfood can be used to run on cellphones because not only it is fast, but it also requires only a small amount of storage.

Note that our analysis is not limited to translation invariant kernels but it also includes inner product formulations. This means that for most practical kernels our tools offer an easy means of making kernel methods scalable beyond simple subspace decomposition strategies. Extending our work to other symmetry groups is subject to future research. Also note that fast multiplications with near-Gaussian matrices are a key building block of many randomized algorithms. It remains to be seen whether one could use the proposed methods as a substitute and reap significant computational savings.

References