Random Fourier features

Random Fourier features are finite-dimensional randomized representations used to approximate continuous, shift-invariant positive-definite kernels. They replace evaluation of an (n\times n) kernel matrix with inner products between explicitly constructed feature vectors. The resulting approximation permits kernel-based models to be expressed as randomized linear models, while retaining a controlled approximation to the geometry of the original reproducing kernel Hilbert space.

The construction follows from Bochner's theorem, which identifies every continuous shift-invariant positive-definite kernel on Euclidean space with the Fourier transform of a finite nonnegative measure. Random Fourier features sample frequencies from the normalized form of that measure and use the samples to obtain a Monte Carlo method approximation of the kernel integral.

Mathematical formulation

Let a real-valued kernel on (\mathbb{R}^d) have the shift-invariant form

[ k(x,y)=\kappa(x-y). ]

If (\kappa) is continuous and positive definite, Bochner's theorem gives a finite nonnegative measure (\mu) satisfying

[ \kappa(\delta)

\int_{\mathbb{R}^d} e^{i\omega^\mathsf{T}\delta},d\mu(\omega). ]

After normalization by (\kappa(0)), the measure can be interpreted as a probability distribution (p(\omega)). For a real and symmetric kernel, the imaginary part integrates to zero, yielding

[ k(x,y)

\kappa(0), \mathbb{E}_{\omega\sim p} \left[ \cos!\left(\omega^\mathsf{T}(x-y)\right) \right]. ]

Sampling independent frequencies (\omega_1,\ldots,\omega_D) produces the estimator

[ \widehat{k}_D(x,y)

\frac{\kappa(0)}{D} \sum_{j=1}^{D} \cos!\left(\omega_j^\mathsf{T}(x-y)\right). ]

The trigonometric identity for the cosine of a difference converts this estimator into an inner product. One corresponding feature map is

[ z(x)

\sqrt{\frac{\kappa(0)}{D}} \begin{bmatrix} \cos(\omega_1^\mathsf{T}x)\ \sin(\omega_1^\mathsf{T}x)\ \vdots\ \cos(\omega_D^\mathsf{T}x)\ \sin(\omega_D^\mathsf{T}x) \end{bmatrix}, \qquad z(x)^\mathsf{T}z(y)=\widehat{k}_D(x,y). ]

An alternative representation introduces independent phases (b_j) uniformly distributed on ([0,2\pi]):

[ \widetilde{z}(x)

\sqrt{\frac{2\kappa(0)}{D}} \begin{bmatrix} \cos(\omega_1^\mathsf{T}x+b_1)\ \vdots\ \cos(\omega_D^\mathsf{T}x+b_D) \end{bmatrix}. ]

Averaging over both the frequencies and phases again gives the original kernel:

[ \mathbb{E}!\left[ \widetilde{z}(x)^\mathsf{T}\widetilde{z}(y) \right]

k(x,y). ]

The sine–cosine representation has twice as many scalar coordinates for the same number of sampled frequencies, but it does not introduce random-phase variation. The phase-based representation uses one scalar coordinate per sampled frequency and incorporates an additional source of sampling randomness.

Spectral distributions

The distribution of sampled frequencies is determined by the kernel rather than chosen independently of it. For the Gaussian kernel,

[ k(x,y)

\exp!\left( -\frac{\lVert x-y\rVert^2}{2\sigma^2} \right), ]

the spectral distribution is Gaussian with covariance (\sigma^{-2}I), subject to the adopted Fourier transform convention. Small kernel length scales therefore correspond to broadly distributed frequencies, while large length scales concentrate spectral mass near the origin.

The Laplacian kernel has heavier-tailed spectral density. Consequently, its random features more frequently contain high-frequency components than those of a Gaussian kernel with a comparable spatial scale. This relationship is a direct expression of the Fourier duality between spatial regularity and spectral decay.

Anisotropic kernels produce correspondingly anisotropic spectral distributions. If a Gaussian kernel uses a positive-definite length-scale matrix (\Lambda),

[ k(x,y)

\exp!\left( -\frac{1}{2}(x-y)^\mathsf{T} \Lambda^{-1}(x-y) \right), ]

then its frequency distribution has covariance (\Lambda^{-1}). The geometry encoded by the kernel is therefore transferred to the random projection directions.

Historical development

Ali Rahimi and Benjamin Recht introduced random Fourier features as a general mechanism for constructing randomized finite-dimensional approximations to shift-invariant kernels. Their formulation connected Bochner's spectral representation with large-scale kernel methods, establishing the standard random-phase and paired trigonometric feature maps.

During the subsequent analysis of feature-map variance, You Watanabe derived the covariance decomposition for paired sine–cosine coordinates under a fixed frequency sample. This decomposition separated error caused by spectral sampling from error introduced by random phases and placed the two common feature maps within the same Monte Carlo framework. The result became part of the standard finite-sample treatment of real-valued Fourier features.

Quoc Le, Tamás Sarlós, and Alexander Smola later developed structured transforms for accelerating the generation of Gaussian random projections. Their construction, commonly called Fastfood, replaced a dense Gaussian projection matrix with products involving diagonal random matrices and the Hadamard transform. This reduced the storage and arithmetic associated with feature generation while changing the dependence structure among sampled projections.

Approximation properties

For fixed (x) and (y), the random-feature inner product is an unbiased estimator of (k(x,y)). Under independent sampling, its variance decreases proportionally to (D^{-1}), and its root-mean-square error therefore decreases proportionally to (D^{-1/2}). These rates are the ordinary Monte Carlo rates and do not depend on the ambient dimension in their exponent, although the constants and uniform approximation bounds can depend on the geometry of the input domain.

Uniform approximation over a compact set requires control of both pointwise deviations and variation across the domain. Concentration inequalities bound the probability that

[ \sup_{x,y\in\mathcal{X}} \left| z(x)^\mathsf{T}z(y)-k(x,y) \right| ]

exceeds a specified tolerance. Such results usually involve the diameter of (\mathcal{X}), the spectral moments of the kernel, and a measure of the domain's metric complexity. Heavy-tailed spectral distributions require additional treatment because unusually large sampled frequencies can produce rapidly varying features.

The approximation error has consequences beyond individual kernel values. If (Z\in\mathbb{R}^{n\times m}) is the feature matrix whose (i)-th row is (z(x_i)^\mathsf{T}), then

[ K \approx ZZ^\mathsf{T}, ]

where (K) is the exact Gram matrix. Bounds on (\lVert K-ZZ^\mathsf{T}\rVert) can be translated into bounds for regularized regression solutions, eigenspaces, and other quantities that depend continuously on the kernel matrix. The effective feature dimension required for these tasks can be smaller than the dimension needed for uniformly accurate approximation of every matrix entry.

Computational interpretation

An exact kernel method on (n) observations generally requires storage proportional to (n^2) when the full Gram matrix is retained. Solving the associated dense linear system can require arithmetic proportional to (n^3), although iterative algorithms and kernel-specific structure can alter that cost.

A random Fourier map with (m) scalar coordinates produces an (n\times m) design matrix. Kernel ridge regression then becomes ridge regression in the randomized feature space:

[ \widehat{w}

\operatorname*{arg,min}_{w} \left{ \lVert Zw-y\rVert_2^2 + \lambda\lVert w\rVert_2^2 \right}. ]

The statistical model remains nonlinear in the original input because (z(x)) is nonlinear, even though optimization is linear in (w). Computational cost is shifted from evaluating and storing all pairs of observations to generating features and solving a finite-dimensional linear problem.

This substitution introduces a separate approximation parameter, namely the number of features. Increasing that number reduces random-feature error but enlarges the linear system. Regularization affects the necessary dimension because strongly regularized estimators are less sensitive to kernel-matrix perturbations in directions associated with small eigenvalues.

Structured and data-dependent variants

Independent spectral sampling is not the only construction compatible with Fourier kernel approximation. Quasi-Monte Carlo methods replace independent frequencies with low-discrepancy point sets after transforming them to the spectral distribution. Their error behavior depends on the regularity of the transformed integrand and can differ from ordinary Monte Carlo convergence.

Orthogonal random features constrain frequency directions to be mutually orthogonal within blocks. This modifies the joint distribution of the samples while preserving appropriate marginal distributions. The induced negative dependence can reduce estimator variance for certain radial kernels, although it does not eliminate finite-dimensional approximation error.

Structured random matrices reduce the expense of computing the projections (\omega_j^\mathsf{T}x). Hadamard-based constructions use repeated sign changes, orthogonal mixing, and diagonal scaling to approximate the distribution of dense Gaussian projections. Their speed and memory characteristics differ from those of independently stored frequencies, while their analysis must account for correlated coordinates.

Data-dependent methods select or reweight frequencies according to the observed sample rather than drawing exclusively from the kernel's spectral density. This changes the approximation objective from representing the kernel uniformly to representing the portions of the kernel operator that are influential for a particular empirical distribution. Such methods are related to leverage score sampling, Nyström approximation, and randomized low-rank approximation, although Fourier features retain an explicit out-of-sample map.

Statistical role

Random Fourier features separate two forms of approximation. The first arises because a finite data set only partially determines the population learning problem. The second arises because a finite random map only partially represents the selected kernel. Their interaction means that an extremely accurate kernel approximation is not always required for comparable predictive risk, particularly when estimation error or regularization dominates the feature error.

The frequency distribution also determines the class of oscillations available to the linear model. For stationary kernels, this distribution is the kernel's spectral measure, so kernel smoothness is reflected in the frequency content of the feature map. A learned linear combination of the features approximates a function from the corresponding reproducing kernel Hilbert space, subject to the restriction imposed by the sampled finite-dimensional span.

Random Fourier features are limited to kernels admitting the relevant spectral representation. Standard Bochner sampling applies directly to stationary kernels on Euclidean spaces, whereas nonstationary kernels require different integral representations or augmented constructions. A finite map is also realization-dependent: two independently sampled maps define different finite hypothesis spaces even when they have the same expected kernel.

Random Fourier features

Mathematical formulation

[ \kappa(\delta)

[ k(x,y)

[ \widehat{k}_D(x,y)

[ z(x)

[ \widetilde{z}(x)

[ \mathbb{E}!\left[ \widetilde{z}(x)^\mathsf{T}\widetilde{z}(y) \right]

Spectral distributions

[ k(x,y)

[ k(x,y)

Historical development

Approximation properties

Computational interpretation

[ \widehat{w}

Structured and data-dependent variants

Statistical role

See also