Structured Kernel Interpolation (SKI)

• Structured Kernel Interpolation (SKI) is a scalable GP inference method that approximates dense kernel matrices by interpolating from inducing points on grids.
• It exploits structured covariance properties (e.g., Toeplitz, Kronecker) to achieve near-linear time and memory complexity for massive datasets.
• Extensions like Product SKI, sparse grids, and SoftKI address high-dimensional challenges, broadening SKI’s applicability in real-world large-scale problems.

Structured Kernel Interpolation (SKI) is a scalable methodology for approximate Gaussian process (GP) inference, enabling near-linear time and memory complexity by leveraging structured interpolation of covariance kernels at inducing points placed on grids. SKI generalizes and unifies classical inducing-point GP approximations, producing highly efficient, accurate kernel matrices for use with both stationary and nonstationary kernels in high-data regimes. The core principle is to approximate the full $n \times n$ kernel Gram matrix by interpolating between values at a small set of $m \ll n$ inducing points while exploiting structure (Toeplitz, Kronecker, lattice) for rapid computation (Wilson et al., 2015).

1. Foundations of Structured Kernel Interpolation

SKI is founded on the observation that standard inducing-point GP approximations—such as Subset of Regressors (SoR) or FITC—can be viewed as interpolation from a set of auxiliary locations $U = \{u_j\}_{j=1}^m$, but with typically dense $n \times m$ cross-covariances. SKI introduces a sparse interpolation weight matrix $W \in \mathbb{R}^{n \times m}$, constructed by local schemes (linear, cubic, inverse-distance), to interpolate the cross-covariance: $K_{XU} \approx W K_{UU}$ This leads to the low-rank, structured approximate kernel: $K_{XX} \approx W K_{UU} W^\top$ where $K_{UU}$ is the $m \times m$ covariance matrix among the inducing points (Wilson et al., 2015, Moreno et al., 1 Feb 2025).

The interpolation weights are determined independently of kernel hyperparameters and are extremely sparse: with cubic interpolation in $d$ dimensions, each row of $W$ has $4^d$ nonzeros at the nearest grid neighbors.

Key trade-offs arise in the number and layout of inducing points, interpolation order, and the regularity of the underlying kernel. For sufficiently smooth kernels (e.g., RBF), local cubic interpolation achieves fourth-order accuracy; lower-order interpolation is computationally cheaper but less accurate (Wilson et al., 2015).

2. Computational Complexity, Structure, and Linear-Time Inference

SKI achieves scalability through two core mechanisms: sparsity in the interpolation matrix and algebraic structure in the grid kernel matrix. Placing inducing points on Cartesian grids enables $K_{UU}$ to inherit Kronecker or Toeplitz structure. This structuring yields rapid matrix-vector multiplies (MVMs) for iterative solvers:

• Individual steps: Applying $W^\top$ or $W$ to a vector is $O(c n)$ (where $c$ is number of nonzeros per row).
• Applying $K_{UU}$: $O(m \log m)$ with Toeplitz structure via FFT in 1D, or $O(P m^{1+1/P})$ with Kronecker algebra in $P$ dimensions.
• Total per-iteration complexity: $O(n + m\log m)$, and memory $O(n + m)$ (Wilson et al., 2015, Yadav et al., 2021).

For massive datasets, per-iteration cost can be reduced to $O(m \log m)$ after an $O(n)$ preprocessing pass computing $W^\top W$ and $W^\top y$, reframing inference as Bayesian linear regression in the $m$-dimensional basis of interpolation functions (Yadav et al., 2021). This yields dramatic memory and runtime reductions, enabling GP inference with $n \gtrsim 10^8$.

SKI also supports constant-time online updates: the sufficient statistics required for posterior and marginal likelihood computation are of size $m \times m$ and can be updated in $O(m^2)$ per new observation, independent of the growing data size (Stanton et al., 2021).

3. Interpolation Strategies and Extensions for High Dimensions

Standard SKI suffers from exponential scaling in $d$ due to $4^d$ nonzeros per row for cubic interpolation and the $m = k^d$ growth of full grids ($k$ per axis). Several extensions address this limitation:

• Product-Kernel Interpolation ("SKIP"): For kernels factorizing across input dimensions, SKI is applied to each 1D factor. The cross-dimensional interpolation and Kronecker product deliver $O(d n + d m \log m)$ complexity, reducing the exponential dependence on $d$ to linear (Gardner et al., 2018).
• Sparse Grid SKI: Sparse grids replace dense Cartesian grids, providing $m = O(2^\ell \ell^{d-1})$ points for grid level $\ell$, with simplex-based interpolation achieving $O(d \binom{\ell+d-1}{d-1})$ sparsity per row in $W$. Nearly linear-time MVMs in $O(m\,\mathrm{polylog}\,m)$ are attainable, and theoretical error bounds guarantee geometric convergence under kernel smoothness assumptions (Yadav et al., 2023).
• Simplex-GP/Permutohedral Lattice SKI: The permutohedral lattice dramatically reduces neighbor count per data point to $(d+1)$, replacing Cartesian grids by a simplicial tiling, and admits $O(d^2 (n + m))$ MVMs. This enables SKI to scale exponentially faster with $d$ and leverage highly parallel GPU acceleration (Kapoor et al., 2021).
• Soft Kernel Interpolation (SoftKI): SKI with softmax-based, learned inducing points generalizes interpolation to non-grid settings, with inference costs $O(m^2 n)$ and no exponential dependence on $d$. SoftKI remains efficient for $d$ up to hundreds, outperforming conventional Nyström approaches in kernel fidelity for moderate $m$ (Camaño et al., 2024).

4. Accuracy, Error Analysis, and Theoretical Guarantees

Convolutional cubic interpolation with a regular grid achieves a pointwise kernel approximation error decaying as $O(m^{-3/d})$ in $d$ dimensions (Moreno et al., 1 Feb 2025). The spectral norm error of the SKI Gram matrix scales as $O(n\,m^{-3/d} c^{2d})$ ($c$ is the interpolation weight norm), with rigorous bounds established for hyperparameter estimation and posterior mean/covariance prediction.

There are two scaling regimes:

• Low-dimensional regime ($d \leq 3$): Any fixed spectral error tolerance $\epsilon$ can be achieved in linear time by taking $m = O(n^{d/3})$.
• High-dimensional regime ($d > 3$): Keeping linear time costs requires relaxing the error target: $m$ must still grow as $O((n/\epsilon)^{d/3})$, but $m\log m = O(n)$ is violated unless $\epsilon$ increases with $n$.

Posterior mean and covariance errors propagate the SKI kernel approximation error linearly or polynomially in $n$, subject to the chosen $m$. For fixed $d$, increasing $m$ reduces error at the cost of increased computation (Moreno et al., 1 Feb 2025).

5. Methodological Extensions: Derivatives and Krylov Solvers

SKI adapts to GP regression problems requiring derivatives—such as in modeling fields as gradients of scalar potentials—by differentiating the interpolation operators. In D-SKI, partial derivative weights are constructed for each observation, and the resulting approximate covariance blocks are built via these differentiated interpolation matrices (Menzen et al., 2023). Inference is performed using preconditioned conjugate gradients (PCG) and the Lanczos tridiagonalization (LOVE) method for fast predictive variance evaluation.

This strategy achieves nearly linear scaling in $N$, subquadratic scaling in the number of inducing grid points, and enables large-scale modeling (e.g., 40,000 3D vector-valued observations in minutes on commodity hardware) (Menzen et al., 2023).

6. Applications and Empirical Performance

SKI and its extensions are effective in diverse settings:

• Kernel learning: When placed on large grids ($m \gg n$), SKI recovers ground-truth spectral or compositional structure more accurately and at lower computational cost than SoR/FITC, even with expressive nonseparable kernels (Wilson et al., 2015).
• Spatiotemporal and environmental datasets: SKI enables large-scale weather radar and magnetic field mapping, supporting tens of thousands to over $10^8$ data points with linear scaling (Yadav et al., 2021, Menzen et al., 2023).
• Audio and time-series modeling: SKI delivers orders-of-magnitude speedups and improved or matched accuracy over traditional MVM methods (Wilson et al., 2015, Yadav et al., 2021).
• High-dimensional regression: Simplex-GP, sparse-grid SKI, and SoftKI scale accurate GP inference to moderate and high $d$ ($d \lesssim 10$ for sparse-grid, $d \lesssim 1000$ for SoftKI) (Yadav et al., 2023, Camaño et al., 2024).

The table summarizes methods addressing the curse of dimensionality in SKI:

Extension	Dimensionality Range	Complexity per MVM
Product SKI	$d \lesssim 12$	$O(d n + d m \log m)$
Sparse-grid SKI	$d \lesssim 10$	$O(m \,\mathrm{polylog}\,m)$
Simplex/Permutohedral Lattice	$d \lesssim 20$	$O(d^{2}(n + m))$
SoftKI	$d \gg 10$	$O(m^2 n)$

7. Limitations and Open Challenges

Despite its scalability, classical SKI faces inherent limits in ambient dimension due to exponential scaling of standard grid-based interpolation. Overcoming these limits is the major focus of recent research via sparse grids, simplicial lattices, product kernel decompositions, and flexible learned interpolations.

SKI’s accuracy is contingent on kernel smoothness and grid density, with theoretical and empirical guidance on grid size selection (e.g., $m \sim n^{d/3}$ for cubic) (Moreno et al., 1 Feb 2025). For nonstationary, highly structured, or high-$d$ data, advanced SKI variants or integration with deep kernel learning may be necessary.

Recent advances provide constant-time online updates and very fast inference after initial preprocessing, but storage of large $m \times m$ structures (e.g., for covariance matrices) can remain a bottleneck in ultra-large $m$ regimes (Yadav et al., 2021, Stanton et al., 2021). Accelerated GPU implementations and Krylov subspace methods further push SKI’s applicability to real-time, streaming, and high-throughput domains.

• (Wilson et al., 2015) Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP)
• (Menzen et al., 2023) Large-scale magnetic field maps using structured kernel interpolation for Gaussian process regression
• (Moreno et al., 1 Feb 2025) The Price of Linear Time: Error Analysis of Structured Kernel Interpolation
• (Yadav et al., 2021) Faster Kernel Interpolation for Gaussian Processes
• (Stanton et al., 2021) Kernel Interpolation for Scalable Online Gaussian Processes
• (Yadav et al., 2023) Kernel Interpolation with Sparse Grids
• (Kapoor et al., 2021) SKIing on Simplices: Kernel Interpolation on the Permutohedral Lattice for Scalable Gaussian Processes
• (Camaño et al., 2024) High-Dimensional Gaussian Process Regression with Soft Kernel Interpolation
• (Gardner et al., 2018) Product Kernel Interpolation for Scalable Gaussian Processes