Exact Gaussian Processes in High Dimensions: Escaping the Curse of Dimensionality Through Structured Sparsity

Gaussian process regression stands as one of the most theoretically principled approaches in machine learning, offering not merely predictions but full probabilistic distributions over functions. Yet this elegance comes at a steep computational price. Conventional implementations demand O(N³) operations and O(N²) storage for N training points, rendering them impractical for large datasets. When confronted with high-dimensional input spaces, the situation deteriorates further; standard approaches to scalability often rely on inducing point approximations that sacrifice exactness for tractability.

In their recent work, "Don't Get Your Kroneckers in a Twist: Gaussian Processes on High-Dimensional Incomplete Grids," Højlund, Lykke-Møller, Moss, and Christiansen introduce CUTS-GPR, a method that challenges the assumption that exact Gaussian process inference must remain confined to modestly sized problems. By exploiting the interplay between additive kernel structures and carefully constructed incomplete grids, the authors demonstrate exact inference with billions of training points in thousands of dimensions, completing full hyperparameter optimization in hours for problems that would otherwise be computationally inaccessible.

The Geometry of Incomplete Grids

The foundation of CUTS-GPR rests upon a specific geometric structure for the training data. Traditional grid-based Gaussian processes utilize complete Cartesian product grids, where the input space X equals X^(1) × ... × X^(D), yielding N = n₁n₂...n_D points. This structure permits fast inference through Kronecker product decompositions, but N grows exponentially with dimensionality D, effectively limiting such methods to low-dimensional settings.

Højlund et al. instead employ incomplete grids defined by a cut level α. Rather than sampling the full factorial combination of all dimensions, the method includes only low-dimensional subgrids or "cuts": a reference point (0D), one-dimensional displacements along each axis (1D cuts), two-dimensional combinations (2D cuts), up to order α. Formally, the grid consists of mode combinations (MCs) where each subgrid is indexed by the specific dimensions displaced from the reference. The total number of points scales as O(D^α) rather than O(n^D), a polynomial reduction that fundamentally alters the tractability landscape for high-dimensional problems.

This construction resonates deeply with the physics of molecular systems, where potential energy surfaces (the target application) often exhibit dominant low-order interactions. The energy landscape of a molecule with D vibrational modes may be well-approximated by examining individual modes and pairwise couplings (α=2), rendering the exponential growth of full factorial sampling unnecessary. The incomplete grid thus represents not merely a computational convenience, but a physically informed sparsity pattern that aligns with the intrinsic structure of many scientific problems.

Algorithmic Architecture and the Fast Matrix-Vector Product

The computational breakthrough of CUTS-GPR emerges from the recognition that incomplete grids combined with additive kernels yield kernel matrices possessing exploitable algebraic structure. When the kernel decomposes as a sum of lower-dimensional contributions, the full kernel matrix K becomes a sum of Kronecker products restricted to the incomplete grid indices.

Rather than computing the Cholesky decomposition of the covariance matrix C = K + σ²I, which requires O(N³) operations, CUTS-GPR employs iterative methods requiring only matrix-vector