Kernel Methods

A kernel method is a learning algorithm that uses the data only through a kernel $k(\mathbf{x}, \mathbf{x}')$ , a similarity between two inputs that equals an inner product $\phi(\mathbf{x})^\top \phi(\mathbf{x}')$ after some feature map $\phi$ . The algorithm is linear in the feature space, so it keeps the convex training and closed-form solutions of linear models, while its predictions are nonlinear in the input. The feature space can have very many or infinitely many dimensions, because it is never built: the algorithm computes the kernel from the original inputs. Goodfellow, Bengio and Courville (§5.7.2) call the algorithms that work this way “kernel machines, or kernel methods”. The support vector machine is one; its article covers the kernel SVM, Mercer’s condition as Cortes and Vapnik print it, and the table of kernels used in vision.

This article covers the framework those methods share. It explains positive-definite kernels and how to build them, and the kernel trick as a general device. It then treats kernel ridge regression and its relation to linear regression and to Gaussian processes, kernel PCA with centring in feature space, and the $n \times n$ Gram matrix that limits every exact kernel method together with the two standard ways around it: the Nyström approximation and random Fourier features.

In computer vision the word “kernel” has other meanings. A convolution or filter kernel (see convolution), a blur or point-spread kernel, a structuring element in morphology, and the density kernels of kernel density estimation and mean shift are all different objects; Murphy (Chapter 17) lists density, Markov-transition and convolution kernels as other meanings of the word. The “kernel methods” in the background-modelling literature that background subtraction describes are density estimators. GPU and operating-system kernels, and the kernel (null space) of a matrix, are unrelated. This article is about positive-definite kernels only.

Definition

A positive-definite kernel, also called a Mercer kernel, is a symmetric function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ such that, for any finite set of inputs $\mathbf{x}_1, \dots, \mathbf{x}_N$ and any real numbers $c_1, \dots, c_N$ ,

\sum_{i=1}^{N} \sum_{j=1}^{N} k(\mathbf{x}_i, \mathbf{x}_j)\, c_i c_j \ge 0 .

Equivalently, every Gram matrix $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$ is symmetric and has no negative eigenvalues. Terms vary between sources. Murphy (§17.1, eq. 17.1) prints the condition with $\ge 0$ , adds an assumption under which equality needs $c_i = 0$ for all $i$ , and then calls $k$ a Mercer kernel if and only if the Gram matrix is “positive definite” for any set of distinct inputs. The support vector machine article gives Bishop’s finite form, a positive semidefinite Gram matrix. A kernel method is any algorithm that can be written so that the inputs enter only through such kernel values.

Intuition

A linear model can only draw straight boundaries or fit planes. A classic way to get curved ones is to map each input to a longer vector of nonlinear features, such as products of pixel values, and fit a linear model to those. The cost is the length of the vector: all products of two pixels of a $16 \times 16$ patch already make tens of thousands of features.

The kernel trick avoids the long vector. Many linear algorithms, including ridge regression, the SVM and PCA, can be rewritten so that the training inputs appear only in inner products $\mathbf{x}_i^\top \mathbf{x}_j$ . If those inner products are replaced by $k(\mathbf{x}_i, \mathbf{x}_j) = \phi(\mathbf{x}_i)^\top \phi(\mathbf{x}_j)$ , the algorithm runs in the feature space of $\phi$ without ever computing $\phi$ . For the quadratic kernel $(\mathbf{x}^\top \mathbf{x}')^2$ on two-dimensional inputs, Murphy (§17.1.1) writes the map out: $\phi(x_1, x_2) = [x_1^2, \sqrt{2}\, x_1 x_2, x_2^2]$ , three features whose inner product the kernel computes with one dot product and one square.

A second reading treats the kernel as a similarity. Goodfellow and colleagues describe the Gaussian kernel as template matching: each training example is a template, and a test input that lies close to a template in Euclidean distance takes much of that example’s label. A kernel model’s prediction $f(\mathbf{x}) = \sum_i \alpha_i k(\mathbf{x}, \mathbf{x}_i)$ is a weighted vote of the training examples, with weights set by how similar each one is to $\mathbf{x}$ . Murphy opens his kernel chapter in the same way: the model keeps the whole training set and predicts by comparing a new input with each training point.

Formal Definition

Feature spaces and Mercer’s theorem

For a finite set of points the link between a kernel and a feature map follows from linear algebra. Murphy (§17.1.1, eqs. 17.4–17.5) diagonalizes a positive-definite Gram matrix as $K = U^\top \Lambda U$ and sets $\phi(\mathbf{x}_i) = \Lambda^{1/2} U_{:i}$ , so that $K_{ij} = \phi(\mathbf{x}_i)^\top \phi(\mathbf{x}_j)$ . Mercer’s theorem extends this from matrices to kernel functions. For the RBF (Gaussian) kernel

k(\mathbf{x}, \mathbf{x}') = \exp\!\big(-\gamma \lVert \mathbf{x} - \mathbf{x}' \rVert^2\big)

the feature space is infinite-dimensional (Murphy, §17.1.1). Murphy writes the RBF kernel with a length scale, $\exp(-\lVert \mathbf{x} - \mathbf{x}' \rVert^2 / 2\ell^2)$ , so $\gamma = 1/(2\ell^2)$ . scikit-learn and this article use $\gamma$ , and Rahimi and Recht use $\gamma$ for the spherical Gaussian kernel. Goodfellow and colleagues give a simpler infinite example on the non-negative integers: if $\phi(x)$ is $x$ ones followed by zeros, then $k(x, x') = \min(x, x')$ is exactly the infinite-dimensional inner product.

Building kernels

New kernels can be built from valid ones. Murphy (§17.1.2.2–§17.1.2.3) lists the rules. If $k_1$ and $k_2$ are valid kernels, so are $c\,k_1$ for $c > 0$ , $f(\mathbf{x})\, k_1(\mathbf{x}, \mathbf{x}')\, f(\mathbf{x}')$ for any function $f$ , $q(k_1)$ for a polynomial $q$ with non-negative coefficients, $\exp(k_1)$ , $\mathbf{x}^\top A \mathbf{x}'$ for a positive semidefinite matrix $A$ , and the sum $k_1 + k_2$ and product $k_1 k_2$ . Starting from the linear kernel $\mathbf{x}^\top \mathbf{x}'$ , these rules give the polynomial kernel $(\mathbf{x}^\top \mathbf{x}' + c)^M$ , which contains all monomials up to degree $M$ . They also give the Gaussian kernel, written as a product $\exp(-\lVert\mathbf{x}\rVert^2/2\sigma^2)\exp(\mathbf{x}^\top\mathbf{x}'/\sigma^2)\exp(-\lVert\mathbf{x}'\rVert^2/2\sigma^2)$ (eq. 17.24). The kernels used on image features, histogram intersection and χ² among them, and the warning that one common form of the χ² kernel is only conditionally positive definite, are collected in the support vector machine article. Similarity and distance measures compares the matching histogram distances. Murphy (§17.1.2.4) also notes that kernels can be defined on structured inputs, such as strings and graphs, which are hard to turn into fixed-length feature vectors.

The kernel trick

Goodfellow and colleagues (eqs. 5.82–5.83) state the trick for any model whose linear function can be rewritten as $\mathbf{w}^\top \mathbf{x} + b = b + \sum_i \alpha_i \mathbf{x}^\top \mathbf{x}^{(i)}$ . Replacing the inner product by $k$ gives

f(\mathbf{x}) = b + \sum_{i} \alpha_i\, k(\mathbf{x}, \mathbf{x}^{(i)}),

which is nonlinear in $\mathbf{x}$ but linear in $\phi(\mathbf{x})$ and in $\boldsymbol{\alpha}$ . They give two reasons why this is powerful. First, the model is nonlinear in $\mathbf{x}$ , but it is fitted by convex optimization, because $\phi$ is fixed and only $\boldsymbol{\alpha}$ is learned. Second, the kernel is often much cheaper to evaluate than the explicit inner product, and it stays finite when $\phi$ does not. Murphy (§17.3.4) introduces the trick through the SVM dual, which has $N$ unknowns and in general takes $O(N^3)$ time to solve, and notes that the trick also lets classifiers work on strings and graphs.

Kernel ridge regression

Ridge regression shows the step most directly, because its kernel form follows by exact algebra. The linear regression article gives the ridge estimate $\hat{\mathbf{w}} = (X^\top X + \lambda I_D)^{-1} X^\top \mathbf{y}$ . Murphy (§17.3.9) rewrites it with the matrix inversion lemma as

\hat{\mathbf{w}} = X^\top (X X^\top + \lambda I_N)^{-1} \mathbf{y} = X^\top \boldsymbol{\alpha}, \qquad \boldsymbol{\alpha} = (X X^\top + \lambda I_N)^{-1} \mathbf{y} .

The weight vector is a combination of the training inputs, and a prediction needs only inner products, $f(\mathbf{x}) = \sum_n \alpha_n \mathbf{x}_n^\top \mathbf{x}$ . Replacing $XX^\top$ by the Gram matrix $K$ and $\mathbf{x}_n^\top \mathbf{x}$ by $k(\mathbf{x}_n, \mathbf{x})$ gives kernel ridge regression (KRR, Murphy eqs. 17.106–17.108):

f(\mathbf{x}) = \mathbf{k}(\mathbf{x})^\top (K + \lambda I_N)^{-1} \mathbf{y}, \qquad \mathbf{k}(\mathbf{x}) = [k(\mathbf{x}, \mathbf{x}_1), \dots, k(\mathbf{x}, \mathbf{x}_N)]^\top .

The primal solve is $D \times D$ in the number of features. The dual solve is $N \times N$ in the number of examples, whatever the dimension of the feature space. As Murphy notes, $\boldsymbol{\alpha}$ is not sparse, so each prediction costs $O(N)$ kernel evaluations. The SVM’s hinge loss makes its multipliers sparse. Squared loss does not. scikit-learn’s KernelRidge solves the dual system with a Cholesky-based routine (scipy.linalg.solve with assume_a="pos", falling back to least squares with a warning if the matrix is singular) and fits no intercept (kernel_ridge.py, _ridge.py). Its defaults are alpha=1, kernel="linear", gamma=None, degree=3 and coef0=1. KernelRidge also passes gamma to the pairwise kernels, where gamma=None means $1/d$ for the RBF and polynomial kernels (metrics/pairwise.py). That differs from the SVC default 'scale', which the SVM article describes.

The linear regression article notes that scikit-learn’s Ridge already uses this dual route when there are more features than samples. A kernel turns that trick into a nonlinear model.

Gaussian processes

The same prediction appears in Bayesian form. A Gaussian process (GP) puts a prior on functions under which the values $f(\mathbf{x}_1), \dots, f(\mathbf{x}_M)$ at any finite set of inputs are jointly Gaussian, with covariance $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$ for a Mercer kernel $k$ (Murphy, §17.2). With Gaussian observation noise of variance $\sigma_y^2$ and zero prior mean, the posterior mean at a test input is $\mathbf{k}_*^\top (K + \sigma_y^2 I)^{-1} \mathbf{y}$ . Murphy (§17.2.2) notes that this is identical to the KRR prediction; comparing the two formulas, the noise variance takes the place of $\lambda$ . The GP also gives a predictive variance, $k_{**} - \mathbf{k}_*^\top (K + \sigma_y^2 I)^{-1} \mathbf{k}_*$ , and a marginal likelihood for choosing the kernel’s parameters. Murphy (§17.2.4) shows that Bayesian linear regression on features $\phi$ is the GP with kernel $\phi(\mathbf{x})^\top \Sigma_w \phi(\mathbf{x}')$ . He recommends a Cholesky factorization rather than an explicit inverse (§17.2.5), at $O(N^3)$ cost. He also contrasts GP regression with Nadaraya–Watson kernel regression, which uses density kernels, takes $O(N)$ time to weight the training points and does not interpolate the data (§17.2.3). Uncertainty estimation covers other ways to obtain predictive variances.

Kernel PCA

Principal component analysis can be computed from the $N \times N$ Gram matrix $XX^\top$ instead of the $D \times D$ covariance, the route that eigenfaces uses for images with more pixels than examples. Replacing $XX^\top$ by a kernel matrix gives kernel PCA, which Murphy (§20.4.6) attributes to Schölkopf, Smola and Müller (1998). If $U$ and $\Lambda$ hold the top eigenvectors and eigenvalues of the kernel matrix, the projection of a test input onto the components is $\mathbf{k}_*^\top U \Lambda^{-1/2}$ (eq. 20.114).

PCA needs centred data, and the mean in feature space cannot be subtracted directly when $\phi$ is never built. The centred Gram matrix can still be computed: with $C_N = I_N - \tfrac{1}{N}\mathbf{1}\mathbf{1}^\top$ , the Gram matrix of the centred features is $\tilde{K} = C_N K C_N$ (Murphy, double centring), that is, $K$ minus its row means and column means plus its overall mean. A test input is centred with the training statistics. scikit-learn’s KernelCenterer subtracts the training column means of $K$ , subtracts the test row’s own mean and adds the overall training mean (preprocessing/_data.py). With a linear kernel, kernel PCA is ordinary PCA and gives at most $D$ components. A non-degenerate kernel gives up to $N$ (Murphy). Murphy adds, citing later work, that kernel PCA with an RBF kernel expands the feature space rather than reducing it, and he calls its results on his running example “arguably not very useful”. The spectral clustering article describes how Ng, Jordan and Weiss relate their algorithm to kernel PCA.

KernelPCA defaults to kernel="linear", gamma=None (meaning $1/d$ ) and coef0=1. Its eigen_solver="auto" switches to ARPACK only when there are more than 200 samples and fewer than 10 components. It fixes the eigenvector signs with svd_flip and returns the training projections as eigenvectors times the square roots of the eigenvalues (_kernel_pca.py).

Properties

Cost grows with the number of examples, not with the features. The Gram matrix has $N^2$ entries, and an exact solve or eigendecomposition takes $O(N^3)$ time (Murphy, §17.2.9). Goodfellow and colleagues note that building the $m \times m$ matrix alone costs $O(m^2)$ (§5.9) and that evaluating a kernel machine costs time linear in the number of training examples, because each one contributes a term $\alpha_i k(\mathbf{x}, \mathbf{x}^{(i)})$ (§5.7.2). The SVM reduces the evaluation cost by making most $\alpha_i$ zero.
Convex fitting. For a fixed kernel, KRR has a closed-form solution and the SVM a convex objective, Goodfellow and colleagues’ first reason for the trick. Kernel PCA is an eigenproblem whose components, for distinct eigenvalues, are unique up to sign (the example matches signs before comparing). None of them has an initialization to tune.
The kernel is a prior, chosen in advance. The kernel and its parameters ( $\gamma$ , degree, offset, $\lambda$ ) decide what the model can represent and how smoothly it interpolates. Non-probabilistic methods such as SVMs choose them by validation loss (Murphy, §17.2.6.1), for instance by cross-validation, which describes how selecting KRR’s hyperparameters can itself overfit. In a GP they can be chosen by the marginal likelihood (Murphy, §17.2.6).
Local kernels interpolate between nearby examples. Goodfellow and colleagues (§5.11.2) describe kernels such as the Gaussian as local, similarity-based template matching that interpolates between the outputs of nearby training examples. They argue that such methods need on the order of $k$ examples to distinguish $k$ regions of input space, which is part of why they write that kernel machines with generic kernels “struggle to generalize well” (§5.7.2). The curse of dimensionality and k-nearest neighbors articles describe the same limit. The support vector machine article reports Bishop’s warning that a high-dimensional feature space is not a way around it.
Shift-invariant kernels are Fourier transforms of distributions. Bochner’s theorem, as Rahimi and Recht (2007, Theorem 1) state it, says that a continuous shift-invariant kernel $k(\mathbf{x} - \mathbf{y})$ on $\mathbb{R}^d$ is positive definite if and only if it is the Fourier transform of a non-negative measure. This is the basis of random Fourier features below.

Scaling: Nyström and random features

Rahimi and Recht contrast two ways to train a kernel machine with fast linear methods: approximately factor the kernel matrix and use the columns of the factor as features, or factor the kernel function itself.

The Nyström method factors the matrix. Williams and Seeger (2000) choose $m$ of the $n$ rows and columns of $K$ at random without replacement and set

\tilde{K} = K_{n,m} K_{m,m}^{-1} K_{m,n},

a reduced-rank approximation that reduces the cost of a predictor from $O(n^3)$ to $O(m^2 n)$ . On the USPS and abalone datasets they report that $m \ll n$ caused no significant loss of accuracy. With $K_{m,m} = V \operatorname{diag}(s) V^\top$ , the features $\mathbf{z}(\mathbf{x}) = \operatorname{diag}(s)^{-1/2} V^\top \mathbf{k}_m(\mathbf{x})$ have inner products $\tilde{K}$ , so any linear method trained on them approximates the kernel method. scikit-learn’s Nystroem (default kernel RBF, 100 components) draws the $m$ basis points as the first $m$ of a random permutation. It computes $K_{m,m}^{-1/2}$ from a singular value decomposition with singular values floored at $10^{-12}$ (kernel_approximation.py).

Random Fourier features factor the function. By Bochner’s theorem, a suitably scaled shift-invariant kernel is an expectation, $k(\mathbf{x} - \mathbf{y}) = \mathbb{E}_{\boldsymbol{\omega} \sim p}[\cos \boldsymbol{\omega}^\top(\mathbf{x} - \mathbf{y})]$ , where $p$ is the kernel’s Fourier transform. For $\exp(-\gamma\lVert\mathbf{x}-\mathbf{y}\rVert^2)$ , $p$ is a Gaussian with covariance $2\gamma I$ (Rahimi and Recht’s table lists the Gaussian $e^{-\lVert\Delta\rVert^2/2}$ with $p(\boldsymbol{\omega}) = (2\pi)^{-D/2} e^{-\lVert\boldsymbol{\omega}\rVert^2/2}$ , where the table’s $D$ is the input dimension). Drawing $D$ frequencies and averaging gives a Monte Carlo estimate. Algorithm 1 of the paper uses the map

\mathbf{z}(\mathbf{x}) = \sqrt{\tfrac{1}{D}}\,\big[\cos(\boldsymbol{\omega}_1^\top\mathbf{x}), \dots, \cos(\boldsymbol{\omega}_D^\top\mathbf{x}), \sin(\boldsymbol{\omega}_1^\top\mathbf{x}), \dots, \sin(\boldsymbol{\omega}_D^\top\mathbf{x})\big]^\top \in \mathbb{R}^{2D},

for which $\mathbf{z}(\mathbf{x})^\top\mathbf{z}(\mathbf{y}) = \tfrac{1}{D}\sum_j \cos\boldsymbol{\omega}_j^\top(\mathbf{x}-\mathbf{y})$ . (The algorithm box writes the Fourier transform with a factor $1/2\pi$ , the table with $(2\pi)^{-D/2}$ ; only the distribution’s shape matters for sampling.) The paper also names an alternative map, $\sqrt{2}\cos(\boldsymbol{\omega}^\top\mathbf{x} + b)$ with $b$ uniform on $[0, 2\pi]$ , which is the form scikit-learn’s RBFSampler implements as $\sqrt{2/D}\cos(XW + \mathbf{b})$ with $W = \sqrt{2\gamma}\,\mathcal{N}(0, 1)$ . Its docstring calls it a variant of Random Kitchen Sinks, and its gamma defaults to 1.0, not $1/d$ . Because each term $\cos\boldsymbol{\omega}^\top(\mathbf{x}-\mathbf{y})$ lies in $[-1, 1]$ , Hoeffding’s inequality bounds the error for a fixed pair of points:

\Pr\big[\,\lvert \mathbf{z}(\mathbf{x})^\top \mathbf{z}(\mathbf{y}) - k(\mathbf{x}, \mathbf{y}) \rvert \ge \epsilon\,\big] \le 2 \exp(-D\epsilon^2/2) .

Their Claim 1 strengthens this to a uniform bound over all pairs in a compact set, with $D = \Omega\big(\tfrac{d}{\epsilon^2}\log\tfrac{\sigma_p \operatorname{diam}(\mathcal{M})}{\epsilon}\big)$ , where $\sigma_p^2 = 2d\gamma$ for the spherical Gaussian kernel. A linear model on $\mathbf{z}$ is then trained by any fast linear solver and evaluated in $O(D + d)$ operations, compared with $O(Nd)$ for the kernel expansion. The paper’s second construction, random binning features, partitions the input space with randomly shifted grids and suits kernels based on the $L_1$ distance, such as the Laplacian. Murphy (§17.2.9.3) gives the cost of random-feature methods as $O(NM + M^3)$ for $M$ features. He also mentions positive random features, preferable in some applications such as attention, and orthogonal random features, which give a lower-variance estimate. His §17.2.9 adds inducing-point (sparse GP) approximations and Krylov methods based on matrix–vector products, which parallelize well on GPUs.

The two families differ in what they sample. Nyström samples data points, so its basis adapts to where the data lie. Random features sample frequencies from a distribution fixed by the kernel alone, independent of the data. Which is more accurate at a given size depends on the data, and the example below compares them on one dataset.

Examples

The example regresses the orientation of an edge from image patches with known ground truth. Each $16 \times 16$ patch is a Gaussian-blurred step edge written in closed form, $I(x, y) = b + c\,(\Phi(u/\sigma_b) - \tfrac12) + n$ with $u = x\cos\theta + y\sin\theta - o$ and $\Phi$ the standard normal cumulative distribution function. Pixel centres run from $-7.5$ to $7.5$ , with $x$ along columns to the right and $y$ along rows downward. $\theta$ , uniform on $[0, 2\pi)$ , is the direction of the edge normal from the dark to the bright side. The offset $o$ ( $\pm 3$ px), blur $\sigma_b$ (0.7–1.5 px), level $b$ (0.3–0.7) and contrast $c$ (0.3–0.6) are random, and the pixel noise $n$ has standard deviation 0.05. The target is $(\cos 2\theta, \sin 2\theta)$ , the orientation without polarity: an edge and its reversed-polarity twin at $\theta + \pi$ share a target. The predicted orientation is $\tfrac12\operatorname{atan2}$ of the two outputs, and the measure is the mean error in degrees, wrapped into $[0°, 90°]$ , over 1,000 test patches; 2,000 patches are used for training. Inputs are the patches divided by 16. All methods use $\lambda = 0.1$ . The RBF kernel’s $\gamma$ is set by the median heuristic, one over the median squared distance between training patches. The polynomial kernel is $(\mathbf{x}^\top\mathbf{x}' + 1)^2$ .

The construction predicts one result before anything is fitted. $\Phi(u/\sigma_b) - \tfrac12$ is odd in $u$ , and turning $\theta$ by $\pi$ while negating the offset flips the sign of the edge term, which gives an equally likely patch with the same target. A linear function of the pixels without an intercept therefore has no correlation with the target to exploit, and its best population coefficients are zero. This is our derivation. A linear kernel should do no better than chance, a mean error of 45°. A degree-2 kernel contains products of pixels, which are even under the flip, so it can learn the orientation.

The code uses NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 with one BLAS thread. The four blocks form one script, which takes about 20 seconds on a laptop. The first block builds the data, checks it, and gives the orientation from the structure tensor of each test patch (summed products of image gradients) as an estimate that needs no training. This confirms that the orientation can be read from the data.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[v] = "1"                       # one BLAS thread: reproducible sums
import numpy as np
from scipy.linalg import cho_factor, cho_solve, solve, eigh
from scipy.special import ndtr                # standard normal CDF
from sklearn.kernel_ridge import KernelRidge
from sklearn.kernel_approximation import Nystroem, RBFSampler
from sklearn.decomposition import KernelPCA

# ---------- Data: 16x16 patches of a Gaussian-blurred step edge -------------------
P = 16
c = np.arange(P) - (P - 1) / 2                # pixel centres -7.5 .. 7.5
xx, yy = np.meshgrid(c, c)                    # x along columns, y down the rows
rng = np.random.default_rng(0)

def edges(n):
    """Edge normal angle theta in [0, 2 pi) (dark to bright), random offset,
    blur, mean level and contrast, plus pixel noise of s.d. 0.05."""
    th = rng.uniform(0, 2 * np.pi, n)
    off, blur = rng.uniform(-3, 3, n), rng.uniform(0.7, 1.5, n)
    lev, con = rng.uniform(0.3, 0.7, n), rng.uniform(0.3, 0.6, n)
    u = (xx[None] * np.cos(th)[:, None, None] + yy[None] * np.sin(th)[:, None, None]
         - off[:, None, None])
    img = lev[:, None, None] + con[:, None, None] * (ndtr(u / blur[:, None, None]) - 0.5)
    img += rng.normal(0, 0.05, img.shape)
    return img, th

img_tr, th_tr = edges(2000)
img_te, th_te = edges(1000)
X, Xt = img_tr.reshape(2000, -1) / 16, img_te.reshape(1000, -1) / 16   # / sqrt(256)
T = np.column_stack([np.cos(2 * th_tr), np.sin(2 * th_tr)])           # orientation target
lam = 0.1

def angle_error(pred, th):
    """Mean wrapped error (degrees) between 0.5*atan2 of (cos 2t, sin 2t) and theta mod pi."""
    d = np.degrees(0.5 * np.arctan2(pred[:, 1], pred[:, 0]) - th) % 180
    return np.minimum(d, 180 - d).mean()

def sqdist(A, B):
    return np.maximum((A**2).sum(1)[:, None] + (B**2).sum(1)[None] - 2 * A @ B.T, 0)

D2 = sqdist(X, X)
gamma = 1 / np.median(D2[np.triu_indices(2000, 1)])      # median heuristic
rbf = lambda A, B: np.exp(-gamma * sqdist(A, B))
kernels = {"linear": lambda A, B: A @ B.T,
           "poly-2": lambda A, B: (A @ B.T + 1) ** 2,
           "RBF": rbf}

print(f"data: |t| - 1 max {np.abs(np.hypot(*T.T) - 1).max():.1e}; theta "
      f"{np.degrees(th_tr.min()):.1f}..{np.degrees(th_tr.max()):.1f} deg; pixels "
      f"{img_tr.min():.3f}..{img_tr.max():.3f}; RBF gamma = {gamma:.4f}")

# Label check: orientation from the structure tensor, no learning.
gy, gx = np.gradient(img_te, axis=(1, 2))
Jxx, Jyy, Jxy = (gx * gx).sum((1, 2)), (gy * gy).sum((1, 2)), (gx * gy).sum((1, 2))
st = np.column_stack([Jxx - Jyy, 2 * Jxy])                # angle 2*phi of the gradient
print(f"structure tensor (no training): mean test angle error {angle_error(st, th_te):.2f} deg")

data: |t| - 1 max 1.1e-16; theta 0.1..359.8 deg; pixels -0.136..1.161; RBF gamma = 9.0307
structure tensor (no training): mean test angle error 2.35 deg

The targets lie on the unit circle to rounding, the angles cover the full circle and the noise carries the extreme pixel values slightly outside $[0, 1]$ , which the noise-free edge model cannot leave. The data check held. The structure tensor recovers the orientation to 2.35° without any training, so the label check held: the target can be learned from the patches.

The second block shows the kernel trick exactly and then fits KRR with three kernels. For (a) the patches are averaged over $2 \times 2$ blocks to $8 \times 8$ , because the explicit degree-2 map of 256 pixels would have 33,153 features. With 64 pixels it has $1 + 64 + 64 + 2{,}016 = 2{,}145$ : a constant, the pixels times $\sqrt{2}$ , their squares, and every product of two pixels times $\sqrt{2}$ . These are the terms of $(\mathbf{x}^\top\mathbf{x}' + 1)^2$ .

# ---------- (a) The kernel trick: explicit degree-2 map on 8x8 patches -----------
x8 = img_tr.reshape(2000, 8, 2, 8, 2).mean((2, 4)).reshape(2000, -1) / 8
x8t = img_te.reshape(1000, 8, 2, 8, 2).mean((2, 4)).reshape(1000, -1) / 8
iu = np.triu_indices(64, 1)
phi = lambda A: np.hstack([np.ones((len(A), 1)), np.sqrt(2) * A, A**2,
                           np.sqrt(2) * A[:, iu[0]] * A[:, iu[1]]])
F, Ft = phi(x8), phi(x8t)
K8 = (x8 @ x8.T + 1) ** 2
w = solve(F.T @ F + lam * np.eye(F.shape[1]), F.T @ T, assume_a="pos")   # primal
a8 = cho_solve(cho_factor(K8 + lam * np.eye(2000)), T)                   # dual
print(f"\n(a) explicit map: {F.shape[1]} features; max|Phi Phi^T - K|/max|K| = "
      f"{np.abs(F @ F.T - K8).max() / np.abs(K8).max():.1e}; "
      f"max|primal - dual| on test = {np.abs(Ft @ w - ((x8t @ x8.T + 1) ** 2) @ a8).max():.1e}")

# ---------- (b) Kernel ridge regression with three kernels ----------------------
print("\n(b) kernel   mean test angle error   max|own - KernelRidge|")
sk_args = {"linear": dict(kernel="linear"),
           "poly-2": dict(kernel="poly", degree=2, gamma=1.0, coef0=1),
           "RBF": dict(kernel="rbf", gamma=gamma)}
for name, k in kernels.items():
    alpha = cho_solve(cho_factor(k(X, X) + lam * np.eye(2000)), T)
    pred = k(Xt, X) @ alpha
    ref = KernelRidge(alpha=lam, **sk_args[name]).fit(X, T).predict(Xt)
    print(f"    {name:7s} {angle_error(pred, th_te):10.2f} deg {np.abs(pred - ref).max():19.1e}")
    if name == "RBF":
        K, exact_err = k(X, X), angle_error(pred, th_te)

(a) explicit map: 2145 features; max|Phi Phi^T - K|/max|K| = 2.9e-15; max|primal - dual| on test = 3.6e-13

(b) kernel   mean test angle error   max|own - KernelRidge|
    linear       44.58 deg             0.0e+00
    poly-2        1.70 deg             0.0e+00
    RBF           1.16 deg             6.5e-14

In (a) the explicit features reproduce the kernel matrix to $2.9 \times 10^{-15}$ of its largest entry. Primal ridge on the 2,145 features (a $2{,}145 \times 2{,}145$ solve) and dual KRR on the Gram matrix (a $2{,}000 \times 2{,}000$ solve) give the same test predictions to $3.6 \times 10^{-13}$ . Both checks held. The kernel did the work of the feature map with one $64$ -dimensional dot product per pair.

In (b) our Cholesky solve matches KernelRidge to rounding for all three kernels, so the agreement check held for each. The orientation hypothesis held. The linear kernel’s 44.58° is at chance, inside the planned band of 40° to 50°, as the symmetry argument requires. The degree-2 kernel (1.70°) and the RBF kernel (1.16°) are each far below the planned bar of one third of the linear error. Both also beat the structure tensor on the same test patches. Nothing in the inputs changed between the three rows: the kernel alone decides whether the even part of the signal is reachable.

The third block replaces the exact $2{,}000 \times 2{,}000$ RBF Gram matrix by three approximations with $m$ features each, for $m$ from 25 to 1,600 and random seeds 0, 1 and 2. The first is our Nyström method on $m$ random columns, using the same columns as scikit-learn’s Nystroem and checked against its approximate Gram matrix. The second is random Fourier features in the offset form, using RBFSampler’s own weights and offsets and checked against its output. The third is Rahimi and Recht’s Algorithm 1, with $m/2$ frequencies and $m$ features. For each it prints the relative Frobenius error of the approximate Gram matrix, $\lVert\tilde{K} - K\rVert_F/\lVert K\rVert_F$ , and the test angle error of ridge regression on the features, which is KRR with $\tilde{K}$ . It then checks the fixed-pair Hoeffding bound with Algorithm 1 for 50 fixed pairs of training patches, using 500 independent draws of $D$ frequencies for each $D$ .

# ---------- (c) Nystrom and random Fourier features against exact RBF KRR --------
def nystroem(m, seed):
    idx = np.random.RandomState(seed).permutation(2000)[:m]   # Williams & Seeger
    s, V = eigh(K[np.ix_(idx, idx)])
    keep = s > 1e-10 * s.max()
    M = V[:, keep] / np.sqrt(s[keep])                       # K_mm^(-1/2) on its range
    sk = Nystroem(kernel="rbf", gamma=gamma, n_components=m, random_state=seed).fit(X)
    Zs = sk.transform(X)
    Z = K[:, idx] @ M
    chk = np.abs(Z @ Z.T - Zs @ Zs.T).max()
    return Z, rbf(Xt, X[idx]) @ M, chk

def rff_offset(m, seed):
    sk = RBFSampler(gamma=gamma, n_components=m, random_state=seed).fit(X)
    W, b = sk.random_weights_, sk.random_offset_
    z = lambda A: np.sqrt(2 / m) * np.cos(A @ W + b)
    return z(X), z(Xt), np.abs(z(X) - sk.transform(X)).max()

def rff_alg1(m, seed):                                       # Rahimi & Recht, Algorithm 1
    Om = np.random.default_rng(seed).normal(0, np.sqrt(2 * gamma), (m // 2, X.shape[1]))
    z = lambda A: np.sqrt(2 / m) * np.hstack([np.cos(A @ Om.T), np.sin(A @ Om.T)])
    return z(X), z(Xt), 0.0

Knorm, sizes = np.linalg.norm(K), [25, 50, 100, 200, 400, 800, 1600]
print(f"\n(c) exact RBF kernel ridge: mean test angle error {exact_err:.2f} deg")
print("method        m   Gram error mean [min, max]       test angle error mean [min, max]")
worst = {"Nystroem": 0.0, "RBFSampler": 0.0}
slopes = {}
for name, f in [("Nystroem", nystroem), ("RFF offset", rff_offset), ("RFF Alg. 1", rff_alg1)]:
    means = []
    for m in sizes:
        g, e = [], []
        for seed in (0, 1, 2):
            Z, Zt, chk = f(m, seed)
            key = "Nystroem" if name == "Nystroem" else "RBFSampler"
            if name != "RFF Alg. 1":
                worst[key] = max(worst[key], chk)
            g.append(np.linalg.norm(Z @ Z.T - K) / Knorm)
            wz = solve(Z.T @ Z + lam * np.eye(Z.shape[1]), Z.T @ T, assume_a="pos")
            e.append(angle_error(Zt @ wz, th_te))
        means.append(np.mean(g))
        print(f"{name:11s} {m:5d}   {np.mean(g):.4g} [{min(g):.4g}, {max(g):.4g}]"
              f"       {np.mean(e):6.2f} [{min(e):.2f}, {max(e):.2f}]")
    slopes[name] = np.polyfit(np.log(sizes), np.log(means), 1)[0]
print("log-log slope of Gram error against m: "
      + ", ".join(f"{k} {v:.3f}" for k, v in slopes.items()))
print(f"checks against scikit-learn (same draws): max|Gram own - Nystroem| {worst['Nystroem']:.1e}, "
      f"max|features own - RBFSampler| {worst['RBFSampler']:.1e}")

# Hoeffding's bound for a fixed pair (Rahimi & Recht, Section 3), Algorithm 1 map
pairs = np.random.default_rng(5).choice(2000, size=(50, 2), replace=False)
A, B = X[pairs[:, 0]], X[pairs[:, 1]]
kAB = np.exp(-gamma * ((A - B) ** 2).sum(1))
print(f"\nHoeffding, 50 fixed pairs (k from {kAB.min():.3f} to {kAB.max():.3f}), 500 draws each:")
print("   D   eps   largest fraction |z'z - k| >= eps   bound 2 exp(-D eps^2 / 2)")
hr = np.random.default_rng(7)
for D in (25, 100, 400):
    dev = np.empty((500, 50))
    for r in range(500):
        Om = hr.normal(0, np.sqrt(2 * gamma), (D, X.shape[1]))
        za = np.hstack([np.cos(A @ Om.T), np.sin(A @ Om.T)]) / np.sqrt(D)
        zb = np.hstack([np.cos(B @ Om.T), np.sin(B @ Om.T)]) / np.sqrt(D)
        dev[r] = (za * zb).sum(1) - kAB
    for eps in (0.1, 0.2, 0.3):
        print(f"{D:5d}  {eps:.1f}   {(np.abs(dev) >= eps).mean(0).max():25.3f}"
              f"   {min(1.0, 2 * np.exp(-D * eps**2 / 2)):22.3f}")

(c) exact RBF kernel ridge: mean test angle error 1.16 deg
method        m   Gram error mean [min, max]       test angle error mean [min, max]
Nystroem       25   0.0606 [0.05363, 0.06649]         2.90 [2.34, 3.45]
Nystroem       50   0.03214 [0.02906, 0.03403]         1.70 [1.66, 1.74]
Nystroem      100   0.01597 [0.0148, 0.01753]         1.61 [1.55, 1.66]
Nystroem      200   0.007753 [0.007544, 0.007973]         1.57 [1.53, 1.62]
Nystroem      400   0.003934 [0.00385, 0.004015]         1.45 [1.45, 1.46]
Nystroem      800   0.001835 [0.001819, 0.001844]         1.27 [1.26, 1.29]
Nystroem     1600   0.0005435 [0.0005381, 0.0005464]         1.17 [1.17, 1.17]
RFF offset     25   0.3737 [0.3316, 0.442]        15.82 [14.84, 16.92]
RFF offset     50   0.2933 [0.238, 0.375]         7.19 [5.87, 7.98]
RFF offset    100   0.2118 [0.1778, 0.2628]         5.32 [4.73, 5.93]
RFF offset    200   0.1797 [0.1485, 0.23]         3.51 [3.28, 3.75]
RFF offset    400   0.1122 [0.1002, 0.1203]         2.85 [2.82, 2.89]
RFF offset    800   0.06745 [0.06456, 0.06918]         2.37 [2.34, 2.43]
RFF offset   1600   0.04925 [0.04411, 0.05503]         2.05 [1.98, 2.14]
RFF Alg. 1     25   0.3724 [0.2901, 0.4343]        13.16 [8.96, 16.83]
RFF Alg. 1     50   0.2649 [0.224, 0.3467]         7.28 [6.75, 8.18]
RFF Alg. 1    100   0.1758 [0.16, 0.1949]         4.44 [3.83, 5.02]
RFF Alg. 1    200   0.1234 [0.1129, 0.1391]         3.45 [3.30, 3.67]
RFF Alg. 1    400   0.08997 [0.07221, 0.1158]         2.76 [2.69, 2.84]
RFF Alg. 1    800   0.06352 [0.0549, 0.07947]         2.37 [2.34, 2.40]
RFF Alg. 1   1600   0.04531 [0.03682, 0.05824]         1.98 [1.94, 2.06]
log-log slope of Gram error against m: Nystroem -1.096, RFF offset -0.497, RFF Alg. 1 -0.507
checks against scikit-learn (same draws): max|Gram own - Nystroem| 5.1e-13, max|features own - RBFSampler| 0.0e+00

Hoeffding, 50 fixed pairs (k from 0.074 to 0.794), 500 draws each:
   D   eps   largest fraction |z'z - k| >= eps   bound 2 exp(-D eps^2 / 2)
   25  0.1                       0.508                    1.000
   25  0.2                       0.180                    1.000
   25  0.3                       0.044                    0.649
  100  0.1                       0.176                    1.000
  100  0.2                       0.010                    0.271
  100  0.3                       0.000                    0.022
  400  0.1                       0.012                    0.271
  400  0.2                       0.000                    0.001
  400  0.3                       0.000                    0.000

Both scikit-learn checks held: our Nyström matrices agree with Nystroem’s to $5.1 \times 10^{-13}$ and our offset features equal RBFSampler’s exactly.

For random features, each entry of $\tilde{K}$ is an unbiased average of independent bounded terms, so the expected squared Frobenius error is proportional to $1/m$ and the error itself should fall as $m^{-1/2}$ . The planned bar was a log–log slope between −0.55 and −0.45. The hypothesis held for both forms, with slopes −0.497 (offset) and −0.507 (Algorithm 1). We also derived the per-entry variances before the run. For a pair with kernel value $k$ it is $(1 - k^2)^2/m$ for Algorithm 1 and $(1 - k^2 + k^4/2)/m$ for the offset form, because the offset form carries an extra term $\cos(\boldsymbol{\omega}^\top(\mathbf{x}+\mathbf{y}) + 2b)$ of variance ½. The diagonal is exact for Algorithm 1. The prediction that Algorithm 1 has the lower seed-mean Gram error at every size held at all seven sizes. At $m = 25$ the margin is small (0.3724 against 0.3737) and the seed ranges overlap, so that row alone would not separate the two forms. The test error follows the Gram error only loosely: Algorithm 1 has the lower Gram error at 50 features but the higher test error (7.28° against 7.19°).

The comparison of Nyström with random features was planned without a predicted direction and is reported as printed. At every size Nyström has both the lower Gram error and the lower test error, and its Gram error falls with a slope of −1.096 instead of −0.5. In our reading, two features of the construction point the same way, and the printed output does not separate them. Nyström’s basis consists of the kernel columns of actual training patches, while random features draw their frequencies without looking at the data and estimate every entry, the diagonal included for the offset form, with a variance proportional to $1/m$ . And $\tilde{K}$ reproduces the sampled columns $K_{n,m}$ exactly ( $K_{n,m}K_{m,m}^{-1}K_{m,m} = K_{n,m}$ ), so the part of the matrix that is approximated shrinks as $m$ approaches $n = 2{,}000$ . That is consistent with the steeper slope. Because Nyström uses the data’s own columns, this comparison says nothing about larger datasets, where $m$ is a small fraction of $n$ .

Both planned test-error hypotheses held for all three methods. The error at 1,600 features is below that at 25: 1.17° against 2.90° (Nyström), 2.05° against 15.82° (offset) and 1.98° against 13.16° (Algorithm 1). At 1,600 features each method is within 2° of exact KRR’s 1.16°. Nyström matches it to 0.01°. The random-feature models remain 0.89° (offset form, 2.05°) and 0.82° (Algorithm 1, 1.98°) worse, although their relative Gram errors are only 0.04925 and 0.04531.

The Hoeffding check held in all nine cells. In three cells the bound is 1 and says nothing. In the others the observed fractions are far below it: for $D = 100$ and $\epsilon = 0.2$ , at most 1.0% of draws exceed the tolerance for any pair, against a bound of 27.1%. For $D = 400$ with $\epsilon = 0.2$ and $\epsilon = 0.3$ , the largest fraction prints as 0.000. With 500 draws per pair a single exceedance would print as 0.002, so no draw exceeded $\epsilon$ for any pair. The bound $2e^{-D\epsilon^2/2}$ is positive, so the check holds even where the bound rounds to 0.000. The gap is expected: Hoeffding’s inequality uses only the range $[-1, 1]$ of each term, while the actual variance of $\cos\boldsymbol{\omega}^\top\Delta$ is $(1 - k^2)^2/2$ , at most ½.

The last block runs kernel PCA with the RBF kernel and checks it against KernelPCA. Our version centres the Gram matrix in feature space, takes the top 10 eigenpairs and centres the test kernel rows with the training means. It then prints how strongly the first four components follow the edge’s direction and orientation.

# ---------- (d) Kernel PCA with centring in feature space -------------------------
n, q = 2000, 10
Kc = K - K.mean(0) - K.mean(1)[:, None] + K.mean()            # C K C
lam_k, U = eigh(Kc, subset_by_index=(n - q, n - 1))
lam_k, U = lam_k[::-1], U[:, ::-1]
Kt = rbf(Xt, X)
Ktc = Kt - K.mean(0) - Kt.mean(1)[:, None] + K.mean()          # test rows, training means
Y, Yt = U * np.sqrt(lam_k), Ktc @ U / np.sqrt(lam_k)
kp = KernelPCA(n_components=q, kernel="rbf", gamma=gamma, eigen_solver="dense")
Ys = kp.fit_transform(X)
sgn = np.sign((Y * Ys).sum(0))
Yts = kp.transform(Xt)
print(f"\n(d) kernel PCA: max eigenvalue rel. diff {np.abs(lam_k / kp.eigenvalues_ - 1).max():.1e}; "
      f"projections max diff / scale: train {np.abs(Y * sgn - Ys).max() / np.abs(Ys).max():.1e}, "
      f"test {np.abs(Yt * sgn - Yts).max() / np.abs(Yts).max():.1e}")
print(f"first {q} eigenvalues hold {lam_k.sum() / np.trace(Kc):.3f} of the centred Gram's trace")
print("component  |corr| with cos t, sin t, cos 2t, sin 2t (test patches)")
for j in range(4):
    r = [abs(np.corrcoef(Yt[:, j], f(g * th_te))[0, 1]) for g in (1, 2) for f in (np.cos, np.sin)]
    print(f"    {j + 1}       " + "  ".join(f"{v:.2f}" for v in r))

(d) kernel PCA: max eigenvalue rel. diff 3.6e-15; projections max diff / scale: train 3.6e-14, test 3.7e-14
first 10 eigenvalues hold 0.735 of the centred Gram's trace
component  |corr| with cos t, sin t, cos 2t, sin 2t (test patches)
    1       0.97  0.19  0.02  0.02
    2       0.14  0.97  0.03  0.02
    3       0.10  0.10  0.02  0.03
    4       0.02  0.06  0.25  0.05

Eigenvalues and projections agree with scikit-learn to about $10^{-14}$ after each component’s sign is matched, for training and for test patches, so the check held. Note that the test patches must be centred with the training statistics. The first two kernel principal components track $\cos\theta$ and $\sin\theta$ (|r| = 0.97), the direction of the edge normal including its polarity, which the regression target deliberately discards. Only the fourth component shows a weak relation to $\cos 2\theta$ (0.25). Kernel PCA is unsupervised: it finds the directions of largest variance in feature space, and here, in our reading, polarity is the largest source of variation between patches. Ten components hold 73.5% of the centred Gram matrix’s trace.

The example was planned before it was run. The hypotheses and bars above were written first. A planning script then ran on the same data and printed only the floors of the agreement checks, with the prediction sizes and eigenvalue ratios they depend on, and the data and label checks, from which every agreement tolerance was set to $10^{-10}$ . It printed no hypothesis value before the bars were fixed. Every planned check and hypothesis held, and the script was not changed after its first run. Compared with the original plan for this article, the example (1) runs the explicit-map check on $8 \times 8$ block averages, for the memory reason above; (2) adds Rahimi and Recht’s Algorithm 1 next to the RBFSampler form, because the paper states its Hoeffding bound for that map, and compares methods at equal numbers of features; (3) checks our Nyström method through the approximate Gram matrix, because feature bases can differ by a rotation, and checks our offset features with the sampler’s own draws; (4) adds the structure-tensor label check; (5) adds the orientation hypothesis for the three kernels, the derived comparison of the two random-feature forms and the two test-error hypotheses, and states the $1/\sqrt{m}$ expectation as a slope bar; (6) adds the informational correlation print to kernel PCA; (7) fixes $\lambda$ , $\gamma$ , the polynomial kernel and the seeds in advance; (8) defines the held-out measure as the wrapped angle error.

Common Misconceptions

“Kernel” means one thing. Convolution filters, blur kernels and the density kernels of mean shift and kernel density estimation do not need to be positive definite, and they are not inner products in a feature space. Only Mercer kernels support the kernel trick.
Any similarity function is a kernel. A function that is not positive definite can give a Gram matrix with negative eigenvalues. The SVM’s dual then loses its guarantees, and $K + \lambda I$ may not be positive definite. Use the construction rules above, or check the Gram matrix’s eigenvalues, as the support vector machine article does for its histogram kernels.
The infinite feature space removes the curse of dimensionality. A Gaussian kernel interpolates between nearby training examples, and Goodfellow and colleagues tie the need for many examples to that locality. The mapped features of a quadratic kernel lie on a low-dimensional surface, as the SVM article reports from Bishop.
Kernel methods are only SVMs. Ridge regression and PCA have kernel forms, Murphy (§17.1) mentions linear discriminant analysis in an implicit feature space, Gaussian processes use kernels as covariances, and the maximum mean discrepancy that domain adaptation uses is defined through a kernel. Of the methods here, the SVM is the one whose solution is sparse.
Random features approximate the kernel method well once the Gram matrix looks close. In the example, Algorithm 1 with 1,600 features had a relative Gram error of 0.04531 but a test error of 1.98° against exact KRR’s 1.16°, a much larger relative gap. Check the approximation on the downstream measure, not only on the kernel matrix.
Kernel PCA components are the features you want. They are the directions of largest variance in feature space. In the example, the leading components encode the edge’s polarity, which the task ignores.
Library defaults agree. KernelRidge and KernelPCA default to the linear kernel; rbf_kernel and polynomial_kernel take $\gamma = 1/d$ when it is not given; RBFSampler takes $\gamma = 1$ ; and SVC uses 'scale' and coef0=0, while the other estimators use coef0=1. Set the kernel parameters explicitly.

Where It Is Used

Classification with SVMs. From the mid-2000s to the mid-2010s the SVM was the usual classifier in vision pipelines, linear on HOG features and with intersection or χ² kernels on histograms of visual words and spatial pyramids. The support vector machine, bag of visual words and Caltech-101 articles describe them. The SVM article also describes the explicit feature maps of Vedaldi and Zisserman for additive kernels, which in our reading play for histogram kernels the role that random Fourier features play here.
Tracking. Kernelized correlation filters train kernel ridge regression on all cyclic shifts of a patch and solve it in the Fourier domain. Correlation filter tracking derives the closed form and checks it against an explicit kernel matrix.
Comparing distributions. The maximum mean discrepancy, a distance between the mean embeddings of two samples in a reproducing-kernel Hilbert space, is used in domain adaptation.
Dimensionality reduction and clustering. Kernel PCA is listed among the variants of principal component analysis, and spectral clustering relates the Ng–Jordan–Weiss algorithm to it.
Coordinate networks. Random Fourier features of pixel coordinates are used as a positional encoding, for instance in SAM’s prompt encoder and in neural implicit representations. There the features feed a network rather than approximating a fixed kernel for a linear model.
Regression with uncertainty. Gaussian processes give KRR’s prediction together with a variance, and Murphy (§17.2.8) notes that a network with one infinitely wide layer of RBF units is equivalent to a GP with an RBF kernel, and that many kinds of deep networks in the infinite limit correspond to a GP with the neural tangent kernel.

Goodfellow and colleagues (§5.7.2) write that the modern incarnation of deep learning was designed to overcome the limitations of kernel machines listed above: evaluation and training costs that grow with the training set, and the poor generalization of generic kernels. They date the current renaissance to Hinton and colleagues’ 2006 network, which outperformed an RBF-kernel SVM on MNIST. The deep learning article summarizes their account of the years from the mid-1990s to 2007, when kernel machines and graphical models achieved good results on many tasks and neural networks fell out of favour.