Kernels

The support vector machines note ended with two observations. The dual problem and the resulting classifier use the training points only through inner products $\langle x_i, x_j\rangle$, and the margin bound does not depend on the dimension. A linear classifier on the raw features still draws a hyperplane, and two concentric rings defeat it. This note replaces the raw features by a richer representation $\phi(x)$ and the inner product by a kernel $k(x, z) = \langle \phi(x), \phi(z)\rangle$ that can be computed without ever forming $\phi(x)$. The same quadratic program then draws curved boundaries, and the same idea turns ridge regression into a flexible nonparametric regressor. The price is paid in the number of training points instead of the number of features, and the end of the note shows how to pay less.

Notation. Training points $x_1, \dots, x_n \in \mathbb R^d$ with targets $y_i$. A feature map $\phi: \mathbb R^d \to \mathcal H$ sends a point to a vector in a space $\mathcal H$ with an inner product, possibly of infinite dimension. A kernel is a function $k(x, z)$; the Gram matrix $K$ is the $n \times n$ matrix with entries $K_{ij} = k(x_i, x_j)$, and $k(x) = (k(x_1, x), \dots, k(x_n, x))^\top$ collects the kernel values of a new point against the training points. $\sigma$ is the RBF bandwidth and $\lambda$ the ridge penalty.

Feature maps#

A linear model becomes nonlinear in $x$ when it is linear in nonlinear features of $x$. Points on a line with one class in the middle cannot be split by one threshold, but after $\phi(x) = (x, x^2)$ they lie on a parabola, where a line splits them:

Figure

A feature map makes the classes linearly separable

The line $w_1 x + w_2 x^2 + b = 0$ in the lifted plane is a quadratic inequality in $x$, which holds between two thresholds. In general one would add all products of features up to some degree $p$. The number of monomials of degree at most $p$ in $d$ variables is $\binom{d + p}{p}$: for $d = 100$ and $p = 3$ it is $176{,}851$. Writing these features down for every training point, and fitting a weight to each, quickly becomes the bottleneck.

The kernel trick in the dual#

The SVM never needs the features themselves. Its dual, with $\phi(x_i)$ in place of $x_i$, reads

$$\begin{aligned} &\max_{\alpha}\ \sum_i \alpha_i - \tfrac12 \sum_{i,j} \alpha_i \alpha_j y_i y_j\, k(x_i, x_j)\\ &\text{subject to}\ \ 0 \le \alpha_i \le C,\quad \sum_i \alpha_i y_i = 0, \end{aligned}$$

and the classifier is $f(x) = \sum_i \alpha_i y_i\, k(x_i, x) + b$, with $b = y_i - \sum_j \alpha_j y_j k(x_j, x_i)$ for any support vector on the margin. Both need only $k(x, z) = \langle \phi(x), \phi(z)\rangle$. Replacing every inner product of a linear method by a kernel is the kernel trick (Boser, Guyon and Vapnik, 1992); it works for any method that can be written in terms of inner products, from the SVM and ridge regression to PCA and $k$-means.

The saving comes from kernels that are cheap to evaluate although $\phi$ is large. For $x \in \mathbb R^2$ and $\phi(x) = (x_1^2,\ \sqrt2\, x_1 x_2,\ x_2^2)$,

$$\langle \phi(x), \phi(z)\rangle = x_1^2 z_1^2 + 2 x_1 x_2 z_1 z_2 + x_2^2 z_2^2 = (x_1 z_1 + x_2 z_2)^2 = \langle x, z\rangle^2 .$$

Worked example. For $x = (1, 2)$ and $z = (3, 1)$: $\phi(x) = (1, 2\sqrt2, 4)$, $\phi(z) = (9, 3\sqrt2, 1)$, and $\langle\phi(x), \phi(z)\rangle = 9 + 12 + 4 = 25 = (1 \cdot 3 + 2 \cdot 1)^2$. The kernel costs one inner product in $\mathbb R^2$ and a square; in $d$ dimensions $\langle x, z\rangle^2$ costs $d$ multiplications against $d^2$ for the explicit features. The feature map of the figure, $\phi(x) = (x, x^2)$, has the kernel $k(x, z) = xz + x^2 z^2$.

Which functions are kernels#

Not every similarity function is an inner product in some space. If $k(x, z) = \langle\phi(x), \phi(z)\rangle$, then for any points $x_1, \dots, x_m$ and any real numbers $c_1, \dots, c_m$,

$$\sum_{i,j} c_i c_j\, k(x_i, x_j) = \Big\lVert \sum_i c_i \phi(x_i)\Big\rVert^2 \ge 0 .$$

So the Gram matrix of a kernel on any finite set of points is symmetric and positive semidefinite. The converse also holds: a symmetric function whose Gram matrices are all positive semidefinite is the inner product of some feature map into a Hilbert space, the reproducing kernel Hilbert space (RKHS) of $k$ (Aronszajn, 1950; Shalev-Shwartz and Ben-David, 2014, Lemma 16.2). Such functions are called positive semidefinite kernels, and they are exactly the functions the kernel trick may use. Mercer’s theorem gives the feature map explicitly for a continuous kernel on a compact domain: $k(x, z) = \sum_j \lambda_j e_j(x) e_j(z)$ with eigenvalues $\lambda_j \ge 0$ and orthonormal eigenfunctions $e_j$, so $\phi(x) = (\sqrt{\lambda_j}\, e_j(x))_j$ (Rasmussen and Williams, 2006, §4.3).

The condition is familiar from probability. A matrix is positive semidefinite exactly when it is the covariance matrix of some random vector, because $\sum_{i,j} c_i c_j \operatorname{Cov}(Z_i, Z_j) = \operatorname{Var}(\sum_i c_i Z_i) \ge 0$ (see covariance). In the same way, $k$ is a positive semidefinite kernel exactly when it is the covariance function of a random function: $k(x, z) = \operatorname{Cov}(F(x), F(z))$ for a Gaussian process $F$. A kernel says how strongly the values of an unknown function at two inputs are correlated; this reading returns in kernel ridge regression.

Checking the definition directly is often hard, so kernels are usually built from simpler ones by rules that preserve positive semidefiniteness. If $k_1$ and $k_2$ are kernels, so are:

• $k_1 + k_2$ and $c\, k_1$ for $c \gt 0$: sums and positive multiples of positive semidefinite matrices are positive semidefinite;
• $k_1 k_2$: the elementwise product of positive semidefinite matrices is positive semidefinite (the Schur product theorem), or, with features, $\phi_1 \otimes \phi_2$;
• $g(x) g(z)$ for any real function $g$: a rank-one Gram matrix $gg^\top$;
• $k_1(\psi(x), \psi(z))$ for any map $\psi$;
• the pointwise limit of kernels, because each inequality $\sum c_i c_j k(x_i, x_j) \ge 0$ survives the limit.

Polynomials with non-negative coefficients of a kernel, and its exponential (a limit of such polynomials), are therefore kernels too (Bishop, 2006, §6.2).

Some natural similarities fail. The distance $k(x, z) = \lVert x - z\rVert$ has zeros on the diagonal of its Gram matrix, so the eigenvalues sum to the trace, 0, and unless all points coincide one of them is negative. For two points at 0 and 1 the Gram matrix is $\begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix}$ with eigenvalues $\pm 1$. On 50 random points in the plane its smallest eigenvalue is $-27.5$; the “sigmoid kernel” $\tanh(\langle x, z\rangle + 1)$ gives $-6.26$, while the RBF kernel stays at $0$ up to rounding. With such a function the SVM dual is no longer concave and its solution loses its meaning.

Common kernels: polynomial and RBF#

The polynomial kernel $k(x, z) = (\langle x, z\rangle + c)^p$ with $c \ge 0$ expands by the binomial formula into $\sum_{j=0}^p \binom{p}{j} c^{p-j} \langle x, z\rangle^j$, a non-negative combination of the kernels $\langle x, z\rangle^j$. Its features are all monomials of degree at most $p$, the monomial $x_1^{a_1} \cdots x_d^{a_d}$ of degree $j = a_1 + \dots + a_d$ weighted by $\sqrt{\binom{p}{j} \binom{j}{a_1, \dots, a_d} c^{p - j}}$, where the multinomial coefficient comes from expanding $\langle x, z\rangle^j$; $c$ controls how much the low degrees count relative to the high ones.

The RBF (Gaussian) kernel is

$$k(x, z) = \exp\Big(-\frac{\lVert x - z\rVert^2}{2\sigma^2}\Big) = e^{-\lVert x\rVert^2/(2\sigma^2)}\, e^{-\lVert z\rVert^2/(2\sigma^2)}\, e^{\langle x, z\rangle/\sigma^2} .$$

The first two factors are $g(x) g(z)$ and the last is the exponential of a kernel, so by the rules above it is a kernel. Its Taylor series contains $\langle x, z\rangle^j$ for every $j$, so its feature space has infinite dimension: no explicit $\phi$ could be stored, yet $k$ costs one distance. Every point has $\lVert\phi(x)\rVert^2 = k(x, x) = 1$, and the squared distance between two points in feature space is $k(x, x) + k(z, z) - 2k(x, z) = 2 - 2k(x, z)$: the RBF kernel is a similarity that decays with distance on the scale $\sigma$. A classifier $f(x) = \sum_i \alpha_i y_i\, k(x_i, x) + b$ with this kernel is a weighted vote of the support vectors in which nearby ones count most, a smooth relative of nearest neighbours with weights fitted by margin maximization.

The bandwidth as bias–variance. A large $\sigma$ makes every training point influence the whole space, and the decision function is nearly linear; a small $\sigma$ makes each point influence only its immediate neighbourhood, and the function can wrap around every training point. On two noisy moons (600 points, noise 0.3, half for testing, $C = 1$; scikit-learn writes the kernel as $e^{-\gamma\lVert x - z\rVert^2}$, so $\gamma = 1/(2\sigma^2)$):

$\gamma$	$\sigma$	train accuracy	test accuracy	support vectors
0.01	7.07	0.843	0.827	198
0.1	2.24	0.883	0.850	125
1	0.71	0.937	0.880	79
10	0.22	0.940	0.883	126
100	0.07	0.990	0.857	256
1000	0.02	1.000	0.653	299

The linear SVM reaches $0.853$. With a very wide kernel the RBF SVM is no better than the line; around $\sigma \approx 0.2$–$0.7$ it gains three points; with a very narrow kernel it memorizes the training set (accuracy 1, 299 of 300 points are support vectors) and drops to $0.653$ on the test set. The number of support vectors is lowest where the test accuracy is highest. $\sigma$ and $C$ are tuned together by cross-validation, on logarithmic grids, after standardizing the features so that one $\sigma$ suits all of them.

The representer theorem#

Why does the solution of the SVM have the form $\sum_i \alpha_i y_i \phi(x_i)$, with one coefficient per training point, even when $\phi$ is infinite-dimensional? The reason is not specific to the hinge loss. Consider any regularized problem

$$\min_{w}\ \sum_{i=1}^n L\big(y_i,\ \langle w, \phi(x_i)\rangle\big) + \lambda \lVert w\rVert^2 ,\qquad \lambda \gt 0 .$$

Split $w = w_\parallel + w_\perp$, where $w_\parallel$ lies in the span of $\phi(x_1), \dots, \phi(x_n)$ and $w_\perp$ is orthogonal to all of them. The predictions on the training points depend only on $w_\parallel$, because $\langle w_\perp, \phi(x_i)\rangle = 0$, while $\lVert w\rVert^2 = \lVert w_\parallel\rVert^2 + \lVert w_\perp\rVert^2$. Dropping $w_\perp$ keeps the loss and strictly lowers the penalty when $w_\perp \ne 0$, so every minimizer has $w_\perp = 0$:

$$w = \sum_{i=1}^n \alpha_i \phi(x_i), \qquad f(x) = \langle w, \phi(x)\rangle = \sum_{i=1}^n \alpha_i\, k(x_i, x) .$$

This is the representer theorem (Kimeldorf and Wahba, 1971; Schölkopf, Herbrich and Smola, 2001, who allow any strictly increasing function of $\lVert w\rVert$ as the penalty; Shalev-Shwartz and Ben-David, 2014, Theorem 16.1). It turns a search over an infinite-dimensional space into a search over $n$ coefficients, and it holds for every loss: squared, logistic, hinge. Substituting the representation, $\lVert w\rVert^2 = \sum_{i,j} \alpha_i \alpha_j k(x_i, x_j) = \alpha^\top K \alpha$ and the vector of training predictions is $K\alpha$.

Kernel ridge regression#

With the squared loss, the problem is ridge regression in feature space, $\min_w \lVert y - \Phi w\rVert^2 + \lambda\lVert w\rVert^2$, where the rows of $\Phi$ are $\phi(x_i)^\top$. Its solution from the regularization note can be rewritten with the identity $\Phi^\top(\Phi\Phi^\top + \lambda I) = (\Phi^\top\Phi + \lambda I)\Phi^\top$: multiplying it by $(\Phi^\top\Phi + \lambda I)^{-1}$ on the left and by $(\Phi\Phi^\top + \lambda I)^{-1}$ on the right, both invertible for $\lambda \gt 0$, gives

$$w = (\Phi^\top \Phi + \lambda I)^{-1}\Phi^\top y = \Phi^\top (\Phi\Phi^\top + \lambda I)^{-1} y .$$

The left form inverts a matrix of the size of the feature space, the right form one of size $n \times n$, and $\Phi\Phi^\top = K$. So $w = \sum_i \alpha_i \phi(x_i)$ with

$$\alpha = (K + \lambda I)^{-1} y, \qquad f(x) = k(x)^\top (K + \lambda I)^{-1} y ,$$

the kernel ridge regression predictor, exactly the form the representer theorem promised (Bishop, 2006, §6.1).

Worked example. Two training points $x = 0$ and $x = 1$ with targets $1$ and $2$, the RBF kernel with $\sigma = 1$ and $\lambda = 0.1$. The kernel between the points is $e^{-1/2} = 0.6065$, so

$$K + \lambda I = \begin{pmatrix} 1.1 & 0.6065 \\ 0.6065 & 1.1 \end{pmatrix}, \qquad \alpha = (K + \lambda I)^{-1}\begin{pmatrix} 1 \\ 2\end{pmatrix} = \begin{pmatrix} -0.134 \\ 1.892\end{pmatrix} .$$

At $x = 0.5$ both kernel values are $e^{-1/8} = 0.8825$ and $f(0.5) = 0.8825 \cdot (-0.134 + 1.892) = 1.551$. At $x = 3$ the kernel values are $e^{-4.5} = 0.011$ and $e^{-2} = 0.135$, and $f(3) = 0.255$: far from the data an RBF regressor returns to 0, its prior mean, whatever the trend. At the training points themselves the fit is $f(0) = 1.013$ and $f(1) = 1.811$, not the targets $1$ and $2$: with $\lambda \gt 0$ the fit smooths rather than interpolates. Exact interpolation needs $\lambda = 0$ and an invertible $K$. Either way it does not extrapolate.

Figure

The RBF bandwidth as model complexity

The bandwidth trades bias for variance here as well. The effective number of parameters of a ridge fit, $\mathrm{df} = \operatorname{tr}\big(K(K + \lambda I)^{-1}\big)$, counts how many directions the fit uses (Hastie, Tibshirani and Friedman, 2009, §3.4.1); it is near $n$ when the fit interpolates and small when it is stiff. For the 30 points of the figure:

$\sigma$	1	0.5	0.3	0.2	0.1	0.05	0.02	0.01
df	2.36	3.31	4.43	5.77	9.25	14.69	21.61	24.26
training error	0.231	0.098	0.039	0.031	0.024	0.019	0.006	0.004
mean squared distance to the true curve	0.181	0.051	0.007	0.015	0.031	0.060	0.116	0.197

The training error falls monotonically as the kernel narrows; the mean squared distance to the true curve has the U shape of the bias–variance decomposition, with its minimum near $\sigma = 0.3$.

The Bayesian reading. For a finite feature map, put a Gaussian prior $w \sim N(0, \tau^2 I)$ on the weights and Gaussian noise of variance $s^2$ on the targets. Then $F(x) = \langle w, \phi(x)\rangle$ is a Gaussian process with covariance $\tau^2 k(x, z)$. For the RBF kernel the feature space is infinite-dimensional, and an isotropic Gaussian prior on infinitely many weights is only a formal limit of this picture; the model is then stated directly: $F$ is a Gaussian process with covariance $\tau^2 k(x, z)$. In both cases the posterior mean of $F(x)$ given the data follows from the conditional of a multivariate Gaussian: $\tau^2 k(x)^\top(\tau^2 K + s^2 I)^{-1} y = k(x)^\top (K + \lambda I)^{-1} y$ with $\lambda = s^2/\tau^2$. Kernel ridge regression is the posterior mean of Gaussian process regression, just as ridge regression is the MAP estimate under a Gaussian prior (Rasmussen and Williams, 2006, §2.2 and §6.2). The Gaussian process adds what the ridge view lacks: a posterior variance, an error bar at every $x$.

The cost, and random Fourier features#

Kernel methods trade the dimension of the features for the number of examples. Computing the Gram matrix costs $O(n^2 d)$ time and $O(n^2)$ memory, solving $(K + \lambda I)\alpha = y$ costs $O(n^3)$, and SVM solvers scale between $O(n^2)$ and $O(n^3)$ in practice. Every prediction costs one kernel evaluation per support vector or training point. For $n = 100{,}000$ the Gram matrix alone has $10^{10}$ entries, 80 GB in double precision. Kernel methods therefore rule small and medium data and lose to linear models, trees and networks on large data.

The remedy of Rahimi and Recht (2007) is to go back to explicit features, but few and random. It rests on Bochner’s theorem: a continuous shift-invariant kernel $k(x, z) = k(x - z)$ with $k(0) = 1$ is positive semidefinite if and only if it is the characteristic function of a probability distribution $p$,

$$k(x - z) = \mathbb E_{\omega \sim p}\, e^{i\langle \omega, x - z\rangle} = \mathbb E_{\omega \sim p} \cos\langle \omega, x - z\rangle$$

(Rahimi and Recht, §3, Theorem 1; the second form holds because $k$ is real, so the imaginary part averages to zero). Characteristic functions appeared in the expectation note. For the RBF kernel, $p$ is the Gaussian $N(0, \sigma^{-2} I)$, whose characteristic function is $e^{-\lVert t\rVert^2/(2\sigma^2)}$. Now draw $\omega_1, \dots, \omega_D$ from $p$ and $b_1, \dots, b_D$ uniformly from $[0, 2\pi]$, and set

$$z(x) = \sqrt{2/D}\,\big(\cos(\langle\omega_1, x\rangle + b_1), \dots, \cos(\langle\omega_D, x\rangle + b_D)\big) .$$

By $2\cos A \cos B = \cos(A - B) + \cos(A + B)$, each term of $\langle z(x), z(x')\rangle$ is $\tfrac1D\big(\cos\langle\omega_j, x - x'\rangle + \cos(\langle\omega_j, x + x'\rangle + 2b_j)\big)$, and the second cosine averages to zero over the uniform $b_j$. So $\langle z(x), z(x')\rangle$ is an average of $D$ independent terms with mean $k(x - x')$: a Monte Carlo estimate of the kernel. Each term lies in $[-2, 2]$, and Hoeffding’s inequality gives $P(\lvert\langle z(x), z(x')\rangle - k(x, x')\rvert \ge \varepsilon) \le 2e^{-D\varepsilon^2/8}$ for each pair; Rahimi and Recht (Claim 1) extend the bound to all pairs in a bounded set at once, at the cost of a factor that grows with $d$.

The error shrinks like $1/\sqrt D$. On 300 points in $\mathbb R^5$ the largest error over the whole $300 \times 300$ Gram matrix is $1.12$, $0.34$, $0.10$ and $0.035$ for $D = 10$, $100$, $1{,}000$ and $10{,}000$, a factor of about $\sqrt{10}$ per step. Any linear method trained on $z(x)$ then approximates the kernel method, at a cost of $O(nDd)$ for the features and $O(nD^2)$ for a ridge fit, linear in $n$. On the moons a linear SVM on $D = 10$ random features already reaches the test accuracy of the exact RBF SVM with the same $\gamma = 1$, $0.88$.

Kernels compare distributions#

Optional depth: kernel mean embeddings and the maximum mean discrepancy, used by the generative models course. The rest of the note does not depend on it.

The mean of the features of a distribution, $\mu_P = \mathbb E_{X \sim P}\phi(X)$, summarizes $P$ as one point of the feature space, and the distance between two such points can again be computed with the kernel alone:

$$\lVert \mu_P - \mu_Q\rVert^2 = \mathbb E\, k(X, X') + \mathbb E\, k(Y, Y') - 2\,\mathbb E\, k(X, Y), \qquad X, X' \sim P,\ Y, Y' \sim Q .$$

This is the squared maximum mean discrepancy. With a polynomial kernel it compares moments up to the degree of the kernel; with the RBF kernel it is zero only when $P = Q$ (Gretton et al., 2012). The generative models course uses it to compare generated images with real ones: the Kernel Inception Distance is this quantity with a cubic polynomial kernel on network features.

What changed#

	Linear SVM and ridge regression	Kernel methods
Decision function	$\langle w, x\rangle + b$, a hyperplane	$\sum_i \alpha_i k(x_i, x) + b$, curved in $x$
Parameters	$d$ weights	$n$ coefficients, one per training point (representer theorem)
What the data enter as	feature vectors	the Gram matrix $K_{ij} = k(x_i, x_j)$
Feature space	$\mathbb R^d$	any space with an inner product, infinite for the RBF kernel
Complexity knob	$C$ or $\lambda$	$C$ or $\lambda$ plus the kernel and its bandwidth
Training cost	linear in $n$	$O(n^2)$ memory, $O(n^2)$–$O(n^3)$ time; random features bring it back to linear

Limitations#

• Quadratic memory, cubic time. Exact kernel methods do not scale past a few tens of thousands of examples; approximations (random features, the Nyström method) trade accuracy for speed.
• The kernel is a modelling choice. It encodes which inputs count as similar. The RBF kernel treats all directions alike after standardization; with many irrelevant features the distances are dominated by noise, as in the curse of dimensionality of the nearest neighbours note.
• Hyperparameters. The bandwidth and $C$ or $\lambda$ interact and need a joint search; the moons table shows test accuracy between $0.65$ and $0.88$ depending on $\sigma$ alone.
• No extrapolation with local kernels. An RBF fit returns to zero away from the data, as the worked example showed.
• Fixed features. The feature map is chosen before seeing the data. Neural networks learn the features instead, which is what lets them scale to images and text.

What you should remember#

• A kernel is an inner product after a feature map, $k(x, z) = \langle\phi(x), \phi(z)\rangle$. Any method written with inner products, such as the SVM dual, works with $k$ in place of $\langle x, z\rangle$ without forming $\phi$.
• $k$ is a valid kernel if and only if all its Gram matrices are symmetric positive semidefinite: equivalently, $k$ is a covariance function. Sums, positive multiples, products, limits and $g(x)g(z)$ preserve this; distances and $\tanh$ do not have it.
• Polynomial kernel: all monomials up to degree $p$. RBF kernel: infinite-dimensional features, similarity on the scale $\sigma$; large $\sigma$ underfits, small $\sigma$ overfits.
• Representer theorem: for any loss plus $\lambda\lVert w\rVert^2$, the minimizer is $\sum_i \alpha_i \phi(x_i)$. Kernel ridge regression: $\alpha = (K + \lambda I)^{-1}y$, $f(x) = k(x)^\top \alpha$, the posterior mean of a Gaussian process.
• Exact kernel methods cost $O(n^2)$ memory and up to $O(n^3)$ time. Random Fourier features estimate a shift-invariant kernel by Monte Carlo with error $O(1/\sqrt D)$ and make training linear in $n$.

Code#

A polynomial kernel as an inner product, Gram matrices that are and are not positive semidefinite, kernel ridge regression by hand, the RBF bandwidth on the moons, and random Fourier features:

import numpy as np, math
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC, LinearSVC
from sklearn.kernel_ridge import KernelRidge
from sklearn.kernel_approximation import RBFSampler

# a kernel is an inner product after a feature map: (x.z)^2 with phi(x) = (x1^2, sqrt(2) x1 x2, x2^2)
x, z = np.array([1.0, 2.0]), np.array([3.0, 1.0])
phi = lambda v: np.array([v[0]**2, math.sqrt(2) * v[0] * v[1], v[1]**2])
print((x @ z)**2, phi(x) @ phi(z), math.comb(100 + 3, 3))   # 25.0 25.0 176851: monomials of degree <= 3 in 100 features

# which functions are kernels: the smallest eigenvalue of the Gram matrix on 50 random points
rng = np.random.default_rng(0)
P = rng.normal(size=(50, 2)); D2 = ((P[:, None] - P[None]) ** 2).sum(-1)
for name, G in (('rbf', np.exp(-D2 / 2)), ('distance', np.sqrt(D2)), ('tanh', np.tanh(P @ P.T + 1))):
    print(name, f'{np.linalg.eigvalsh(G).min():.3g}')
# rbf 4.24e-10 | distance -27.5 | tanh -6.26: only the RBF Gram matrix is positive semidefinite

# kernel ridge regression by hand: x = (0, 1), y = (1, 2), RBF with sigma = 1, lambda = 0.1
X = np.array([[0.0], [1.0]]); y = np.array([1.0, 2.0])
k = lambda A, B: np.exp(-(A - B.T) ** 2 / 2)
alpha = np.linalg.solve(k(X, X) + 0.1 * np.eye(2), y)
Xn = np.array([[0.5], [3.0]])
print(alpha.round(4), (k(Xn, X) @ alpha).round(4))     # [-0.1343  1.8922] [1.5514 0.2546]
print((k(X, X) @ alpha).round(3))                      # [1.013 1.811]: the fit at the training points, not (1, 2)
print(KernelRidge(alpha=0.1, kernel='rbf', gamma=0.5).fit(X, y).predict(Xn).round(4))   # [1.5514 0.2546]

# the RBF bandwidth as bias-variance: an SVM on two noisy moons, gamma = 1/(2 sigma^2)
Xm, ym = make_moons(n_samples=600, noise=0.3, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(Xm, ym, test_size=0.5, random_state=0)
print('linear', round(SVC(kernel='linear').fit(Xtr, ytr).score(Xte, yte), 3))   # linear 0.853
for g in (0.01, 0.1, 1, 10, 100, 1000):
    m = SVC(kernel='rbf', gamma=g, C=1).fit(Xtr, ytr)
    print(g, round(1 / math.sqrt(2 * g), 3), round(m.score(Xtr, ytr), 3), round(m.score(Xte, yte), 3), len(m.support_))
# gamma, sigma, train and test accuracy, support vectors
# 0.01 7.071 0.843 0.827 198 | 0.1 2.236 0.883 0.85 125 | 1 0.707 0.937 0.88 79
# 10 0.224 0.94 0.883 126    | 100 0.071 0.99 0.857 256 | 1000 0.022 1.0 0.653 299

# random Fourier features: the largest error over a 300 x 300 Gram matrix (5 features, gamma = 0.1), mean of 5 draws
Q = rng.normal(size=(300, 5)); K = np.exp(-0.1 * ((Q[:, None] - Q[None]) ** 2).sum(-1))
for D in (10, 100, 1000, 10000):
    err = [np.abs((Z := RBFSampler(gamma=0.1, n_components=D, random_state=s).fit_transform(Q)) @ Z.T - K).max() for s in range(5)]
    print(D, round(float(np.mean(err)), 3))
# 10 1.118 | 100 0.343 | 1000 0.101 | 10000 0.035: about 1/sqrt(D)

# a linear SVM on D random features approaches the RBF SVM (gamma = 1: test accuracy 0.88)
for D in (2, 10, 100, 1000):
    f = RBFSampler(gamma=1, n_components=D, random_state=0).fit(Xtr)
    print(D, round(LinearSVC(C=1).fit(f.transform(Xtr), ytr).score(f.transform(Xte), yte), 3))
# 2 0.663 | 10 0.883 | 100 0.88 | 1000 0.877

Next#

The SVM and its kernels separate two classes. Many problems have more: ten digits, thousands of topics, a vocabulary of tens of thousands of words, or several labels at once for the same object. The next note reduces such problems to binary ones, extends logistic regression to many classes with the softmax, and adapts the metrics of the classification metrics note to several classes and to several labels per example.

References#

• B. E. Boser, I. M. Guyon, V. N. Vapnik. A Training Algorithm for Optimal Margin Classifiers. COLT, 1992, pp. 144–152. The kernel trick applied to the maximum-margin classifier.
• N. Aronszajn. Theory of Reproducing Kernels. Transactions of the American Mathematical Society, 68(3), pp. 337–404, 1950.
• G. Kimeldorf, G. Wahba. Some Results on Tchebycheffian Spline Functions. Journal of Mathematical Analysis and Applications, 33(1), pp. 82–95, 1971. The representer theorem for the squared loss.
• B. Schölkopf, R. Herbrich, A. J. Smola. A Generalized Representer Theorem. COLT, LNCS 2111, pp. 416–426, 2001.
• A. Rahimi, B. Recht. Random Features for Large-Scale Kernel Machines. NIPS, 2007. §3: Bochner’s theorem (Theorem 1), random Fourier features (Algorithm 1), uniform convergence (Claim 1).
• A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, A. Smola. A kernel two-sample test. JMLR, 13, pp. 723–773, 2012. The maximum mean discrepancy.
• C. E. Rasmussen, C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006. §2.1–2.2 (weight-space and function-space views of regression), §4.3 (eigenfunction analysis of kernels, Mercer’s theorem), §6.1–6.2 (reproducing kernel Hilbert spaces, regularization).
• S. Shalev-Shwartz, S. Ben-David. Understanding Machine Learning: From Theory to Algorithms. Cambridge University Press, 2014. Chapter 16: embeddings into feature spaces, the kernel trick, Theorem 16.1 (representer theorem), Lemma 16.2 (positive semidefinite kernels).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. §6.1 (dual representations), §6.2 (constructing kernels), §6.4 (Gaussian processes).
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. §3.4.1 (ridge regression and its effective degrees of freedom), §5.8 (regularization and reproducing kernel Hilbert spaces), §12.3 (support vector machines and kernels).
• E. Sokolov et al. Машинное обучение (HSE), lecture notes: lecture 13 (kernel methods, constructing kernels, the polynomial and Gaussian kernels), lecture 14 (the kernel SVM, random Fourier features).