Kernel Methods & Reproducing Kernel Hilbert Space
Kernel methods provide a principled framework for extending linear algorithms to nonlinear settings by implicitly mapping data to high-dimensional feature spaces. The mathematical foundation is the theory of reproducing kernel Hilbert spaces (RKHS).
1. Positive Definite Kernels
1.1 Definition
Definition: A function is a positive definite kernel if for all , all , and all :
Equivalently, the kernel matrix with is positive semi-definite for all finite sets of points.
1.2 Mercer's Theorem
Theorem (Mercer, 1909): If is continuous, symmetric, and positive definite on a compact set , then there exist an orthonormal basis of and non-negative eigenvalues such that:
The convergence is absolute and uniform. This decomposition defines the feature map:
such that .
1.3 Common Kernels
RBF (Gaussian) Kernel:
This corresponds to an infinite-dimensional feature space. The eigenfunctions are Hermite polynomials.
Polynomial Kernel:
The feature space consists of all monomials of degree up to .
Matérn Kernel:
where is the modified Bessel function. Controls smoothness via .
2. Reproducing Kernel Hilbert Space (RKHS)
2.1 Definition
Definition: A Hilbert space of functions is a reproducing kernel Hilbert space if there exists a kernel such that:
- For all , (containment)
- For all and : (reproducing property)
2.2 Moore-Aronszajn Theorem
Theorem: For every positive definite kernel , there exists a unique RKHS with as its reproducing kernel.
Construction: Define the feature map . The RKHS is:
with inner product defined by .
2.3 Representer Theorem
Theorem (Kimeldorf & Wahba, 1970): Consider the optimization problem:
where is any loss function and . Then the optimal has the form:
for some .
Significance: The solution lies in a finite-dimensional subspace spanned by kernel evaluations at training points, regardless of the (possibly infinite) dimension of .
3. Kernel Ridge Regression
3.1 Formulation
Given data , solve:
By the representer theorem, . Substituting:
where is the kernel matrix.
3.2 Closed-Form Solution
Taking derivatives and setting to zero:
The prediction at a new point :
where .
3.3 Computational Considerations
- Naive: for matrix inversion
- Cholesky: but numerically stable;
- Nyström approximation: where is the number of inducing points
- Random Fourier features: where is the number of random features
4. Support Vector Machines
4.1 Maximum Margin Classification
For binary classification with labels , the SVM solves:
subject to:
4.2 Dual Formulation
Using Lagrangian duality:
subject to:
The decision function:
4.3 Support Vectors
Points with are support vectors. For the hard margin case (), support vectors lie exactly on the margin:
The margin width is .
4.4 Kernel Trick
The SVM dual only involves inner products . This is the kernel trick: we can work in infinite-dimensional feature spaces without explicitly computing .
Example: The RBF kernel corresponds to an infinite-dimensional feature space, but computing the kernel only requires operations.
5. Kernel PCA
5.1 Formulation
Standard PCA finds directions of maximum variance in the original feature space. Kernel PCA performs PCA in the RKHS .
The centered kernel matrix is:
Eigendecomposition:
5.2 Projection
The projection of a new point onto the -th principal component:
5.3 Feature Space Interpretation
In the feature space, the principal components are:
The projection is:
6. Gaussian Processes as Kernel Methods
6.1 Definition
A Gaussian process is a collection of random variables such that any finite subset has a multivariate Gaussian distribution:
where and .
6.2 Connection to RKHS
The RKHS norm of a GP sample:
where are the expansion coefficients in the Mercer basis. This connects the measure-theoretic GP to the function-theoretic RKHS.
6.3 Posterior Distribution
Given training data with , :
where:
with .
7. Kernel Methods Theory
7.1 Kernel Matrix Properties
Theorem: For a positive definite kernel and points :
- is symmetric positive semi-definite
- All eigenvalues
- (with equality if is strictly positive definite)
7.2 Universal Kernels
Definition: A kernel is universal on if the RKHS is dense in (continuous functions with sup norm).
Examples:
- RBF kernel is universal on any compact
- Polynomial kernels are not universal (only approximate polynomials)
7.3 Approximation Bounds
Theorem: For the RBF kernel with bandwidth and training points, any function can be approximated by:
where is the kernel ridge regression estimator.
8. String and Graph Kernels
8.1 String Kernels
The spectrum kernel for strings :
where counts occurrences of substring in .
8.2 Graph Kernels
The random walk kernel:
where are adjacency matrices and is the Kronecker product.
9. Code: Kernel Methods Implementation
import numpy as np
from scipy.optimize import minimize
from scipy.linalg import cho_solve, cho_factor
class KernelRidgeRegression:
"""Kernel Ridge Regression with various kernels."""
def __init__(self, kernel='rbf', gamma=1.0, degree=3, coef0=1.0, lambda_reg=1.0):
self.kernel = kernel
self.gamma = gamma
self.degree = degree
self.coef0 = coef0
self.lambda_reg = lambda_reg
self.alpha_ = None
self.X_train_ = None
self.y_train_ = None
def _compute_kernel(self, X1, X2):
"""Compute kernel matrix between X1 and X2."""
if self.kernel == 'rbf':
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return np.exp(-self.gamma * sq_dists)
elif self.kernel == 'polynomial':
return (self.gamma * X1 @ X2.T + self.coef0) ** self.degree
elif self.kernel == 'linear':
return X1 @ X2.T
elif self.kernel == 'sigmoid':
return np.tanh(self.gamma * X1 @ X2.T + self.coef0)
else:
raise ValueError(f"Unknown kernel: {self.kernel}")
def fit(self, X, y):
"""Fit kernel ridge regression."""
self.X_train_ = X
self.y_train_ = y
K = self._compute_kernel(X, X)
# Add regularization to diagonal
K_reg = K + self.lambda_reg * np.eye(len(X))
# Solve using Cholesky decomposition
L, low = cho_factor(K_reg)
self.alpha_ = cho_solve((L, low), y)
return self
def predict(self, X):
"""Predict at new points."""
K = self._compute_kernel(X, self.X_train_)
return K @ self.alpha_
def score(self, X, y):
"""Compute R^2 score."""
y_pred = self.predict(X)
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
return 1 - ss_res / ss_tot
class SupportVectorMachine:
"""Support Vector Machine with kernel trick."""
def __init__(self, kernel='rbf', C=1.0, gamma=1.0, degree=3):
self.kernel = kernel
self.C = C
self.gamma = gamma
self.degree = degree
self.alpha_ = None
self.b_ = None
self.X_train_ = None
self.y_train_ = None
self.support_vectors_ = None
def _compute_kernel(self, X1, X2):
if self.kernel == 'rbf':
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return np.exp(-self.gamma * sq_dists)
elif self.kernel == 'polynomial':
return (X1 @ X2.T + 1) ** self.degree
else:
return X1 @ X2.T
def fit(self, X, y):
"""Fit SVM using SMO algorithm (simplified)."""
n = len(X)
self.X_train_ = X
self.y_train_ = y
K = self._compute_kernel(X, X)
# Simplified: use cvxopt or scipy if available
# For demonstration, use a simple QP solver approach
from scipy.optimize import minimize
def objective(alpha):
return 0.5 * alpha @ (np.outer(y, y) * K) @ alpha - np.sum(alpha)
constraints = [{'type': 'eq', 'fun': lambda a: np.sum(a * y)}]
bounds = [(0, self.C) for _ in range(n)]
alpha0 = np.zeros(n)
result = minimize(objective, alpha0, method='SLSQP',
bounds=bounds, constraints=constraints)
self.alpha_ = result.x
# Find support vectors
sv_mask = self.alpha_ > 1e-6
self.support_vectors_ = X[sv_mask]
# Compute bias
sv_alpha = self.alpha_[sv_mask]
sv_y = y[sv_mask]
self.b_ = np.mean(sv_y - sv_alpha @ K[sv_mask][:, sv_mask])
return self
def predict(self, X):
K = self._compute_kernel(X, self.X_train_)
return np.sign(K @ (self.alpha_ * self.y_train_) + self.b_)
class KernelPCA:
"""Kernel Principal Component Analysis."""
def __init__(self, kernel='rbf', gamma=1.0, n_components=2):
self.kernel = kernel
self.gamma = gamma
self.n_components = n_components
self.eigenvectors_ = None
self.eigenvalues_ = None
self.X_train_ = None
self.mean_K_ = None
def _compute_kernel(self, X1, X2):
if self.kernel == 'rbf':
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return np.exp(-self.gamma * sq_dists)
return X1 @ X2.T
def fit(self, X):
"""Fit kernel PCA."""
self.X_train_ = X
n = len(X)
# Compute centered kernel matrix
K = self._compute_kernel(X, X)
# Center in feature space
one_n = np.ones((n, n)) / n
K_centered = K - one_n @ K - K @ one_n + one_n @ K @ one_n
# Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(K_centered)
# Sort and select top components
idx = np.argsort(eigenvalues)[::-1]
self.eigenvalues_ = eigenvalues[idx[:self.n_components]]
self.eigenvectors_ = eigenvectors[:, idx[:self.n_components]]
return self
def transform(self, X):
"""Transform new data."""
K = self._compute_kernel(X, self.X_train_)
# Center test kernel matrix
n = len(self.X_train_)
one_n = np.ones((n, n)) / n
K_centered = K - one_n @ K.T # Approximate centering
return K_centered @ self.eigenvectors_ / np.sqrt(self.eigenvalues_)
class GaussianProcess:
"""Gaussian Process Regression."""
def __init__(self, kernel='rbf', noise=1e-2, length_scale=1.0, signal_variance=1.0):
self.kernel = kernel
self.noise = noise
self.length_scale = length_scale
self.signal_variance = signal_variance
self.X_train_ = None
self.y_train_ = None
self.L_ = None
self.alpha_ = None
def _compute_kernel(self, X1, X2):
if self.kernel == 'rbf':
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return self.signal_variance * np.exp(-0.5 * sq_dists / self.length_scale**2)
return X1 @ X2.T
def fit(self, X, y):
"""Fit GP by computing Cholesky decomposition."""
self.X_train_ = X
self.y_train_ = y
K = self._compute_kernel(X, X)
K += self.noise * np.eye(len(X))
# Cholesky decomposition
self.L_ = np.linalg.cholesky(K)
self.alpha_ = np.linalg.solve(self.L_.T, np.linalg.solve(self.L_, y))
return self
def predict(self, X, return_std=False):
"""Predict with uncertainty."""
K_s = self._compute_kernel(X, self.X_train_)
K_ss = self._compute_kernel(X, X)
mu = K_s @ self.alpha_
if return_std:
v = np.linalg.solve(self.L_, K_s.T)
var = np.diag(K_ss) - np.sum(v**2, axis=0)
var = np.maximum(var, 0) # Numerical stability
return mu, np.sqrt(var)
return mu
def log_marginal_likelihood(self):
"""Compute log marginal likelihood."""
y = self.y_train_
alpha = self.alpha_
# -0.5 * y^T * alpha - sum(log(diag(L))) - n/2 * log(2*pi)
lml = -0.5 * y @ alpha - np.sum(np.log(np.diag(self.L_))) - len(y) / 2 * np.log(2 * np.pi)
return lml
def kernel_alignment(K1, K2):
"""Compute kernel alignment between two kernel matrices."""
# Frobenius inner product
frobenius = np.sum(K1 * K2)
# Norms
norm1 = np.sqrt(np.sum(K1 * K1))
norm2 = np.sqrt(np.sum(K2 * K2))
return frobenius / (norm1 * norm2)
def kernel_target_alignment(X, y, kernel_func):
"""
Compute kernel-target alignment.
Measures how well the kernel matches the classification task.
"""
K = kernel_func(X, X)
# Target kernel (outer product of labels)
K_target = np.outer(y, y)
return kernel_alignment(K, K_target)
10. Summary
Kernel methods provide a unified framework for nonlinear learning:
- Mercer's theorem decomposes kernels into inner products in feature space
- RKHS theory provides the function space foundation
- Representer theorem ensures finite-dimensional solutions
- Kernel trick enables infinite-dimensional computation in finite time
- Gaussian processes connect kernel methods to Bayesian inference
The key insight is that the kernel function implicitly computes inner products in potentially infinite-dimensional spaces, enabling nonlinear extensions of linear algorithms while maintaining computational tractability.