Random Features for Large-Scale Kernel Machines

Support Vector Machines and other kernel methods are powerful because they can learn any decision boundary, given enough data. The kernel trick lets you work in a potentially infinite-dimensional feature space without ever computing the coordinates — you only need the inner product between every pair of points, packaged in a kernel matrix $K$ where $K_{ij} = k(x_i, x_j)$.

The catch: that matrix has $N^2$ entries. Storing it for $N = 500{,}000$ takes about 1 TB of RAM; computing it takes $O(N^2 d)$ time. Decomposition solvers like SMO or SVMlight nibble at this matrix iteratively, but they still don't scale past a few hundred thousand points. Meanwhile, linear SVMs — which simply find a hyperplane $w^\top x + b = 0$ — run in $O(Nd)$ time, but they can only draw straight decision boundaries.

The question this paper asks: can we get kernel-quality nonlinearity at linear-method speed?

The kernel trick says $k(x, y) = \langle \phi(x), \phi(y) \rangle$ for some (often infinite-dimensional) lifting $\phi$. Rahimi and Recht flip this: instead of working implicitly in the huge $\phi$-space, build a short explicit map $z: \mathbb{R}^d \to \mathbb{R}^D$ (with $D \ll N$) so that:

$k(x, y) \approx z(x)^\top z(y)$

Once every training point is transformed to $z(x_i)$, you store an $N \times D$ matrix instead of $N \times N$. Any linear learner — ridge regression, linear SVM, logistic regression — now solves the problem in $O(ND^2)$ time rather than $O(N^2 d)$.

The key insight is that for shift-invariant kernels — kernels that depend only on the difference $x - y$, such as the Gaussian RBF — Bochner's theorem from harmonic analysis tells us exactly how to build $z$.

A shift-invariant kernel depends only on the displacement: $k(x, y) = k(x - y)$. Bochner's theorem (1933) says that any such continuous positive-definite kernel is the Fourier transform of a non-negative measure. If $k$ is properly scaled so that $k(0) = 1$, its Fourier transform $p(\omega)$ is a proper probability distribution.

This means we can write the kernel as an expectation:

$k(x - y) = \mathbb{E}_{\omega \sim p}\left[ e^{j\omega^\top(x-y)} \right]$

Since both $p(\omega)$ and $k$ are real, the imaginary parts cancel, and we can replace the complex exponential with a cosine. Drawing $\omega$ from $p$ and $b$ uniformly from $[0, 2\pi]$, the scalar function $z_\omega(x) = \sqrt{2}\cos(\omega^\top x + b)$ satisfies:

$\mathbb{E}[z_\omega(x)\, z_\omega(y)] = k(x - y)$

Each random cosine projection is an unbiased estimator of the kernel evaluation. Stack $D$ of them into a vector, normalize by $\sqrt{D}$, and the inner product $z(x)^\top z(y)$ becomes a low-variance approximation. Hoeffding's inequality guarantees exponentially fast convergence in $D$.

Think of it as a frequency scanner: each random cosine projects the input onto a randomly chosen line in the data space, then wraps that projection onto a circle. Two points that land near each other on many circles must have been close in the original kernel sense — the same way two radio stations that produce similar readings on many randomly tuned receivers must be broadcasting similar signals.

Bochner's theorem tells us that every shift-invariant kernel has a spectral distribution $p(\omega)$. The choice of kernel determines where the random frequencies $\omega$ come from:

• Gaussian RBF $k(\Delta) = \exp(-\|\Delta\|^2 / 2)$: draw $\omega$ from a Gaussian $\mathcal{N}(0, I)$.
• Laplacian $k(\Delta) = \exp(-\|\Delta\|_1)$: draw each $\omega_d$ from a Cauchy distribution $\frac{1}{\pi(1 + \omega_d^2)}$.
• Cauchy kernel $k(\Delta) = \prod_d \frac{2}{1 + \Delta_d^2}$: draw $\omega$ from a Laplacian (double-exponential) $e^{-\|\omega\|_1}$.

The pattern is elegant: kernels and their spectral distributions are Fourier duals. A narrow kernel (local influence) has a wide spectral distribution, and vice versa — exactly the uncertainty principle from signal processing.

Not all kernels are smooth. The Laplacian kernel $e^{-|x-y|}$ and other L1-based kernels have sharp corners. For these, Rahimi and Recht propose a second map: Random Binning Features.

The idea is beautifully spatial: partition the input space with a randomly shifted grid at a randomly chosen resolution. Each point is encoded as a binary indicator vector marking which bin it lands in. If two points fall in the same bin, their dot product is 1; otherwise 0. The grid resolution is sampled from a distribution derived from the kernel, so the probability of co-binning equals the kernel value: $\Pr[\hat{x} = \hat{y}] = k(x, y)$.

Stack $P$ independent random binnings and average — you get a low-variance, unbiased estimator of the kernel. Unlike Fourier features, binning is piecewise constant — it explicitly preserves locality, which makes it superior for non-smooth decision boundaries like the Forest Cover dataset where tens of thousands of support vectors are needed.

The two feature maps have complementary strengths:

Random Fourier Features produce smooth, continuous maps well-suited for interpolation tasks — problems where the decision boundary is smooth and the kernel is differentiable (e.g., Gaussian RBF). They excel on the CPU and Census regression benchmarks.

Random Binning Features produce piecewise-constant maps that explicitly preserve locality. They shine on memorization tasks — problems that need many support vectors and have rough, complex boundaries. The Forest Cover classification dataset is the dramatic example: binning achieves 2.2% error (matching exact SVM), while Fourier features struggle at 11.6%.

A practical rule of thumb: if the exact SVM needs few support vectors, Fourier features will work well. If it needs many, prefer binning. And since the two can be mixed freely, you can concatenate both types of features for the best of both worlds.

A natural concern: how good is this approximation? The paper provides uniform convergence bounds — not just for a single pair $(x, y)$, but for all pairs simultaneously over a compact set $M$.

For Random Fourier Features: with high probability, the worst-case approximation error across all pairs is bounded by $\epsilon$ when the number of features is $D = \Omega\!\left(\frac{d}{\epsilon^2} \log \frac{\sigma_p \cdot \text{diam}(M)}{\epsilon} \right)$, where $\sigma_p^2$ is the second moment of the spectral distribution (the trace of the kernel's Hessian at the origin).

The key takeaway: $D$ depends on the data dimension $d$ and the desired accuracy $\epsilon$, not on the number of training points $N$. Even for a million training examples, a few hundred random features can achieve sub-percent approximation error on standard kernels.

The paper evaluates Random Features + ridge regression against Core Vector Machines (CVM) and exact SVMs on five large-scale benchmarks. The results are striking:

• On CPU (6,500 instances, regression): Fourier features achieve 3.6% error in 20 seconds — exact SVM takes 31 seconds for 11% error.
• On Adult (32K instances, classification): 14.9% error in 9 seconds vs. SVMlight's 15.1% in 7 minutes.
• On Forest Cover (522K instances, classification): Binning features achieve 2.2% error in 25 minutes vs. libSVM's 2.2% in 44 hours.

Two additional findings stand out. First, accuracy keeps improving as the training set grows — doubling data reduces error by up to 40%, a luxury that exact kernel machines cannot afford. Second, good performance appears from a remarkably modest number of features: $D = 300$ to $500$ for most datasets.

The complete Random Features workflow is refreshingly simple:

1. Choose your kernel — Gaussian RBF, Laplacian, or any shift-invariant kernel.
2. Compute the spectral distribution $p(\omega)$ — the Fourier transform of the kernel.
3. Sample D random frequencies $\omega_1, \dots, \omega_D$ from $p(\omega)$ and D phases $b_1, \dots, b_D$ from $\text{Uniform}[0, 2\pi]$.
4. Transform every data point: $z(x) = \sqrt{2/D}\,[\cos(\omega_1^\top x + b_1), \dots, \cos(\omega_D^\top x + b_D)]^\top$.
5. Train a linear model on the transformed data (ridge regression, linear SVM, etc.).
6. Predict: $f(x) = w^\top z(x)$ — requires only $O(D + d)$ operations.

The random frequencies are drawn once and reused for all data points. The entire pipeline replaces $O(N^2)$ kernel computation with $O(ND)$ feature computation plus $O(ND^2)$ linear solving.

Look at the Random Fourier Feature map again: $z(x) = \sqrt{2/D}\,\cos(W^\top x + b)$. This is exactly a one-hidden-layer neural network with random, frozen weights $W$ and cosine activation — only the output layer is learned.

This observation is more than a coincidence. The Neural Tangent Kernel (NTK) theory (Jacot et al., 2018) shows that training an infinitely wide neural network with gradient descent is equivalent to kernel regression with a specific kernel determined by the network architecture. Random Features are the finite-width, explicit-map version of this correspondence.

The Random Features paper planted the seed: kernel approximation via random projections is not just a computational trick — it reveals a deep structural connection between kernel methods and neural networks that reshaped how we think about both.