Computational Aspects of Kernels

Kernel methods represent one of the most elegant frameworks in machine learning—they enable us to work implicitly in high-dimensional (even infinite-dimensional) feature spaces while performing computations only in the original input space. The kernel trick seems almost magical: we can fit complex, nonlinear models using only pairwise similarity computations.

But there's a catch. A significant one.

The very mechanism that makes kernel methods so powerful—computing and storing pairwise kernel evaluations—becomes their Achilles' heel at scale. When you have a million data points, storing the kernel matrix alone requires petabytes of memory. When you have even ten thousand points, the $O(n^3)$ cost of matrix inversion makes training prohibitively slow.

This computational reality has profound implications for when and how we can deploy kernel methods in practice. Understanding these limitations is not merely academic—it's essential for any practitioner who wants to use kernel methods beyond toy datasets.

At the heart of every kernel method lies the kernel matrix (also called the Gram matrix), denoted $\mathbf{K}$. For a dataset of $n$ training examples ${\mathbf{x}_1, \ldots, \mathbf{x}_n}$ and a kernel function $k(\cdot, \cdot)$, the kernel matrix is defined as:

$$\mathbf{K} = \begin{bmatrix} k(\mathbf{x}_1, \mathbf{x}_1) & k(\mathbf{x}_1, \mathbf{x}_2) & \cdots & k(\mathbf{x}_1, \mathbf{x}_n) \ k(\mathbf{x}_2, \mathbf{x}_1) & k(\mathbf{x}_2, \mathbf{x}_2) & \cdots & k(\mathbf{x}_2, \mathbf{x}_n) \ \vdots & \vdots & \ddots & \vdots \ k(\mathbf{x}_n, \mathbf{x}_1) & k(\mathbf{x}_n, \mathbf{x}_2) & \cdots & k(\mathbf{x}_n, \mathbf{x}_n) \end{bmatrix}$$

This matrix is symmetric ($K_{ij} = K_{ji}$) because valid kernels satisfy $k(\mathbf{x}, \mathbf{x}') = k(\mathbf{x}', \mathbf{x})$. For proper kernels, it is also positive semi-definite (all eigenvalues are non-negative).

Memory requirements:

Each element of $\mathbf{K}$ is a floating-point number. Using 64-bit (8 bytes) double-precision floats—the standard for numerical computation—the memory required to store the full kernel matrix is:

$$\text{Memory} = n^2 \times 8 \text{ bytes}$$

Let's compute concrete numbers for realistic dataset sizes:

Why the full matrix?

One might ask: do we really need the entire kernel matrix? For many kernel methods, unfortunately, yes:

• Kernel Ridge Regression requires solving $(\mathbf{K} + \lambda \mathbf{I})\boldsymbol{\alpha} = \mathbf{y}$, involving the full matrix.
• Gaussian Processes compute the posterior mean and variance using $\mathbf{K}^{-1}$.
• Support Vector Machines need $\mathbf{K}$ entries for the quadratic programming solver.

The kernel matrix is not just stored—it's manipulated, decomposed, and inverted. This compounds the computational difficulty significantly.

Before we can even store or manipulate the kernel matrix, we must compute it. Each entry $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$ requires evaluating the kernel function on a pair of input vectors. The computational cost depends on both the kernel function and the input dimensionality.

Common kernel evaluation costs:

For input vectors $\mathbf{x}, \mathbf{x}' \in \mathbb{R}^d$:

Total cost to compute the kernel matrix:

For $n$ training points and per-evaluation cost $O(d)$ (considering the common RBF kernel):

$$\text{Kernel matrix computation} = O(n^2 \cdot d)$$

This seems manageable—it's just $n^2$ inner products. But for structured kernels (strings, graphs, trees), per-evaluation costs can dominate. A string kernel with $O(|s|^2)$ cost per pair on 10,000 sequences of length 1,000 requires:

$$10,000^2 \times 1,000^2 = 10^{14} \text{ operations}$$

Even at 10 billion operations per second, this takes over an hour just to compute the matrix, before any learning occurs.

Prediction-time costs:

After training, making predictions on new points also incurs kernel evaluation costs. For kernel ridge regression, the prediction for a new point $\mathbf{x}_*$ is:

$$\hat{y}* = \mathbf{k}*^T \boldsymbol{\alpha}$$

where $\mathbf{k}* = [k(\mathbf{x}, \mathbf{x}1), \ldots, k(\mathbf{x}, \mathbf{x}_n)]^T$ is the kernel vector between the test point and all training points.

Cost per prediction: $O(n \cdot d)$ for kernel evaluations plus $O(n)$ for the dot product with $\boldsymbol{\alpha}$.

This means prediction time scales linearly with training set size. A model trained on 100,000 points requires 100,000 kernel evaluations per prediction—orders of magnitude slower than parametric models where prediction depends only on model parameters, not dataset size.

Computing and storing the kernel matrix is expensive, but it's the matrix operations that truly dominate the computational cost of kernel methods. Let's analyze the key operations.

Problem 1: Solving linear systems

Kernel Ridge Regression requires solving:

$$(\mathbf{K} + \lambda \mathbf{I}) \boldsymbol{\alpha} = \mathbf{y}$$

The naive approach—computing the inverse and multiplying—costs $O(n^3)$ for standard matrix inversion. Even with optimized algorithms (LU decomposition, Cholesky for positive definite matrices), the complexity remains $O(n^3)$.

Why $O(n^3)$ is problematic:

The cubic complexity creates severe scaling limitations. Let's compute concrete execution times assuming a modern workstation performing $10^{10}$ floating-point operations per second (10 GFLOPS effective throughput):

Problem 2: Gaussian Process predictive variance

Gaussian Processes provide uncertainty quantification through the predictive variance. For a test point $\mathbf{x}_*$:

$$\text{Var}[f_] = k(\mathbf{x}_, \mathbf{x}*) - \mathbf{k}^T (\mathbf{K} + \sigma^2 \mathbf{I})^{-1} \mathbf{k}_$$

This requires computing $\mathbf{K}^{-1}$ (or more practically, solving $\mathbf{K}^{-1} \mathbf{k}_*$ for each test point). With $m$ test points, the naive approach costs $O(n^3 + m \cdot n^2)$.

Problem 3: Log marginal likelihood for hyperparameter optimization

Hyperparameter optimization in GPs requires evaluating:

$$\log p(\mathbf{y}|\mathbf{X}, \boldsymbol{\theta}) = -\frac{1}{2}\mathbf{y}^T \mathbf{K}{\boldsymbol{\theta}}^{-1} \mathbf{y} - \frac{1}{2}\log|\mathbf{K}{\boldsymbol{\theta}}| - \frac{n}{2}\log(2\pi)$$

Both the $\mathbf{K}^{-1}\mathbf{y}$ term and the log-determinant $\log|\mathbf{K}|$ require $O(n^3)$ computation. During gradient-based optimization with many iterations, this cost is incurred repeatedly.

Let's synthesize all the computational costs into a complete picture. For training a kernel method on $n$ samples with $d$-dimensional inputs and making predictions on $m$ test points:

Training Phase:

Prediction Phase (for $m$ test points):

Hyperparameter Optimization Multiplier:

In practice, we rarely train a single model. Hyperparameter optimization (choosing $\lambda$, kernel parameters, etc.) requires training many models:

• Grid search with $k$ hyperparameter configurations: Multiply training cost by $k$.
• Cross-validation with $v$ folds: Multiply by $v$ (though each fold has $n(1-1/v)$ points).
• Gradient-based optimization for GP hyperparameters: Each gradient step requires computing the log marginal likelihood and its gradients—both $O(n^3)$ operations.

A 5-fold cross-validated grid search over 20 configurations requires 100 model fits. At 50,000 points with 1.2 hours per fit, that's 120 hours (5 days) of computation.

Beyond raw computational cost, kernel matrices present significant numerical stability challenges. Understanding these issues is crucial for implementing reliable kernel methods.

The Condition Number Problem:

The condition number $\kappa(\mathbf{K}) = \lambda_{\max} / \lambda_{\min}$ measures how errors in the input are amplified during matrix operations. For kernel matrices:

• High condition numbers are common, especially with RBF kernels on clustered data
• As $n$ increases, eigenvalue spectra often become increasingly ill-conditioned
• Regularization helps by adding $\lambda$ to all eigenvalues: $\kappa(\mathbf{K} + \lambda\mathbf{I}) = (\lambda_{\max} + \lambda)/(\lambda_{\min} + \lambda)$

Cholesky Failure:

Cholesky decomposition requires the matrix to be positive definite (all eigenvalues strictly positive). In practice, numerical errors can cause:

1. Near-zero eigenvalues becoming slightly negative due to floating-point errors
2. Cholesky failure when the algorithm encounters a non-positive pivot

This is especially common when:

• The regularization parameter is very small
• Data points are nearly duplicate (kernel matrix becomes rank-deficient)
• Single-precision (32-bit) floats are used instead of double-precision (64-bit)

Mitigation strategies:

Let's walk through a concrete scenario to illustrate when exact kernel methods hit their limits.

Scenario: Time-series forecasting with Gaussian Processes

You're building a GP model for electricity demand forecasting. The dataset contains:

• 3 years of hourly measurements: $3 \times 365 \times 24 = 26,280$ data points
• 50 input features: time indicators, weather, lagged demand, etc.
• You want to predict demand 24 hours ahead with uncertainty quantification

Computational requirements (exact GP):

With 26,000 points, exact methods are at the edge of feasibility—possible on a workstation, but requiring careful memory management and overnight hyperparameter tuning.

Now scale up:

Suppose you want 10 years of data (87,600 points) or add minute-level resolution for the most recent year (adding 525,600 more points).

Real-world implications:

This computational cliff explains why kernel methods, despite their theoretical elegance, were largely abandoned in the deep learning era when datasets grew to millions or billions of examples. Neural networks, with their minibatch-based training, scale to virtually unlimited data sizes.

However, the approximation methods we'll study in subsequent pages—Random Fourier Features, Nyström approximation, and Sparse GPs—reclaim kernel methods' utility by reducing complexity to:

• $O(nm^2)$ time and $O(nm)$ space for rank-$m$ approximations
• This makes million-point datasets feasible with $m \approx 1,000-10,000$

These approximations are not second-class solutions; they often achieve near-exact accuracy while enabling kernel methods at previously impossible scales.

We have conducted a thorough analysis of the computational complexity of kernel methods, establishing the fundamental constraints that motivate all scalability techniques.

The core challenges:

The path forward:

The remaining pages in this module present the solutions to these challenges:

• Random Fourier Features (next page): Transform kernel methods into linear models via random sampling in Fourier space
• Nyström Approximation: Low-rank matrix approximation using a subset of training points
• Sparse Gaussian Processes: Introduce inducing points to reduce the effective dataset size
• Scalability Strategies: Combine approximations with distributed computing for massive datasets

These techniques don't just make kernel methods faster—they fundamentally change what's possible, bringing the elegance of kernel learning to modern big data applications.