Gaussian processes constitute a cornerstone of probabilistic machine learning, yet scaling them to large datasets typically forces a trade-off between computational efficiency and model fidelity. This work bridges this gap by presenting an adaptive multi-resolution Gaussian process framework that is both scalable and exact. Our key innovation is constructing a naturally data-sparse covariance matrix with adaptive multi-resolution basis functions. These basis functions are directly anchored to samples, eliminating the need for auxiliary points. By shrinking the support domains of multi-resolution basis, the matrix block sizes are limited, guaranteeing sparsity. The inverse of the data-sparse covariance matrix is computed exactly and efficiently via the sparse Cholesky inverse algorithm. To further improve predictive uncertainties, we construct an augmented basis function. Theoretical analysis and numerical experiments demonstrate that our model achieves exact inference with O(nlog2n) training cost and O(logdn) prediction cost, establishing a principled framework for scalable and high-fidelity Gaussian process regression.
Figures & tables
Figure 1: Examples of our adaptive multi-resolution basis for samples with ρ=2.5 in [0,10] .
Figure 2: For this tree structure, the nodes are collected in nested dissection order as {4,5,2,8,9,6,7,3,1} .
Figure 3: Sparse pattern of matrices in an AMRGP model with 10,000 samples.
Table 1: Synthetic data used in numerical experiments.
Figure 4: Synthetic data.
Figure 5: OpenMP parallel computing with 1 million data points.
Figure 6: Performance of our AMRGP on Case 7 with different ρ ’s and m ’s.
Figure 7: Performance of our AMRGP on Case 8 with different ρ ’s and m ’s.
Case
1
2
3
4
5
6
Time (s)
25.9
20.5
19.9
23.4
20.6
20.2
Memory (GB)
11.5
11.5
11.5
11.5
11.5
11.5
Normalized CRPS
0.004
0.002
0.345
0.009
0.013
0.024
Table 2: Performance of the choice ρ=4.0 and m=100 on Case 1–6 with 1 million training data points.
Figure 8: Complexity of our AMRGP.
Case
Dimension
Metric
SVGP
SKIP
SVecGP
AMRGP
1
1
Time cost (s)
123.4 (2.5)
77.7 (1.5)
58.8 (7.1)
5.1 (3.8)
Normalized CRPS × 100
0.12
0.02
0.01
0.10
2
1
Time cost (s)
125.3 (2.5)
86.1 (1.5)
57.4 (6.9)
3.3 (1.6)
Normalized CRPS × 100
2.11
2.73
0.07
0.21
3
2
Time cost (s)
131.7 (2.6)
176.0 (3.4)
84.0 (8.9)
2.8 (2.3)
Normalized CRPS × 100
0.22
0.59
1.15
0.99
Table 3: Comparison of scalable GP models. The time in parenthesis is the average time cost by one iteration in approximate GPs and the training time in AMRGP.
Figure 9: Predictions for Case 2.
Figure 10: Predictions for Case 4: exact and predicted values (top) and standard deviations (bottom).
Figure 11: Prediction for the Burger’s equation over the computational domain.
Figure 12: Predicted for the wave equation over the computational domain. The black rectangle denotes the test region.
Gaussian processes (GPs) offer appealing properties but are costly to train at scale. Sparse variational GP (SVGP) approximations reduce cost yet still rely on Cholesky decompositions of kernel matrices, ill-suited to low-precision, massively parallel hardware. While one can construct valid variational bounds that rely only on matrix multiplications (matmuls) via an auxiliary matrix parameter, optimising them with off-the-shelf first-order methods is challenging. We make the inverse-free approach practical by proposing a better-conditioned bound and deriving a matmul-only natural-gradient update for the auxiliary parameter, markedly improving stability and convergence. We further provide simple heuristics, such as step-size schedules and stopping criteria, that make the overall optimisation routine fit seamlessly into existing workflows. Across regression and classification benchmarks, we demonstrate that our method 1) serves as a drop-in replacement in SVGP-based models (e.g., deep GPs), 2) recovers similar performance to traditional methods, and 3) can be faster than baselines when well tuned.
We present a theoretically grounded Gaussian process framework that leverages neural feature maps to construct expressive kernels. We show that the learned feature map can be interpreted as an optimal low-rank approximation to a Gram matrix derived from an implied RKHS, from which we establish consistency of the GP posterior. We further analyse the spectral properties of the induced kernels and introduce product feature-map kernels to address oversmoothing. This simple yet powerful approach enables fast, scalable, and accurate exact GP inference with minimal upfront work. The flexibility of kernel design supports seamless application to both regression and classification tasks across diverse data modalities, including tabular inputs and structured domains such as images. On benchmark datasets, this approach surpasses pre-existing methods in terms of accuracy and training and prediction efficiency.
Gaussian processes (GPs) provide principled probabilistic predictions while encoding prior knowledge, including equivariances. Yet, their use in large-scale scientific problems is limited by computational cost. Equivariant neural networks are common but typically lack the uncertainty quantification offered by GPs, which is valuable in applications such as molecular research. High-dimensional inputs and large symmetry groups further demand scalability. We establish results pertaining to the interplay of GP equivariance and conditioning and leverage them to obtain equivariant sparse GPs through suitable mean functions and covariance kernels. We instantiate this framework with a flexible class of integration-free equivariant kernels, yielding scalable and data-efficient GP inference. In particular, we introduce triply scalable equivariant Gaussian processes. We employ equivariant sparse variational Gaussian processes for SO(2)-equivariant vector fields and molecular property prediction. Alongside the SVGP, we develop a matrix-free equivariant full-GP implementation that combines an exact Kronecker reduction with preconditioned conjugate-gradient solves, enabling fast and scalable evaluation of the full joint predictive density. We further compare different approaches for selecting inducing points in the equivariant sparse GP models. Our test cases include synthetic SO(2)-equivariant fields as well as the prediction of electric dipole moments of N-methylformamide based on quantum chemistry simulations, achieving accurate, uncertainty-aware predictions at a fraction of the computational cost of classical GP inference.
Tim Steinert, David Ginsbourger
Institute of Mathematical Statistics and Actuarial Science, University of Bern, Bern, Switzerland