Turbocharging Gaussian Process Inference with Approximate Sketch-and-Project

Gaussian processes (GPs) are fundamental in probabilistic modeling, offering flexible, non-parametric predictions with uncertainty quantification. Since Rasmussen and Williams (2006), efforts have focused on scaling GPs via sparse approximations, variational methods, and iterative solvers. While these methods have extended applicability, they still face computational bottlenecks at large scales, especially when data exceeds millions of samples. Traditional direct methods like Cholesky decomposition are limited to small datasets due to O(n^3) complexity. Iterative methods like PCG improve scalability but suffer from ill-conditioning issues, slowing convergence. SDD offers better scalability but lacks strong theoretical guarantees, especially under ill-conditioned matrices. Consequently, large-scale GP inference remains an open challenge, motivating the development of more robust, scalable algorithms.

The core challenge in large-scale GP inference is solving the linear system (K + λI)w = y efficiently when n is huge. As n grows, the kernel matrix K becomes large and often ill-conditioned, making iterative solvers slow or unstable. Existing methods either require extensive hyperparameter tuning (SDD) or fail to guarantee convergence under poor conditioning (PCG). Moreover, the computational cost of exact matrix inversion or factorization becomes prohibitive at scale. Addressing these issues requires algorithms that are both robust to ill-conditioning and capable of handling massive datasets with minimal tuning, enabling practical Bayesian inference in real-world applications.

The paper's key innovations include: 1) Integrating sketch-and-project with Nyström low-rank kernel approximation to reduce computational complexity; 2) Theoretical analysis based on determinantal point processes that guarantees rapid convergence along principal spectral components, independent of condition number for moderate accuracy; 3) Distributed GPU-based matrix operations to handle datasets with hundreds of millions of samples; 4) Incorporation of Nesterov acceleration to speed up convergence. These innovations collectively enable scalable, robust GP inference, bridging the gap between theoretical guarantees and practical performance in massive data regimes.

• �� Kernel low-rank approximation: Use Nyström method with random sampling to approximate K, reducing the problem size.
• �� Sketch-and-project framework: Randomly select row blocks, perform approximate projections onto subspaces, and iteratively refine solutions.
• �� Approximate preconditioning: Replace exact inverse with low-rank Nyström approximation, enabling faster matrix-vector products.
• �� Distributed computation: Partition data and matrices across multiple GPUs, parallelizing matrix multiplications to handle large datasets.
• �� Acceleration: Apply Nesterov momentum to improve convergence rate, especially in early iterations.
• �� Theoretical analysis: Use determinantal point process properties to establish convergence bounds independent of condition number for moderate accuracy targets.
• �� Hyperparameter tuning: Default settings for block size, damping, and acceleration parameters to ensure out-of-the-box performance.

Experiments involve benchmarking ADASAP against PCG and SDD on large datasets, including UCI regression benchmarks and a 330 million sample transportation dataset. Metrics include RMSE and negative log-likelihood, with hyperparameter tuning performed for fair comparison. The datasets span from hundreds of thousands to hundreds of millions of samples, testing scalability and robustness. GPU clusters with multiple units are used to demonstrate distributed efficiency. Ablation studies analyze the impact of Nyström approximation quality, block size, and acceleration parameters. Results consistently show ADASAP achieves faster convergence, lower errors, and better scalability, validating the theoretical advantages.

ADASAP consistently outperforms PCG and SDD across all benchmarks, reducing RMSE and NLL by significant margins. On the houseelec dataset (1.84 million samples), it converges in less than half the time of competitors, achieving the lowest error metrics. In the transportation dataset (>330 million samples), ADASAP successfully completes inference, demonstrating scalability. Theoretical analysis confirms that the convergence rate along dominant spectral directions is superlinear and condition number-free for moderate accuracy, explaining the empirical efficiency. The combination of low-rank approximation, distributed computing, and acceleration results in a robust, scalable solution for large-scale GP inference.

This method is suitable for large-scale Bayesian optimization, scientific modeling, and real-time data analysis where datasets are massive and computational resources are distributed. It enables practitioners to perform accurate probabilistic inference without prohibitive costs, facilitating advances in fields like genomics, materials science, and industrial process optimization. The approach also supports integration with deep feature extractors, broadening its applicability in complex, high-dimensional settings.

While highly scalable, the method's accuracy depends on the quality of Nyström approximation, which may degrade for certain kernels or data distributions. The distributed implementation requires substantial hardware infrastructure, and communication overhead can limit efficiency at extreme scales. In ultra-high precision tasks, the convergence slows down as the residual error approaches machine precision. Future improvements include adaptive sampling, better approximation schemes, and reduced communication overhead to further enhance robustness and efficiency.