Gaussian Processes and Kernel Methods Rbf Svm Bayesian Inference
# Gaussian Processes and Kernel Methods: RBF, SVM, Bayesian Inference
## 1. Introduction & Motivation
Gaussian Processes (GPs) and kernel methods provide powerful, theoretically principled approaches to regression and classification. Unlike deep learning's black-box predictions, GPs provide:
- Uncertainty quantification: Predict not just point estimates but confidence intervals
- Sample efficiency: Learn from limited data
- Interpretability: Kernel structure reveals data relationships
- Bayesian inference: Principled uncertainty via posterior distributions
Kernel methods including SVMs, GPs, and kernel ridge regression share a mathematical foundation enabling efficient computation in high-dimensional spaces through the "kernel trick."
This article comprehensively covers kernel methods, GP theory, practical implementation, and applications.
## 2. Core Concepts & Theory
### 2.1 Gaussian Processes
A GP is a probability distribution over functions, fully specified by mean and covariance:
$$f(x) \sim \mathcal{GP}(m(x), k(x, x'))$$
Prior beliefs about function smoothness encoded in kernel k.
For regression:
$$y = f(x) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma_n^2)$$
### 2.2 Kernel Functions
Kernel k(x, x') measures similarity between inputs. Common kernels:
RBF (Radial Basis Function):
$$k_{ ext{RBF}}(x, x') = \sigma^2 \exp\left(-\frac{\|x - x'\|^2}{2\ell^2} ight)$$
where
$$ \ell $$
is length scale,
$$ \sigma^2 $$
is variance. Smooth, flexible.
Matérn:
$$k_{ ext{Matern}}(x, x') = \frac{\sigma^2}{2^{ u-1}\Gamma( u)} \left(\sqrt{2 u} \frac{\|x-x'\|}{\ell} ight)^ u K_ u\left(\sqrt{2 u} \frac{\|x-x'\|}{\ell} ight)$$
Generalization of RBF, controls smoothness via
$$ u $$
.
Linear:
$$k_{ ext{linear}}(x, x') = \sigma^2 x^T x'$$
For linear functions.
Periodic:
$$k_{ ext{periodic}}(x, x') = \sigma^2 \exp\left(-\frac{2\sin^2(\pi(x-x')/p)}{\ell^2} ight)$$
For seasonal/periodic patterns.
### 2.3 Posterior Distribution
Given observations (X, y), posterior over functions:
$$p(f | X, y) = \mathcal{N}(f; \mu, \Sigma)$$
Mean (point estimate):
$$\mu = K(X, X)(K(X, X) + \sigma_n^2 I)^{-1} y$$
Covariance (uncertainty):
$$\Sigma = K(X, X) - K(X, X)(K(X, X) + \sigma_n^2 I)^{-1} K(X, X)$$
### 2.4 Support Vector Machines (SVM)
Maximum margin classifier for binary classification:
$$\min_{w, b, \xi} \frac{1}{2}\|w\|^2 + C \sum_i \xi_i$$
$$ ext{s.t.} \quad y_i(w^T \phi(x_i) + b) \geq 1 - \xi_i$$
Hinge loss with L2 regularization. Kernel trick maps to high-dimensional space without explicit
$$ \phi $$
.
## 3. Mathematical Formulation
### 3.1 Kernel Trick
Avoid expensive explicit feature mapping
$$ \phi(x) $$
:
$$K(x, x') = \phi(x)^T \phi(x')$$
Compute inner product in feature space efficiently in original space.
Example: Polynomial kernel
$$k(x, x') = (x^T x' + c)^d = \langle \phi(x), \phi(x')
angle$$
where
$$ \phi $$
implicitly computes all degree-d monomials (exponentially many).
### 3.2 GP Posterior Predictive Distribution
For new point
$$ x_* $$
:
$$p(f_* | X, y, x_*) = \mathcal{N}(\mu_*, \sigma_*^2)$$
Mean:
$$\mu_* = k_*^T K^{-1} y$$
Variance:
$$\sigma_*^2 = k(x_*, x_*) - k_*^T K^{-1} k_*$$
where
$$ k_* = k(x_*, X) $$
(vector of covariances).
### 3.3 SVM Dual Formulation
Dual problem:
$$\max_\alpha \sum_i \alpha_i - \frac{1}{2} \sum_{i,j} \alpha_i \alpha_j y_i y_j k(x_i, x_j)$$
$$ ext{s.t.} \quad 0 \leq \alpha_i \leq C, \quad \sum_i \alpha_i y_i = 0$$
Decision function:
$$f(x) = ext{sign}\left(\sum_i \alpha_i y_i k(x, x_i) + b ight)$$
Only
$$ \alpha_i $$
corresponding to support vectors non-zero.
### 3.4 Kernel Matrix Properties
For valid kernel (Mercer's theorem), kernel matrix
$$ K = [k(x_i, x_j)] $$
must be:
- Symmetric:
$$ K = K^T $$
- Positive semi-definite:
$$ z^T K z \geq 0 $$
for all z
Ensures well-defined feature space.
## 4. Advanced Theory & Extensions
### 4.1 Hyperparameter Optimization
GP hyperparameters: kernel variance
$$ \sigma^2 $$
, length scales
$$ \ell $$
, noise
$$ \sigma_n^2 $$
.
Marginal likelihood (evidence):
$$\log p(y | X) = -\frac{1}{2} y^T K^{-1} y - \frac{1}{2} \log |K| - \frac{n}{2} \log(2\pi)$$
Optimize via gradient descent or Bayesian optimization.
### 4.2 Sparse Approximations
Full GP requires
$$ O(n^3) $$
for matrix inversion. For
$$ n > 10K $$
, expensive.
Sparse GPs use inducing points
$$ Z \subset X $$
:
$$\approx ext{GP}_{ ext{sparse}}(m(x), k(x, x') - k(x, Z)K_z^{-1}k(Z, x'))$$
Reduces complexity to
$$ O(m^3) $$
where
$$ m << n $$
.
### 4.3 Multi-Output and Multi-Task GPs
Extend to vector outputs:
$$\mathbf{f}(x) = [f_1(x), \ldots, f_d(x)]^T \sim \mathcal{GP}$$
Shared structure across tasks improves sample efficiency.
### 4.4 Deep Gaussian Processes
Stack multiple GP layers:
$$f(x) = ext{GP}_2( ext{GP}_1(x))$$
More expressive but harder to train (non-conjugate inference).
## 5. Computational Considerations
### 5.1 Complexity Analysis
GP inference:
- Cholesky decomposition:
$$ O(n^3) $$
- Prediction:
$$ O(n^2) $$
per point
- Hyperparameter gradient:
$$ O(n^3) $$
Infeasible for
$$ n > 10K $$
without approximations.
SVM training:
- Quadratic programming:
$$ O(n^2) $$
to
$$ O(n^3) $$
- Prediction: O(m) where m is support vectors (often
$$ m << n $$
)
Practical for
$$ n < 1M $$
with efficient solvers.
### 5.2 Kernel Computation Optimization
Vectorized computation:
- Compute all pairwise distances efficiently
- Use NumPy broadcasting
- GPU acceleration for large datasets
Caching:
- Store kernel matrix (if n small)
- Avoid recomputing identical pairs
### 5.3 Memory Management
Full kernel matrix:
$$ O(n^2) $$
memory
- For
$$ n = 10K $$
: ~800MB (manageable)
- For
$$ n = 100K $$
: ~80GB (infeasible)
Sparse approximations:
$$ O(n \cdot m) $$
memory
- Much more practical
### 5.4 Parallelization
Embarrassingly parallel:
- Kernel computation across data pairs
- Multiple hyperparameter evaluations
GPU acceleration:
- Matrix operations
- Cholesky decomposition
- 10-100x speedup typical
## 6. Practical Implementation Strategies
### 6.1 Choosing a Kernel
RBF: Default choice
- Smooth, flexible
- Works well for most problems
- Interpretable length scale
Matérn: More control over smoothness
-
$$ u = 1/2 $$
: Non-smooth
-
$$ u = \infty $$
: Infinitely smooth (RBF)
- Middle values: Balances smoothness
Linear: For high-dimensional sparse data
- Interpretable weights
- Fast computation
Periodic: For temporal/seasonal data
- Explicit periodicity
- Better for forecasting
Additive: Combine multiple kernels
-
$$ k = k_1 + k_2 + \cdots $$
- Different structure for different features
### 6.2 Hyperparameter Selection
**Length scale
$$ \ell $$
:**
- Controls how far points influence each other
- Smaller
$$ \ell $$
: Wiggly (overfitting)
- Larger
$$ \ell $$
: Smooth (underfitting)
- Practical: Learn via marginal likelihood maximization
**Signal variance
$$ \sigma^2 $$
:**
- Scale of function
- Usually learned from data
**Noise variance
$$ \sigma_n^2 $$
:**
- Measurement noise
- Often hyperparameter
- Small values: Trust data, risk overfitting
- Large values: Ignore noise, risk underfitting
### 6.3 Data Normalization
Essential for kernel methods:
$$x' = \frac{x - \mu}{\sigma}$$
Ensures all features on similar scale. Length scale then meaningful.
### 6.4 Efficient GP Training
Initialization: Start from reasonable hyperparameters
- Use empirical length scale estimate
- Data-dependent initialization
Optimization: Use L-BFGS or SGD
- Marginal likelihood non-convex; local optima acceptable
- Restarts help escape poor local minima
Checkpointing: Store intermediate results
- Restart training if interrupted
## 7. Benchmark Datasets & Evaluation
### 7.1 Regression Benchmarks
Boston Housing:
- 506 samples, 13 features
- GP RBF RMSE: 3.1-3.5
- Neural network: 3.0-3.5
- Competitive performance, better uncertainty
California Housing:
- 20,640 samples, 8 features
- GP RBF (subset of data): RMSE 0.73
- Neural network: RMSE 0.71
- Similar performance, GP slower due to n
### 7.2 Classification Benchmarks
UCI Datasets (small to medium):
- Breast cancer: GP ~97% accuracy, SVM ~96%
- Iris: GP ~97%, SVM ~97%
- Typically comparable
Large datasets (n > 10K):
- Sparse GP: ~95% accuracy
- Neural network: ~96%+
- Deep learning faster, accuracy slightly better
### 7.3 Uncertainty Quantification
Calibration: Do predicted intervals contain true values at stated rate?
Proper calibration: 68% of points in ±1σ interval, 95% in ±2σ interval.
Results:
- GP: ~95% calibrated
- Neural network ensembles: 60-80% (overconfident)
- Modern Bayesian NN: 85-90%
### 7.4 Evaluation Metrics
Negative log likelihood (NLL):
$$ ext{NLL} = -\log p(y | x, ext{model})$$
Rewards both accuracy and calibration. GP typically ~0.2-0.5 NLL better than non-Bayesian.
## 8. Key Challenges & Limitations
### 8.1 Computational Bottleneck
$$ O(n^3) $$
complexity limits practical scale:
- Sparse approximations necessary for
$$ n > 10K $$
- Approximations reduce predictive quality
- Trade-off between speed and accuracy
### 8.2 Hyperparameter Sensitivity
GP performance sensitive to hyperparameters:
- Wrong length scale: 10-30% accuracy loss
- Non-smooth kernel on smooth function: Slow convergence
- Optimization gets stuck in local optima
Requires careful tuning or Bayesian optimization.
### 8.3 High-Dimensional Curse
RBF kernel struggles in high dimensions:
- Length scale uninformative (all points equidistant)
- Requires exponentially more data
- Feature selection or dimensionality reduction needed
### 8.4 Limited Scalability
- Full GPs: Up to ~10K points practical
- Sparse GPs: Up to ~1M points
- Neural networks: Billions of samples
For big data, neural networks often preferred.
## 9. Hyperparameter Tuning & Optimization
### 9.1 Kernel Hyperparameters
RBF length scale: 0.1-10 typical
- Too small: Noisy fit
- Too large: Underfitting
- Grid search or Bayesian optimization
Variance: Often
$$ \sigma^2 = ext{Var}(y) $$
(data-dependent)
Noise: Depends on measurement error
- Typical: 0.01 to 0.1 times data variance
- Cross-validation if uncertain
### 9.2 Sparse GP Hyperparameters
Number of inducing points: 100-1000 typical
- More: Better accuracy, higher cost
- Typical: min(1000, n/10)
Placement: Grid or learned
- Learned via variational optimization
- Often better than grid
### 9.3 SVM Parameters
Regularization C: Balances fit and margin
- Smaller C: Wider margin, may underfit
- Larger C: Tight margin, risk overfitting
- Typical: 0.1-10, tune via cross-validation
Kernel parameters: Depend on kernel choice
- RBF gamma: 0.001-1, typically 1/n_features
- Polynomial degree: 2-4
### 9.4 Optimization Strategies
Marginal likelihood maximization:
- Gradient-based (L-BFGS)
- Multiple restarts (helps escape local optima)
- Typical: 10-100 iterations
Cross-validation:
- For hyperparameter selection
- Nested CV for robust estimate
- Computationally expensive for GPs
## 10. Real-World Applications & Case Studies
### 10.1 Bayesian Optimization
Problem: Find minimum of expensive black-box function
Setup:
- Expensive function evaluation (e.g., experiment, simulation)
- Limited budget (10-100 evaluations)
- Multiple dimensions (5-20)
Method:
- Surrogate: GP model of function
- Acquisition function: Balances exploration/exploitation
- Iteratively evaluate promising points
Results:
- Finds good solutions in ~50 evaluations
- Neural networks would need 1000s
- Standard methods (grid search) very inefficient
Example: Hyperparameter tuning, drug discovery, engineering design
### 10.2 Time Series Forecasting with GPs
Problem: Forecast electricity demand with uncertainty
Setup:
- Data: Daily demand for 365 days
- Goal: Forecast 7 days ahead with confidence intervals
- Kernel: Combination of periodic + linear trend
Method:
- Kernel:
$$ k = \sigma^2_{p} ext{exp}(-\sin^2(\pi(t-t')/365)/\ell_p^2) + \sigma^2_{l}(1 + t \cdot t') $$
- Combines periodicity (daily) and linear trend
- Marginal likelihood optimization for hyperparameters
Results:
- Point accuracy: RMSE 3.2 (comparable to LSTM)
- Uncertainty: Confidence intervals well-calibrated
- Interpretability: Trend and seasonality separated
Advantage over LSTM:
- Better calibrated uncertainty
- Requires less data
- Interpretable components
### 10.3 Active Learning with GP
Problem: Learn model with minimal labels (expensive annotation)
Setup:
- Unlabeled pool: 10K samples
- Budget: 500 labels
- Goal: Maximize accuracy with minimum labels
Method:
- Initialize: Train GP on small random sample
- Iteratively: Select points of highest predicted uncertainty
- Label and retrain
Results:
- Random sampling: 85% accuracy with 500 labels
- Uncertainty sampling (GP): 92% accuracy with 500 labels
- Efficiency: 60% fewer labels needed
Why GP good:
- Principled uncertainty
- Exploration-exploitation balance
- Sample-efficient learning
### 10.4 Gaussian Process Regression for Calibration
Problem: Calibrate sensor measurements
Setup:
- Sensor: Biased, noisy measurements
- Reference: Expensive true values (only 50 points)
- Goal: Use 50 reference measurements to calibrate 10K sensor readings
Method:
- Train GP on (sensor, true) pairs
- Posterior predictive: Calibrated estimate with uncertainty
- Use learned non-linearity to improve all 10K readings
Results:
- Naive (linear regression): MSE 5.2
- GP (RBF kernel): MSE 2.1 (60% improvement)
- Uncertainty quantification: Critical for downstream use
## 11. Integration with Other Methods
### 11.1 Kernel Methods with Deep Learning
Use GP as top layer (Bayesian uncertainty on deep features):
$$\mathbf{z} = ext{CNN}(x)$$
$$y \sim \mathcal{GP}(\mathbf{z}, k)$$
Benefits: Deep learning features + GP uncertainty.
### 11.2 Multiple Kernel Learning
Learn weighted combination of kernels:
$$k = \sum_i w_i k_i$$
Different kernels for different feature groups. Improves interpretability.
### 11.3 Kernel + Ensemble
Combine GP with other models:
$$y = 0.7 \cdot ext{GP}(x) + 0.3 \cdot ext{RF}(x)$$
Leverages strengths of each.
### 11.4 Kernelized Feature Selection
Use kernel to identify important features:
$$ ext{Importance}_i = \frac{\partial}{\partial \ell_i} \log p(y | X)$$
Gradient of marginal likelihood w.r.t. length scales reveals which features matter.
## 12. Future Research Directions
### 12.1 Scalable Gaussian Processes
Making GPs practical for >1M samples:
- Variational inference
- Inducing points
- Random features approximation
### 12.2 Deep Kernel Learning
Combine deep learning with kernel methods:
- Learn kernel from data
- End-to-end training
- Benefits of both
### 12.3 Non-Euclidean Kernels
Kernels on graphs, manifolds, structured data:
- Graph neural networks
- Kernel on sets/sequences
- Richer structure modeling
### 12.4 Causal Inference with GPs
Learn causal structure:
- Causal kernels
- Interventional data
- Counterfactual prediction
## 13. Summary & Key Takeaways
Kernel Methods:
- SVM: Fast training/prediction, good classification
- Gaussian Processes: Principled uncertainty, sample-efficient
- Kernel Ridge Regression: Simple, interpretable
Kernels:
- RBF: Default, smooth, flexible
- Matérn: Control smoothness
- Linear: High-dimensional sparse
- Periodic: Seasonal patterns
Typical Performance:
- Small datasets (n < 1000): GP or SVM excellent
- Medium datasets (1000 < n < 10K): Sparse GP competitive
- Large datasets (n > 100K): Neural networks usually better
Uncertainty Quantification:
- GPs: Well-calibrated, 95% confidence intervals accurate
- SVM: Confidence via probability calibration or distance-based
- Neural networks: Ensembles or Bayesian NN less calibrated
Hyperparameters:
- RBF kernel: Length scale key (tune via marginal likelihood)
- Noise: Domain-dependent (0.01-0.1 of output variance typical)
- SVM C: 0.1-10 (tune via cross-validation)
- Sparse GP: Inducing points min(1000, n/10)
Computational Cost:
- Full GP:
$$ O(n^3) $$
- practical for n < 10K
- Sparse GP:
$$ O(m^3) $$
- practical for n < 1M (with m ~500-1000)
- SVM:
$$ O(n^2) $$
to
$$ O(n^3) $$
- practical for n < 1M
When to Use:
- Limited data + need uncertainty → GP
- Classification, need interpretability → SVM
- Large dataset, need scalability → Neural networks
- Optimizing expensive function → Bayesian optimization (GP-based)
Kernel methods and GPs remain essential for uncertainty quantification and sample efficiency. Complementary to deep learning; often combined in practice.
---
## Appendix: Practical Implementation Labs
### Lab 1: Gaussian Process Regression
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel as C
def train_gp_regressor(X_train, y_train):
"""Train GP with RBF kernel"""
# Kernel: constant * RBF + white noise
kernel = C(1.0, (1e-3, 1e3)) * RBF(1.0, (1e-2, 1e2)) + WhiteKernel(1e-5)
gp = GaussianProcessRegressor(kernel=kernel, alpha=1e-6,
normalize_y=True, n_restarts_optimizer=5)
gp.fit(X_train, y_train)
return gp
def predict_with_uncertainty(gp, X_test):
"""Predict and get confidence intervals"""
mean, std = gp.predict(X_test, return_std=True)
# 95% confidence interval
ci_lower = mean - 1.96 * std
ci_upper = mean + 1.96 * std
return mean, std, (ci_lower, ci_upper)### Lab 2: Support Vector Machine
from sklearn.svm import SVC, SVR
from sklearn.preprocessing import StandardScaler
def train_svm_classifier(X_train, y_train, C=1.0, kernel='rbf'):
"""Train SVM for classification"""
# Normalize features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_train)
svm = SVC(C=C, kernel=kernel, gamma='scale', probability=True)
svm.fit(X_scaled, y_train)
return svm, scaler
def svm_predict(svm, scaler, X_test):
"""Predict with SVM"""
X_scaled = scaler.transform(X_test)
predictions = svm.predict(X_scaled)
probabilities = svm.predict_proba(X_scaled)
return predictions, probabilities### Lab 3: Kernel Matrix and Kernel Trick
import numpy as np
def rbf_kernel(X1, X2=None, gamma=0.1):
"""Compute RBF kernel matrix"""
if X2 is None:
X2 = X1
# Pairwise squared distances
distances = np.sum(X1**2, axis=1, keepdims=True) - 2*[email protected] + np.sum(X2**2, axis=1)
# RBF kernel
K = np.exp(-gamma * distances)
return K
def polynomial_kernel(X1, X2=None, d=2, c=1):
"""Compute polynomial kernel"""
if X2 is None:
X2 = X1
K = (X1 @ X2.T + c) ** d
return K
def kernel_ridge_regression(X_train, y_train, K_train, lambda_reg=0.01):
"""Kernel ridge regression"""
n = K_train.shape[0]
alpha = np.linalg.solve(K_train + lambda_reg * np.eye(n), y_train)
return alpha
def kernel_predict(alpha, K_test):
"""Predict with kernel method"""
return K_test @ alpha### Lab 4: Bayesian Optimization with GP
from scipy.optimize import minimize
import numpy as np
def expected_improvement(x, gp, X_train, y_best, xi=0.01):
"""Acquisition function: Expected Improvement"""
mu, sigma = gp.predict(x.reshape(1, -1), return_std=True)
if sigma == 0:
return 0
Z = (mu - y_best - xi) / sigma
ei = (mu - y_best - xi) * scipy.stats.norm.cdf(Z) + sigma * scipy.stats.norm.pdf(Z)
return -ei # Negative because minimize
def bayesian_optimization(objective_func, bounds, n_iterations=20):
"""Optimize using Bayesian optimization"""
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
X_train = []
y_train = []
gp = GaussianProcessRegressor(kernel=RBF(), alpha=1e-6)
# Initial random points
for _ in range(3):
x = np.random.uniform(bounds[:, 0], bounds[:, 1])
y = objective_func(x)
X_train.append(x)
y_train.append(y)
X_train = np.array(X_train)
y_train = np.array(y_train)
for i in range(n_iterations):
# Update GP
gp.fit(X_train, y_train)
# Find next point to evaluate
res = minimize(lambda x: expected_improvement(x, gp, X_train, y_train.min()),
x0=X_train[y_train.argmin()],
bounds=bounds)
x_next = res.x
y_next = objective_func(x_next)
X_train = np.vstack([X_train, x_next])
y_train = np.append(y_train, y_next)
return X_train, y_train, gp