Kernel Ridge Regression

Kernel ridge regression is mathematically elegant and statistically principled. But does it scale? The answer is nuanced: kernel methods face fundamental computational challenges that limit their applicability to large datasets.

In this page, we dissect the computational complexity of kernel ridge regression. We'll analyze every step—from kernel matrix construction to prediction—understanding where time and space are consumed. This analysis isn't merely academic: it determines whether kernel methods are viable for your problem.

Understanding these limitations is crucial for:

• Deciding when kernel methods are appropriate
• Estimating computational resources needed
• Choosing between exact and approximate methods (covered in Module 5)
• Designing efficient implementations

Let's systematically analyze the computational cost of training kernel ridge regression. We have $n$ training points with $d$-dimensional inputs.

Step 1: Kernel Matrix Construction

We need to compute $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$ for all pairs $i, j$.

Number of kernel evaluations: $\frac{n(n+1)}{2}$ (exploiting symmetry $K_{ij} = K_{ji}$)

Cost per kernel evaluation: Depends on kernel:

• Linear kernel $k(\mathbf{x}, \mathbf{x}') = \mathbf{x}^\top \mathbf{x}'$: $O(d)$ operations
• Polynomial kernel $(\mathbf{x}^\top \mathbf{x}' + c)^p$: $O(d)$ operations (compute inner product, then power)
• RBF kernel $\exp(-\gamma|\mathbf{x} - \mathbf{x}'|^2)$: $O(d)$ operations

Total for kernel matrix: $O(n^2 d)$ time, $O(n^2)$ space

Step 2: Solving the Linear System

We need to solve $(\mathbf{K} + \lambda\mathbf{I})\boldsymbol{\alpha} = \mathbf{y}$ for $\boldsymbol{\alpha}$.

Direct Methods:

Cholesky Decomposition Process:

1. Compute $\mathbf{L}$ such that $\mathbf{K} + \lambda\mathbf{I} = \mathbf{L}\mathbf{L}^\top$: $O(n^3/3)$
2. Solve $\mathbf{L}\mathbf{z} = \mathbf{y}$ (forward substitution): $O(n^2)$
3. Solve $\mathbf{L}^\top \boldsymbol{\alpha} = \mathbf{z}$ (back substitution): $O(n^2)$

Total: $O(n^3)$ dominated by Cholesky factorization

Iterative Methods:

When $n$ is very large or $\mathbf{K}$ has favorable structure, iterative methods can be faster.

Convergence: CG converges in at most $n$ iterations to exact solution; in practice, $k << n$ iterations often suffice for good approximation, especially when $\lambda$ is not too small.

Break-even: Iterative beats direct when $k < n/3$ approximately.

Space for solving: $O(n)$ for CG (beyond storing $\mathbf{K}$); $O(n^2)$ for Cholesky (storing $\mathbf{L}$)

After training, we have the dual coefficients $\boldsymbol{\alpha}^* \in \mathbb{R}^n$. For a new test point $\mathbf{x}_*$:

Single Point Prediction:

1. Compute kernel vector $\mathbf{k}_* = (k(\mathbf{x}1, \mathbf{x}), \ldots, k(\mathbf{x}n, \mathbf{x}))^\top$

   • Cost: $O(nd)$ for $n$ kernel evaluations, each $O(d)$

2. Compute prediction $f(\mathbf{x}*) = \mathbf{k}^\top \boldsymbol{\alpha}^$

   • Cost: $O(n)$ for dot product

Total per prediction: $O(nd)$ time

Batch Prediction (m test points):

1. Compute test kernel matrix $\mathbf{K}_* \in \mathbb{R}^{n \times m}$: $O(nmd)$ time
2. Matrix-vector product $\mathbf{f}^* = \mathbf{K}_^\top \boldsymbol{\alpha}^$: $O(nm)$ time

Total: $O(nmd)$ time, $O(nm)$ space

Comparison to Other Models:

How does kernel ridge regression compare to simpler models?

Key Insight:

Linear regression has $O(d)$ prediction, constant in $n$—the training data is "compressed" into $d$ parameters. Kernel ridge regression has $O(nd)$ prediction, growing with $n$—we must "consult" all training points for each prediction.

This means:

• Adding more training data makes KRR predictions slower
• For real-time applications with large $n$, exact KRR may be too slow
• Approximations that reduce effective $n$ are crucial (sparse methods)

Memory is often the first bottleneck encountered in practice. Let's quantify requirements precisely.

Storage Components:

1. Training data: $n \times d$ matrix = $O(nd)$ values
2. Kernel matrix $\mathbf{K}$: $n \times n$ = $O(n^2)$ values
3. Cholesky factor $\mathbf{L}$ (if stored): $n \times n$ = $O(n^2)$ values (can overwrite K)
4. Dual coefficients $\boldsymbol{\alpha}$: $n$ values = $O(n)$
5. Targets $\mathbf{y}$: $n$ values = $O(n)$

Dominant Term: The kernel matrix at $O(n^2)$ dominates unless $d$ is huge.

Memory-Efficient Strategies:

1. On-the-fly kernel computation:

Don't store $\mathbf{K}$ explicitly. Recompute $K_{ij}$ when needed.

• Trades memory for time: no $O(n^2)$ storage, but each element computed multiple times
• Suitable for iterative solvers (CG) that access entries in structured patterns
• Cost: training becomes $O(n^3 d)$ instead of $O(n^2 d + n^3)$

2. Exploiting sparsity:

If the kernel "effectively" produces sparse matrices (e.g., compact support kernels), only store non-zeros.

• RBF with small $\sigma$: many entries near-zero but not exactly zero
• Truncated kernels (set $k = 0$ below threshold): explicit sparsity

3. Blocked/out-of-core computation:

Store $\mathbf{K}$ on disk, process in blocks that fit in RAM.

• Dramatically slower (I/O bound)
• Can handle very large $n$ with sufficient disk space

4. GPU memory:

Modern GPUs have 16-80 GB RAM with much higher bandwidth.

• n ≈ 45,000 fits entirely in 16 GB GPU memory
• GPU matrix operations are highly optimized

We've focused on the dual formulation, but the primal exists too. When is each preferable?

Primal Formulation:

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

Requirements:

• Explicit feature map $\phi(\mathbf{x}) \in \mathbb{R}^D$ must be computable
• Store and process $D \times D$ matrix (feature covariance)

Complexity:

• Build $\mathbf{\Phi}$: $O(nD)$ space, $O(nD \cdot t_\phi)$ time where $t_\phi$ is cost to compute $\phi$
• Build $\mathbf{\Phi}^\top \mathbf{\Phi}$: $O(nD^2)$ time, $O(D^2)$ space
• Solve: $O(D^3)$ time
• Store weights: $O(D)$ space
• Predict: $O(D)$ per point

Dual Formulation:

$$\boldsymbol{\alpha}^* = (\mathbf{K} + \lambda\mathbf{I}_n)^{-1}\mathbf{y}$$

Requirements:

• Kernel function $k(\mathbf{x}, \mathbf{x}')$ must be computable
• Store and process $n \times n$ matrix (kernel matrix)

Complexity:

• Build $\mathbf{K}$: $O(n^2 d)$ time, $O(n^2)$ space
• Solve: $O(n^3)$ time
• Store coefficients: $O(n)$ space
• Predict: $O(nd)$ per point

Detailed Comparison:

Practical Guidelines:

1. If using RBF, Matérn, or other infinite-dimensional kernels: must use dual
2. If using linear kernel: primal is often better (reduce to standard ridge regression)
3. If using polynomial kernel with low degree: consider primal with explicit features
4. If $n > 10D$: dual becomes prohibitively expensive relative to primal
5. If $D > 10n$: primal becomes prohibitively expensive relative to dual

Kernel ridge regression has hyperparameters:

• Regularization strength $\lambda$
• Kernel-specific parameters (e.g., RBF bandwidth $\sigma$, polynomial degree $p$)

Selecting these via cross-validation multiplies computational cost.

K-Fold Cross-Validation:

For each hyperparameter configuration:

1. Split data into K folds
2. For each fold: train on K-1 folds, evaluate on held-out fold
3. Average performance across folds

Naive cost: $K \times (\text{training cost})$ per configuration

With $g$ grid points per hyperparameter and $h$ hyperparameters: $g^h \times K \times O(n^3)$

Example: 10 values of $\lambda$, 10 values of $\sigma$, 5-fold CV on 10,000 points: $$10 \times 10 \times 5 \times O(10000^3) = 500 \times O(10^{12}) \approx 5 \times 10^{14} \text{ operations}$$

This is catastrophically expensive!

Efficient Hyperparameter Search Strategies:

1. Eigendecomposition Reuse:

Compute $\mathbf{K} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^\top$ once: $O(n^3)$

For any $\lambda$: $$(\mathbf{K} + \lambda\mathbf{I})^{-1} = \mathbf{U}(\mathbf{\Lambda} + \lambda\mathbf{I})^{-1}\mathbf{U}^\top$$

LOO-CV for many $\lambda$ values: $O(n^3)$ upfront + $O(n^2)$ per $\lambda$

2. Kernel Parameter Sensitivity:

Changing kernel parameters (e.g., $\sigma$ in RBF) requires rebuilding $\mathbf{K}$ entirely.

Strategy: Coarse-to-fine search

1. Rough grid of $(\sigma, \lambda)$: few models
2. Refinement around best: more models in promising region

3. Bayesian Optimization:

Use Gaussian Processes (ironically!) to model the validation error surface as a function of hyperparameters. Intelligently select next configuration to evaluate.

• Far fewer total evaluations than grid search
• Each evaluation still costs $O(n^3)$, but total count is much lower

Modern hardware offers massive parallelism. Which parts of kernel ridge regression can exploit it?

Highly Parallelizable Components:

1. Kernel Matrix Construction:

Each entry $K_{ij}$ can be computed independently.

• Embarrassingly parallel: distribute $(i, j)$ pairs across workers
• GPU-friendly: batch kernel evaluations as matrix operations
• Speedup: effectively linear in number of cores

GPU optimization: For RBF kernel, compute pairwise squared distances first (highly optimized BLAS operations), then apply exponential.

2. Matrix-Matrix Products:

When predicting for multiple test points: $\mathbf{f}^* = \mathbf{K}_^\top \boldsymbol{\alpha}^$

• Standard BLAS Level 3 operation
• Excellent CPU and GPU implementations exist

3. Cross-Validation Folds:

For K-fold CV, each fold is independent and can run on a separate machine/core.

GPU Acceleration:

GPUs excel at kernel ridge regression due to:

• Massive parallelism for kernel matrix construction
• Optimized BLAS for matrix operations
• High memory bandwidth for dense matrix access

Typical speedups:

• Kernel matrix construction: 50-100× over single CPU core
• Cholesky factorization: 5-20× over CPU
• Overall training: 10-30× depending on problem size

Multi-GPU and Distributed:

For very large problems:

• Distribute kernel matrix blocks across GPUs
• Parallel matrix factorization algorithms (ScaLAPACK, MAGMA)
• Communication becomes a bottleneck

Hardware Landscape:

Given the computational reality, when should you give up on exact kernel ridge regression?

Signs It's Time for Approximations:

1. Memory exceeded: If $n^2 \times 8$ bytes > available RAM, you cannot store the kernel matrix

2. Training time unacceptable: If $O(n^3)$ exceeds your time budget

3. Prediction latency too high: If $O(nd)$ per prediction is too slow for real-time requirements

4. Hyperparameter search prohibitive: If $O(n^3)$ per configuration times number of configurations exceeds budget

Approximate Methods Preview (Module 5):

Decision Framework:

IF n ≤ 5000:
    Use exact kernel ridge regression
ELSE IF n ≤ 20000 AND have GPU:
    Use exact with GPU acceleration
ELSE IF n ≤ 50000:
    Consider iterative solver (CG) if well-conditioned
ELSE:
    Use approximate methods (Nyström, Random Features)

The Accuracy-Efficiency Trade-off:

Approximate methods sacrifice some accuracy for computational tractability:

We've comprehensively analyzed the computational landscape of kernel ridge regression. Let's consolidate the key findings:

What's Next:

With computational reality understood, we move to Kernel Selection—the art and science of choosing appropriate kernels for your problem:

• Properties of different kernel families
• How to match kernels to data characteristics
• Hyperparameter tuning strategies
• Domain-specific kernel design
• When different kernels are appropriate