Turbocharging Gaussian Process Inference with Approximate Sketch-and-Project
Authors: Pratik Rathore
Zachary Frangella
Sachin Garg
Shaghayegh Fazliani
Michal Derezinski
Madeleine Udell
Venue: NeurIPS 2025 poster ()
Keywords: Gaussian processes
inference
sketch-and-project
distributed methods
determinantal point processes
Primary area: probabilistic_methods
TL;DR: We develop a distributed, approximate, accelerated sketch-and-project method that scales Gaussian process inference to hundreds of millions of samples.

Abstract: Gaussian processes (GPs) play an essential role in biostatistics, scientific machine learning, and Bayesian optimization for their ability to provide probabilistic predictions and model uncertainty.
However, GP inference struggles to scale to large datasets (which are common in modern applications), since it requires the solution of a linear system whose size scales quadratically with the number of samples in the dataset.
We propose an approximate, distributed, accelerated sketch-and-project algorithm ($\texttt{ADASAP}$) for solving these linear systems, which improves scalability.
We use the theory of determinantal point processes to show that the posterior mean induced by sketch-and-project rapidly converges to the true posterior mean. 
In particular, this yields the first efficient, condition number-free algorithm for estimating the posterior mean along the top spectral basis functions, showing that our approach is principled for GP inference.
$\texttt{ADASAP}$ outperforms state-of-the-art solvers based on conjugate gradient and coordinate descent across several benchmark datasets and a large-scale Bayesian optimization task.
Moreover, $\texttt{ADASAP}$ scales to a dataset with $> 3 \cdot 10^8$ samples, a feat which has not been accomplished in the literature.

Turbocharging Gaussian Process Inference with
Approximate Sketch-and-Project
Pratik Rathore Stanford University pratikr@stanford, edu
Zachary Frangella Stanford University 2fran@stanford. edu
Sachin Garg University of Michigan sachgOumice ch, edu
Shaghayegh Fazliani Stanford University fazliani wawaaniitttafrrro edu
Michal Dereziński University of Michigan derezin@umich edu
Madeleine Udell Stanford University udel1@stanf ord. edu
Abstract
Gaussian processes (GPs) play essential role in biostatisties, scientific machine learning, and Bayesian optimization for their ability to provide probabilistic pre- dictions and model uncertainty. However, GP inference struggles to scale large datasets (which are common in modern applications), since it requires the solution linear system whose size scales quadratically with the number of samples the dataset. We propose an approximate, distributed, accelerated sketch-and-proof algorithm (ADASAP) for solving these linear systems, which improves scalability. We use the theory of determinantal point processes to show that the posterior mean induced by sketch-and-projee rapidly converges to the true posterior mean. In par- ticular, this yields the first efficient, condition number-free algorithm for estimating the posterior mean along the top spectral basis functions, showing that our approach is principled for GP inference. ADASAP outperforms oftte-of-the-ar solvers based on conjugate gradient and coordinate descent across several benchmark datasets and a large-scale Bayesian optimization task. Morcover, ADASAP scales to dataset with > 3 10 samples, a feat which has not been accomplished the literature.
Introduction
Gaussian processes (GPs) are a mainstay of modern machine learning and scientific computing, due to their ability provide probabilistic predictions and handle uncertainty quantification. Indeed, GPs arise in applications spanning Bayesian optimization [Hernández-Lobee
et al., 2017], genetics [McDowell et al., 2018], health care analytics [Cheng et al., 2020]. materials science [Frazier and Wang, 2016], and partial differential equations [Chen et al., 2025].
GP inference for dataset with samples requires the solution of a dense n x linear system, which is challenging to solve at large scale. Direct methods like Cholesky decomposition require O(n²) computation, limiting them to n 10a. Consequently, two main approaches have arisen for large-scale inference: (i) exact inference based on iterative methods solving linear systems [Wang et al., 2019, Lin et al., 2023, 2024] and (ii) approximate inference based on inducing points and variational methods [Titsias, 2009, Hensman et al., 2013]. The oftte-of-thh-ar approaches for large- scale inference preconditioned conjugate gradient (PCG) [Wang et al., 2019] and stochastic dual descent (SDD) [Lin et al., 2024]. PCG offers strong convergence guarantees and good performance on ill-conditioned problems, but is slow when 10 In contrast, SDD scales >>>>>>>>>>>> 10 but lacks strong theoretical guarantees, can slow down in the face of ill-conditioning, and introduces stepsize parameter which be challenging tune. Inducing points and variational methods [Titsias, 2009, Hensman et al., 2013] scale to large datasets, but are often outperformed by exact
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

houseciec (12 189)
houseeloc (n = 1.84- 10")
2.0
20
1.5
15
Toa RMSE
20
0.0
200
400
600
800
1000
1200
200
400
000
800
1000
1200
(4)
Figure 1: ADASAP attains lower root mean square error (RMSE) and mean negative log likelihood (NLL) than start-of-the-art methods SDD [Lin al., 2024] and PCG on the houseelec dataset. SDD-1 and SDD-10 10 correspond to two particular stepsize selections for the SDD method. The solid lines indicate the mean performance of each method, while the shaded regions indicate the range between worst and best performance of each method over five random splits of the data.
inference [Wang et al., 2019, Lin et al., 2023, 2024]. Altogether, the current state of algorithms for large-scale GP inference is unsatisfactory, as practitioners have trade-off between quality of inference, robustness to ill-conditioning, ease of setting hyperparameters, and scalability.
To address this gap in the literature, introduce the Approximate Distributed Accelerated Sketch- and-Project (ADASAP) algorithm. ADASAP rooted in the saetch-and-projee framework Gower and Richtárik [2015], and obtains the robustness to ill-conditioning of PCG and the scalability of SDD, while having reliable, default hyperparameters that work well in practice. For example, ADASAP dramatically outperforms tuned SDD and PCG dataset with 10 samples (Fig. 1).
Our contributions are as follows:
1. We develop ADASAP for large-scale GP inference. ADASAP approximate preconditioning and acceleration to address ill-conditioning, and uses distributed computing to improve the speed of bottleneck operations. ADASAP also comes with effective default hyperparameters that work out of the box.
2. Using the theory of determinantal point processes, we show that sketch-and susich--nnd-rojjet--tttty methods converge faster than stochastic first-order methods like SDD in the presence of ill- conditioning. In particular, we give trst-off-tts-kin condition number-free free time complexity bound for estimating the GP posterior mean moderate accuracy, explaining the excellent test еттог performance of ADASAP.
3. We empirically verify that ADASAP, with its default hyperparameters, outperforms state-of- the-art competitors on benchmark large-scale scale GP inference tasks, and capable of scaling to a dataset with n > 3 10 samples.
4. We show ADASAP yields the best performance on a large-scale Bayesian optimization task from Lin al. [2023].
2 GP Regression and Inference
Let &' be a set. A random function f: X' R is called a Gaussian process if for any finite subset 22..2.......
C X of points (where is any positive integer), the random vector (f(x1), f(22).... f(r... follows multivariate Gaussian distribution [Rasmussen and Williams, 2005]. A Gaussian process is typically denoted as f ~ GP(m, (m,)), where m: - function and k: X' x X - IR is the covariance function (or kernel). Throughout this paper, the kemel function à follows broadcasting conventions, where operations between individual points and sets of points produce vectors matrices of kerel evaluations appropriate.
mean

---

Suppose we have Gaussian process prior f GP k). Given observations (X, E R"xd R"" and likelihood variance À > we would like to perform inference using the GP posterior, i.c., (i) sampling from the posterior and (ii) computing the posterior mean.
For conciseness, define K == xxx;
R Then sample from the GP posterior is random function fa GP(Mm. kn) [Rasmussen and Williams, 2005], where
())+) = m(-) ((; X) + A1) y/ () + \1)-)
(2)
The standard approach for drawing posterior samples at locations X, is perform linear transfor- mation of standard Gaussian random variables Ç _N(0, 1) [Wilson et al., 2020, 2021]:
fn(X.) = m. (X))= + kn (Xx,, X.)
Unfortunately, the approach in (3) does not scale when [X.] is large, since the matrix square root kn (X,, X.) requires O[ ()) computation.
2.1 Pathwise conditioning
To address the scaling challenges in (3), Wilson al. [2020, 2021] develop pathwise conditioning. Pathwise conditioning rewrites (3) as
= =
{~,N00 f g~ k). (4) Equation (4) enables scalable sampling from the posterior as f can be approximated using random features [Rahimi and Recht, 2007, Rudi and Rosasco, 2017] 1 for k, which costs O(q X- |). where q is the number of random features. This yields significant reduction in cost from the standard posterior sampling scheme in (3): by leveraging pathwise conditioning, posterior sampling can be performed by generating via random features [Wilson al., 2020, 2021] and solving (K + - (y f(X,,, - using an iterative method like PCG or SDD. Our experiments (Section 5) use pathwise conditioning to sample from the posterior.
2.2 Drawbacks of state-of-the-art exact inference methods
The current state-of-the-an for large-scale posterior sampling is SDD [Lin et al., 2024], a stochastic first-order algorithm combining coordinate descent [Richtárik and Takáč, 2014, Qu and Richtárik, 2016] with momentum and geometric averaging. Stochastic methods like SDD can better handle single (or lower) precision than PCG [Wu et al., 2024, Rathore et al., 2025], which is essential to maximizing GPU performance and reducing storage costs [Abdelfattah et al., 2021]. While SDD overcomes the scalability challenges associated with PCG, it faces challenges due to the ill- conditioning of + AI, which determines its worst-case convergence rate [Nemirovski and Yudin, 1983]. Therefore, reliable large-scale inference requires method that has cheaper per-iteration costs than PCG and can handle the ill-conditioning of K + XI. However, to the best of our knowledge, no existing method effectively tackles both issues: SGD and SDD offer speed and scalability but suffer from ill-conditioning, while PCG offers robustness to ill-conditioning but struggles to scale beyond 106. In Section 3 propose SAP, which is robust to ill-conditioning: in Section 4, develop ADASAP, an extension SAP that scales to large datasets.
Sketch-and-projet for GP Inference
Section 2 shows the main challenge of large-scale GP inference is the need solve a large linear system with the ill-conditioned matrix K + XT, which cannot be done satisfactorily with PCG or SDD. We introduce the SAP algorithm, which effectively addresses the conditioning challenges of GP inference. However, SAP alone does not address the scalability challenges of GP inference, motivating the development of ADASAP in Section 4.
3.1 Sketch-and-project
We formally present SAP in Algorithm SAP resembles randomized block coordinate descent, except that it employs subspace preconditioning. That is, it preconditions the subspace gradient

---

Algorithm SAP for KAW- = Y
Require: blocksize b, distribution of row indices D over [n]*, number of iterations maxx, initializa- tion Wo E Ruxxa,, averaging boolean tail_ lail_average
for t = 0, ,,..,.... -Tax - do
Sample row indices B of size from [n] according D
D, < 17 (K8B + XI) (Kg, Wt + NgW, - Y8) & Search direction; costs Obbnnts + bo)
W4+1 W Dr
Update parameters; costs (b
if tail average then
W++1 < 2/(t + 1) E)=
end if end for return Wr if taaiallvvvee avil_aveeege else WT
Approximate solution to KgW = Y
KB, W, + XgW, - Ys by the inverse of the subspace Hessian K88 + AI. Thus, unlike SDD, SAP includes second-order information in its search direction. In addition, we have also augmented vanilla SAP from Gower and Richtárik [2015] with tail averaging [Jain et al., 2018, Epperly et al., 2024a], which helps address the inherent noise in the iterates, and plays crucial role in our analysis.
SAP (assuming the number of right-hand sides Bhts 1 for simplicity) only costs Obbn + b²) per- iteration, which is significantly smaller than the O(n²) iteration cost of PCG, while being comparable with SDD, which costs (bn) per-iteration. Furthermore, PCG typically requires memory linear in to store the preconditioner, while SAP only incurs and additional O(b²) storage cost. Thus, SAP strikes a balance between PCG and methods like SDD. incorporates second-order information, like PCG, but does so modest increase in cost relative to SDD and SGD. By leveraging the theory of determinantal point processes [Kulesza and Taskar, 2012], we show that tail-averaged SAP improves dependence upon the condition number in the early phase of the convergence, which establishes it as a potential solution to the conditioning challenges discussed Section 2, particularly in the presence noisy data and generalization error.
3.2 Fast convergence along top- subspace for GP inference
Prior works [Dicker et al., 2017, Lin al., 2020, 2023] have shown that the components of the solution most relevant for learning lie along the dominant eigenvectors of K, which correspond to the dominant spectral basis functions of the reproducing kernel Hilbert space (RKHS) respect to X. Lin et al. [2023] shows that SGD enjoys improved early convergence along these directions, but attains a slow sublinear rate asymptotically. Motivated by these observations, we study the convergence rate of tail-averaged SAP along the dominant spectral basis functions, showing that it obtains a two-phase convergence rate: an initial condition number- frreeerrrfeeeeeeee sublinear rate, followed by an asymptotic linear rate where the condition number dependence is mitigated by the blocksize b. The proofs of the results in this section appear in Appendix B.3.
induced by kernel & with
Following Lin et al. [2023], we define the spectral basis functions of H. as
u (0) argmax [|((X]| | (u, (1) = for j < <1)
HEH.]|ulbu==
The spectral basis functions (1)), _2(2) ,.. defined in such way that each of them takes maximal values at the observation points while being orthogonal to the previous ones. Hence, to accurately estimate the posterior mean in the data-dense regions, is most important to minimize the error along the dominant basis functions.
We focus on posterior mean estimation, which corresponds to solving the linear system K,W W y using averaged SAP. Recall that, given a sequence of weight vectors Wi,.., We produced by sketch-and-Pe mnojjeeecct (which are vectors, not matrices, nhs 1), tail averaging obtains weight estimate by averaging the second half of this sequence. We use these weights to define an estimate of the GP posterior mean m,, which is given by m, = k(: X)W.
Theorem 3.2 our main theoretical result. provides a convergence guarantee for the SAP posterior mean estimate fy along the top-l spectral basis functions, where L is parameter that can be chosen

---

freely. Theorem 3.2 shows that while the estimate converges to the true posterior mean over the entire RKHS H. its initial convergence rate along the dominant basis functions even faster. Our theoretical results rely on the smoothed condition number, which we define below.
Definition 3.1 (Smoothed condition number). Let E Rax be a positive-semidefiite matrix with eigenvalues X 2 ..... 2 An. We define the smoothed condition number (bb,p) p) of A as
X X
((b,
b
The smoothed condition number controls both the SGD-style sublinear convergence rate and the PCG-style linear convergence rate in Theorem 3.2. Note that, for notational convenience, we 2b to denote the SAP blocksize in the below statement.
Theorem 3.2. Suppose we have a kernel matrix K E RR*", observations y E RR from a GP prior with likelihood variance À 2 O, and let c(...) denote the smoothed condition number of K + XI. Given any & E [nl, let proje denote orthogonal projection onto the spun of u()), ...., (l For any 1 S b S n/2, in O(nb²) time we can construct distribution over index subsets B of size 26 such that SAP (Algorithm 1) initialized at zero after t iterations satisfies
E (ma) (m))) S min [o(b,)) 2)
t (1 (1- 1
1/2
lelll
20(b,
This result provides composite convergence rate for SAP: the first rate, Só(b.)) is sublinear rate, 1/2
similar to SGD [Lin et al., 2023], while the second rate, (1 20(b,n)
is linear convergence rate. The linear rate is independent of the subspace norm, and thus, in particular, also applies for é t= n, i.e., when comparing posterior means within the entire space. However, for E < n, the sublinear rate wins out in the early phase of the convergence, since ((b, E) < ((b,
For ill-conditioned kernels with small prior variance X, the linear rate 1 - (/(p((b, n)) is close to 1 (as increases, the preconditioning effect becomes stronger, which improves the rate). On the other hand, the constant in the sublinear rate ob, E) is much smaller, since it takes advantage of the top- subspace norm, which makes it effectively independent of the conditioning of K.
If the SAP blocksize is proportional to l and K + AI exhibits any polynomial spectral decay, then obb E) (1), making the sublinear rate condition number-free free. Most popular kernels exhibit polynomial spectral decay [Ma and Belkin, 2017, Kanagawa et al., 2018]. Furthermore, polynomial spectral decay a standard assumption in the generalization analysis kernel methods in the fixed and random design settings [Caponnetto and Vito, 2007, Bach, 2013, Rudi et al., 2015]. To the best of our knowledge, our result the first of its kind, as previous comparable guarantees for stochastic methods [Lin et al., 2023, Proposition 1] still depend on the condition number. We note that our condition number-frre free rate is attained only when seeking moderately accurate estimates of the posterior mean (which is often the case due to the presence of noise in the data). When secking highly accurate estimates, the asymptotic (condition number dependent) linear rate will eventually take over. We illustrate these claims with the following corollary:
Corollary 3.3. Suppose that the matrix K exhibits polynomial spectral decay, i.e., (xKK) = O(i-9) for some 1. Then for any l [I,....n) and E (0, 1), choosing block size b = 4f can find in that with probability at least 0.99 satisfies [|proje (m proje(mm) 11 S e (|lprojo(mm))
in ( time.
δ 611² nl min {2(1+ (1+m(2,(K))
(()e(K)+ +
+ (1+
log
(K) X)
Remark 3.4. The Ō(n + nl)/e) runtime (phase one, attained for moderate accuracy e) is entirely independent of the condition number of K and the likelihood variance X. In Appendix B, we show that when À is sufficiently small (the highly ill-conditioned case), then the phase one complexity of sketch-anddpee and-project can be further improved by increasing the block size to attain O(nb² + (n² nb²) (e -1)0-//) runtime. In particular, when & & /n, we obtain meaningful convergence in o(n²) time, i.e., before reading all of the data. Thus, the top- subspace convergence of SAP can actually benefit from ill-conditioned kernels as the leading spectral basis functions become more dominant.
5

---

Although our theoretical results require tail averaging, we do not believe it is needed for good practical performance of sketch-and -nd--eeeectttctttt algorithms. Indeed, we run aaatchhand-projees with and without tail averaging in Appendix C.6, and we find that (i) tail averaging does not improve subspace convergence by substantial amount and (ii) tail averaging performs worse than not using tail averaging at larger blocksizes.
4 ADASAP: Approximate, Distributed, Accelerated SAP
We introduce ADASAP (Algorithm 2), scalable extension of SAP for GP inference. We first intro- duce the modifications made to SAP for scalability (Section 4. 4.1), before presenting ADASAP full
(Section 4.2). For detailed discussion of how ADASAP relates to prior work, please see Appendix A.
4.1
The key ingredients: approximation, distribution, and acceleration
We begin by discussing the essential elements of the ADASAP algorithm: (i) approximate subspace preconditioning, (ii) distribution, (iii) acceleration, and how they enhance performance.
4.1.1 Approximate subspace preconditioning
SAP enjoys significant improvements over PCG, it only requires Obb + b²) computation per- iteration. Unfortunately, SAP is limited in how large a blocksize b it may use, as factoring KBS + X/ to perform subspace preconditioning costs O(b)). This is problematic, Theorem 3.2 shows the convergence rate of SAP improves as increases.
To address this challenge, ADASAP draws inspiration from previous works [Erdogdu and Montanari, 2015, Frangella al., 2024b, Rathore al., 2025], which replace exact linear system solves in iterative algorithms with inexact solves based on low-rank approximations. ADASAP replaces Kgs in SAP with rank-r randomized Nyström approximation KES [Williams and Seeger, 2000, Gittens and Mahoney, 2016, Tropp al., 70. This strategy natural, as kernels exhibit approximate low-rank structure [Bach, 2013, Rudi et al., 2015, Belkin, 201 8). Indeed, Rathore et al. [2025] develops an approximate SAP solver for kernel ridge regression with one right-hand side, based on approximating KBS by KBB, and observes strong empirical performance.
ADASAP computes KE following the numerically stable procedure from Tropp al. [2017], which presented Appendix C. The key benefit is that (KRE + pl) can be applied to vectors (br) time (Appendix C), where >> 0 is damping parameter ensure invertibility. This reduces the cost compared the (b3) exact SAP update and allows ADASAP to use larger blocksizes: on the taxi dataset (n 3.31 10*), ADASAP uses blocksize b 1.65 105.
4.1.2 Distributed matrix- matrix products
SAP with approximate subspace preconditioning allows large blocksize b, but two bottlenecks remain: computing KB, W; and constructing the sketch K88s22 The former costs (bn), while the latter costs O(br). As matrix-mati ation can address these bottlenecks. By partitioning KBn across N workers GPUs, ColDistMatMar significantly reduces the time to compute the product KB. Similarly, RouDistMatMat for Kgg? reduces the time to compute the sketch. our largest dataset, taxi, using GPUs achieves 3.4x speedup (Fig. 4), which corresponds to near-perfect parallelism.
multiplication is embarrassingly parallel, distributed multi-GPP acceler-
4.1.3 Nesterov acceleration
Prior work has shown Nesterov acceleration [Nesterov, 1983] improves the convergence rate of SAP and approximate SAP [Tu al., 2017, Gower et al., 2018. Dereziński al., 2024, Rathore al., 2025]. Hence we use Nesterov acceleration in ADASAP to improve convergence.
4.2 ADASAP algorithm
The pseudocode for ADASAP is presented Algorithm 2. ADASAP uses tail averaging, approximate subspace preconditioning, distributed matrix-matrix products, and Nesterov acceleration. Assuming

---

Algorithm 2 ADASAP for KWW-
Require:
distribution of row indices D over [n]a, distribution Nyström sketching matri-
ces Daystrim over Rxr, number of iterations max. initialization Wo Ruxnn, workers (W...., WN workern ), averaging boolean tail average, acceleration parameters μ4,
Initialize acceleration parameters
+ 1- Va/u, y 7 1/1/pp, a at 1/(1+ Vo Wo. 20 W
for / 0, 1......... T maxx do
# Step I: Compute matrix-matrix product K W Sample row indices B of size from [n] according D
> Costs Obbrittss/ Nworkers)
KBn W + ColDistMatMat B, (W,..., WN. Wiwwaskaru)
b Algorithm
# Step II: Compute Nyström sketch KBB??
b Costs O(rb² /N workers)
Sample Q E R2xx from Dyyatrmm
KB6 RowDistMatMat(I.4 B, (W..... Wwocker
Algorithm
# Step III: Compute Nyström approximation and get stepsize U.S + BandNyssppx((sss] Q,,))
Costs O + br A Algorithm
Ps - USUT + (S., A)I
> Get preconditioner. Never explicitly formed!
78 < GetStepsize (Bs, KB5 XI)
> Algorithm 6
# Step IV: Compute updates using acceleration
Costs O Obrrat + 727
D. < P (P AgW, YYB) W+11, V4+. Ze+1 + NestAcc(W,, V1, Zi, Dr, nB, 7,
D Algorithm
tail average then
W411 < 2/(t + E)
end end for return Wr if tail. average else Wr
Approximate solution to KyW Y
perfect parallelism, ADASAP has iteration runtime of O Wookkers + b2r workers), significant improvment over the Omnaab + 6²) iteration time of SAP. The per-iteration runtime of ADASAP is comparable to distributed SDD-in other words, ADASAP effectively preconditions the problem, reducing the total iterations required, without significantly slowing down each iteration. Moreover, unlike SDD, ADASAP automatically sets the stepsize at each iteration, removing the need for expensive tuning. By default, we the acceleration parameters μ and À and n/b, respectively (as done in Rathore et al. [2025]), and find they yield excellent performance in Section 5,
4.2.1 Theory vs. Practice
Our theory in Section does not cover the approximate preconditioning, acceleration, uniform sampling that is used ADASAP. Despite the theory-practice gap between SAP (Algorithm 1) and ADASAP (Algorithm 2), believe the theory developed for SAP be extended ADASAP. This belief rooted Dereziński and Yang [2024], Dereziński et al. [2024], Rathore et al. [2025], which show that SAP methods still converge when using approximate preconditioning, acceleration, and uniform sampling. We leave this extension as direction for future research.
For simplicity, our experiments run ADASAP without tail averaging. Tail averaging is needed to estab- lish Theorem 3.2, but we not expect to yield significant practical improvements (Appendix C.6). This is in line with the SGD literature, where averaging leads to better theoretical convergence rates, but the last iterate delivers similar performance in practice [Shamir and Zhang, 2013, Johnson and Zhang, 2013].

---

Table RMSE and mean negative log-likelihood (NLL) obtained by ADASAP and competitors on the test set. The results are averaged over five 90%-10% train-test splits of each dataset. We bold a result it gets to within 0.01 of the best found RMSE or mean NL.L.
Dataset
yolanda
song
benzene
malonaldehyde 10
acsincome
houseelec
71 d
3.60 10
4.64
10
5.65 10
8.94
1.50 10
1.84 10 9
100
90
66
36
k
RBF
Matérn-3/2
Matérn-5/2
Matérn-5/2
RBF
Matérn-3/2
ADASAP
0.795
0.752
0.012
0.015
0.789
0.027
Test RMSE
ADASAP-I
0.808
0.782
0.168
0.231
0.795
0.066
SDD-1
0.833
0.808
0.265
0.270
0.801
0.268
SDD-10
0.801
0.767
0.112
Diverged
0.792
0.119
SDD-100
Diverged
Diverged
Diverged 0.141
Diverged 0.273
Diverged 0.875
Diverged
PCG
0.795
0.752
1.278
Test Mcan NLL
ADASAP
1.179
1.121
-2.673
-2,259
1.229
-2.346
A
1.196
1.170
-0.217
0.466
1.235
-2.185
SDD-1
1.225
1.203
0.531
0.903
1.242
-0.281
SDD-10
1.187
1.149
-0.762
Diverged
1.232
-1.804
SDD-100
Diverged
Diverged 1.121
Diverged -0.124
Diverged 0.925
Diverged 1,316
Diverged 2.674
PCG
1.179
Experiments
We present three sets of experiments showing ADASAP outperforms state-of-the-art methods for GP inference: (i) GP inference on large benchmark datasets (Section 5. (ii) GP inference on huge-scale transportation data analysis with 3- 10 samples (Section 5.2), and (iii) the Bayesian optimization task from Lin et al. [2023] (Section 5.3). We evaluate ADASAP against the following competitors:
ADASAP-I: variant of ADASAP where the subspace preconditioner P8 is set to the identity matrix. ADASAP-I is the same as accelerated randomized block coordinate descent.
The coordinate descent method introduced by Lin et al. [2024]. We tune SDD using three different stepsizes, and denote these variants by SDD-1, SDD-10, and SDD-100.
SDD:
PCG: A combination of block CG (O'Leary, 1980] with Nyström preconditioning [Frangella et al., 2023].
Our experiments are run in single precision on 48 GB NVIDIA RTX A6000 GPUs using Python 3.10 3.10, PyTorch 2.6.0 [Paszke al., 2019], and CUDA 12.5. We use 2, 3, and 1 GPU(s) per exper- iment in Sections 5.1 to 5.3, respectively. Code for reproducing our experiments is available at htt.raatrrrrresccccceeeeeneeeeee
Additional details are in Appendix D.
5.1 GP inference on large-scale datasets
We benchmark on six large-scale regression datasets from the UCI repository, OpenML, and sGDML [Chmiela al., 2017]. The results are reported in Table 1. ADASAP achieves the lowest RMSE and mean negative log- likellool (NLL) on each dataset; the NLL is computed using 64 posterior samples (via pathwise conditioning). SDD-10 is competitive with ADASAP yolanda and acsincome, but performs much worse on the other datasets. SDD is also sensitive to the stepsize: SDD-1 attains a larger RMSE and NLL than ADASAP on all datasets, while SDD-100 diverges all datasets. PCG obtains the same RMSF and NLL as ADASAP on yolanda and song, but its performance degrades for larger datasets. Additionally, Fig. 2 shows ADASAP achieves the lowest RMSE and NLLI throughout the optimization process, demonstrating its efficiency with respect to wall-clock time.
5.2 Showcase: Transportation data analysis with > 3 10 samples
To demonstrate the power of ADASAP on huge-scale problems, we perform GP inference on a subset of the taxi dataset
with 3.31 10 samples

---

2201
bruvee in 5:15 (P))
traieaaiiinn (e- 48.110)
210
Trewedee (n- 144 lif)
100
022
000
cankmabbdodee 10%
12-
912
1200
Figure 2: Performance of ADASAP and competitors on RMSE and mean NLI., as function of time, for benzene, malonaldehyde, and houseelee. The solid curve indicates mean performance over random splits of the data; the shaded regions indicate the range between the worst and best performance over random splits of the data, ADASAP outperforms the competition.
tasd (n 3.31 10%)
RM55
0.8
0.7
Thxt
0.6
0.5
L.8% speedup
0.4
SDD
ADAGAP-1
ADABAP
SDD--10
Figure 3: Comparison between ADASAP and competitors transportation data analysis. ADASAP attains the lowest RMSE and obtains 1.8x speed up over the second-best method, SDD-10. SDD-100 diverges and PCG runs out of memory, so they do not appear in the figure.
and dimension d the task predict taxi ride durations in New York City. To the best of knowledge, this is the first time that full GP inference has been scaled to dataset of this size. Due memory constraints, we are unable to compute 64 posterior samples as done in Section 5.1 for computing NLL, so we only report RMSE.
The results shown in Fig. 3. ADASAP obtains the lowest RMSE out of all the methods. Once again, SDD is sensitive to the stepsize: SDD-1 obtains an RMSE of 0.60, as opposed to ADASAP, which obtains RMSE of 0.50, while SDD-100 diverges. SDD-10 obtains an RMSE of 0.52, which is similar to that of ADASAP. However, ADASAP reaches the RMSE attained by SDD-10 in 45% less time than SDD-10, which translates memory, as the memory required for storing the sketch in single precision is 3.31 10 100 ~ 130 GB, which exceeds the 48 GB of memory in the A6000 GPUs used in experiments.
difference of 14 hours of runtime! Furthermore, PCG runs out of
5.3 Large-scale Bayesian optimization
run ADASAP competitors on a variant of the synthetic large-scale Bayesian optimization tasks from Lin et al. [2023]. These Bayesian optimization tasks consist of finding the maximum of black

---

box functions f: [0. _ - R sampled from Matérn-31 3/2 Gaussian process. We use two different lengthscales (2.0 and 3.0) and 5 random functions per lengthscale. To avoid misspecification, we set the kernel of each model to match that of the black box function.
The results shown in Table 2. ADASAP makes the biggest improvement over the random search baseline. As we have seen in the other experiments, SDD is sensitive to the stepsize: SDD-10 and SDD-100 make no progress. PCG also makes less progress than ADASAP across both lengthscales.
Table 2: Percentage improvement over random search for Bayesian optimization, averaged over five seeds. SDDD-10 and SDD-1100 provide improvement over random search because they unstable.
Lengthscale 2.0
Lengthscale = 3.0
ADASAP
10.42
13.86
Improv. (%)
ADASAP-I
7.04
11.27
SDD-1
6,50
11.17
SDD-10
0.00
0.00
SDD-100
0.00
0.00
PCG
0.13
5.54
6 Conclusion
We introduce ADASAP, an approximate, distributed, accelerated sketch-and-project method for GP inference. We demonstrate that ADASAP outperforms state-of-the-art GP inference methods like PCG and SDD on large-scale benchmark datasets, huge-scale dataset with > 3 10 samples, and lcage-scale Bayesian optimization. Moreover, we show that SAP-style methods are theoretically principled for GP inference- we prove that SAP is the first efficient, condition number-frn free algorithm for estimating the posterior mean along the top spectral basis functions. Future work should extend the theoretical results for SAP to ADASAP and investigate ADASAP in lower precision (e.g., float16).
Acknowledgments and Disclosure of Funding
We would like to thank Jihao Andreas Lin for helpful discussions regarding this work. PR, ZF, SF, and MU gratefully acknowledge support from the National Science Foundation (NSF) Award IIS-2233762. the Office of Naval Research (ONR) Awards N000142212825,
0000142412306 and 07001423122000 the Alfred Sloan Foundation, and from IBM Research founding member of Stanford Institute for Human-centered Artificial Intelligence (HAI). MD and SG gratefully acknowledge support from NSF Award CCF-2338655.
References
Ahmad Abdelfattah, Hartwig Anzt, Erik G Boman, Erin Carson, Terry Cojean, Jack Dongarra, Alyson Fox, Mark Gates, Nicholas J Higham, Xiaoye S Li, Jennifer Loe, Piotr Luszczek. Srikara Pranesh, Siva Rajamanickam, Tobias Ribizel, Barry Smith, Kasia Swirydowicz, Stephen Thomas, Stanimire Tomov, Yaohung M Tsai, and Ulrike Meier Yang. A survey of numerical linear algebra methods utilizing mixed-precisionn arithmetic. The International Journal of High Performance Computing Applications, 35(4):344-36) 2021.
Sivaram Ambikasaran, Daniel Foreman Mackey, Leslie Greengard, David W Hogg, and Michael O'Neil. Fast direct methods for Gaussian processes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(2):252-265, 2015.
Nima Anari and Michat Dereziński. Isotropy and log- log-concave polynomials: Accelerated sampling and high-precision counting of matroid bases. In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (F (FOCS), 2020.
Nima Anari, Yang P Liu, and Thuy-Duong Vuong. Optimal sublinear sampling of spanning trees and determinantal point processes via average-case entropic independence. SIAM Journal on Computing, pages FOCS22-93. 2024.
10

---

Haim Avron, Kenneth L Clarkson, and David P Woodruff. Faster kernel ridge regression using sketching and preconditioning. SIAM Journal on Matrix Analysis and Applications, 38(4):1116- 1138, 2017.
Francis Bach. Sharp analysis of low-rank kernel matrix approximations. In Conference on Learning Theory, 2013.
Mikhail Belkin. Approximation beats concentration? An approximation view inference with smooth radial kernels. In Conference On Learning Theory, 2018.
Daniele Calandriello, Michal Dereziński,, and Michal Valko. Sampling from a k-DPP without looking at all items. In Advances in Neural Information Processing Systems, 2020.
Andrea Caponnetto and Emesso DeVito. Optimal rates for the regularized least-squares algorithm. Foundations of Computational Mathematics, 7:331-368, 2007.
Benjamin Charlier, Jean Feydy, Joan Alexis Glaunes, François-David Collin, and Ghislain Durif. Kernel operations on the GPU, with autodiff, without memory overflows. Journal of Machine Learning Research, 22(74):1-6, 2021.
Yifan Chen, Houman Owhadi, and Florian Schäfer. Sparse Cholesky factorization for solving nonlinear PDEs via Gaussian processes. Mathematics of Computation, 94(353): 2025.
Li-Fang Cheng. Bianca Dumitrascu, Gregory Darnell, Corey Chivers, Michael Draugelis, Kai Li, and Barbara E Engelhardt. Sparse multi-output gaussian processes for online medical time series prediction. BMC medical informatics and decision making, 20:1-23, 2020.
Stefan Chmicla, Alexandre Tkatchenko, Huziel E. Sauceda, Igor Poltavsky, Kristof Schütt, and Klaus-Robert Müller. Machine learning of accurate energy-conserving molecular force fields. Science Advances, 365)s1 2017.
Kurt Cutajar, Michael Osborne, John Cunningham, and Maurizio Filippone. Preconditioning kernel matrices. In Proceedings of the 33rd International Conference on Machine Learning, 2016.
Michat Dereziński and Michael W Mahoney. Determinantal point processes in randomized numerical linear algebra. Notices of the American Mathematical Sociery, 68(1):34 45, 2021.
Michat Dereziński and Elizaveta Rebrova. Sharp analysis sketch-and-projee methods via a connection to randomized singular value decomposition. SIAM Journal on Mathematics of Data Science, 6(1):127-153, 2024.
Michat Dereziński and Jiaming Yang. Solving dense linear systems faster than via preconditioning. In Proceedings of the 56th Annual ACM Symposium on Theory of Computing, 2024.
Michat Dereziński, Daniel LeJeune, Deanna Needell, and Elizaveta Rebrova. Fine-graines analysis and faster algorithms for iteratively solving linear systems. arXiv preprint wrXiv:2405.05816 2024.
Michat Dereziński, Deanna Needell, Elizaveta Rebrova, and Jiaming Yang. Randomized Kaczmarz methods with Sayood-KKyloo convergence. arXiv preprint arXiv:2501.1167. 2025.
Mateo Díaz, Ethan N Epperly, Zachary Frangella, Joel A Tropp, and Robert Webber. Robust, randomized preconditioning for kernel ridge regression. arXiv preprint arXiv:2304. 124445, 2023.
Lee H Dicker, Dean P Foster, and Daniel Hsu. Kernel ridge principal component regression: Minimax bounds and the qualification of regularization operators. Electronic Journal of Statistics, 11:1022-1047, 2017.
Ethan N Epperly, Gil Goldshlager, and Robert Webber. Randomized Kaczmarz with tail averaging. arXiv preprint 2rXiv:2411.19877 2024a.
Ethan N. Epperly, Joel A. Tropp, and Robert J. Webber. Embrace rejection: Kernel matrix approxi- mation by accelerated randomly pivoted Cholesky. arXiv preprint 2xX1v:2410.03969, 2024b.

---

Murat A Erdogdu and Andrea Montanari. Convergence rates of sub-sampled Newton methods. In Advances in Neural Information Processing Systems, 2015.
Zachary Frangella, Joel A. Tropp, and Madeleine Udell. Randomized Nyström preconditioning.
SIAM Journal on Matrix Analysis and Applications, 44(2):118-752, 2023.
Zachary Frangella, Pratik Rathore, Shipu Zhao, and Madeleine Udell. Promise: Preconditioned
stochastic optimization methods by incorporating scalable curvature estimates. Journal of Machine 2024a. Learning Rexearch, 25(346):1-57, 20244.
Zachary Frangella, Pratik Rathore, Shipu Zhao, and Madeleine Udell. Sketchysgd: reliable stochastic optimization via randomized curvature estimates. SIAM Journal on Mathematics of Data Science, 6(4):1173-1204 2024b.
Peter Frazier and Jialei Wang. Bayesian optimization for materials design. Information science for materials discovery and design, pages 45-75, 2016.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPyTorch: Blackbox matrix matrix-matrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems, 2018.
Alex Gittens and Michael W Mahoney, Revisiting the Nyström method for improved large-scale machine learning. The Journal of Machine Learning Research, 17(1):3977-4041, 2016.
Robert Gower, Filip Hanzely, Peter Richtarik, and Sebastian U Stich. Accelerated stochastic matrix inversion: General theory and speeding up BFGS rules for faster second-order optimization. In Advances in Neural Information Processing Systems, 2018.
Robert Gower, Dmitry Kovalev, Felix Lieder, and Peter Richtárik. RSN: randomized subspace Newton. In Advances in Neural Information Processing Systems, 2019.
Robert M Gower and Peter Richtárik. Randomized iterative methods for linear systems. SIAM Journal on Matrix Analysis and Applications, 36(4): 1664)):660-1600, 2015.
Philip Greengard, Manas Rachh, and Alex H. Barnett. Equispaced Fourier representations for efficient Gaussian process regression from billion data points. SIAM/ASA Journal on Uncertainty Quanrification, :3(1):63 (11:))63-880 89, 2025.
Filip Hanzely, Nikita Doikov, Yurii Nesterov, and Peter Richtarik. Stochastic subspace cubic Newton method. In Proceedings of the 37th International Conference on Machine Learning, 2020.
James Hensman, Nicolò Fusi, and Neil D Lawrence. Gaussian processes for big data. In Proceedings of the Twetty-Nith Conference on Uncertainty in Artificial Intelligence, 2013.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, 2015.
James Hensman, Nicolas Durrande, and Arno Solin. Variational Fourier features for Gaussian processes. Journal of Machine Learning Research, 18(151):1-52, 2018.
José Miguel Hernández-Loban James Requeima, Edward O. Pyzer-Knapp. and Alán Aspuru- Guzik. Parallel and distributed Thompson sampling for large-scale accelerated exploration of chemical space. In Proceedings of the 34th International Conference on Machine Learning, 2017.
Nicholas J Higham. Accuracy and stability of numerical algorithms. SIAM, 2002.
Prateek Jain, Sham M Kakade, Rahul Kidambi, Praneeth Netrapalli, and Aaron Sidford. Parallelizing stochastic gradient descent for least squares regression: mini-batching, averaging, and model misspecification. Journal of Machine Learning Research, 18(223):1 2018.
Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems, 2013.
12

---

Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur. Gaus- sian processes and kernel methods: review on connections and equivalences. arXiv preprint arXiv:/807.02558, 2018.
Jacek Kuczyński and Henryk Wožniakowski. Estimating the largest eigenvalue by the power and lanczos algorithms with random start. SIAM Journal on Matrix Analysis and Applications, 13(4): 1094-1122 1122, 1992.
Alex Kulesza and Ben Taskar. Determinantal point processes for machine learning, Foundations and Trends@ in Machine Learning, 5(2-3): 2012.
D. Leventhal and A. S. Lewis. Randomized methods for linear constraints: Convergence rates and conditioning. Mathematics of Operations Research, 35(3):641-65) 2010.
Jihao Andreas Lin, Javier Antoran, Shreyas Padhy, David Janz, José Miguel Bermández-Lobate and Alexander Terenin. Sampling from Gaussian process posteriors using stochastic gradient descent. In Advances in Neural Information Processing Systems, 2023.
Jihao Andreas Lin, Shreyas Padhy, Javier Antoran, Austin Tripp. Alexander Terenin, Csaba Szepes- vari, José Miguel Hernández-Lobate and David Janz. Stochastic gradient descent for Gaussian processes done right. In International Conference on Learning Representations, 2024.
Junhong Lin, Alessandro Rudi, Lorenzo Rosasco, and Volkan Cevher. Optimal rates for spectral algorithms with least-squares regression over Hilbert spaces. Applied and Computational Harmonic Analysis, 48(3):868-88) 890, 2020.
Siyuan Ma and Mikhail Belkin. Diving into the shallows: A computational perspective on large-scale shallow learning. In Advances in Neural Information Processing Systems, 2017.
Per-Gunnal Martinsson and Joel A Tropp. Randomized numerical linear algebra: Foundations and algorithms. Acta Numerica, 29:403-572, 2020.
Ian McDowell, Dinesh Manandhar, Christopher M Vockley, Amy K Schmid, Timothy E Reddy, and Barbara E Engelhardt. Clustering gene expression time series data using an infinite gaussian process mixture model. PLOS computational biology, 14(1)ce 1005896, 2018.
Victor Minden, Anil Damle, Kenneth L Ho, and Lexing Ying. Fast spatial Gaussian process maximum likelihood estimation via skeletonization factorizations. Multiscale Modeling & Simulation, 15(4): 1584-1611, 2017.
Mojmir Mutny, Michal Derezinski, and Andreas Krause. Convergence analysis of block coordinate algorithms with determinantal sampling. In International Conference on Artificial Intelligence and Statistics, 2020.
Arkadi S Nemirovski and David B Yudin. Problem complexity and method efficiency in optimization. Wiley-Interscience, 1983.
Nesterov. A method of solving a convex programming problem with convergence rate o (1/k** * 2). Doklady Akademii Nauk SSSR, 269(3), 1983.
Dianne P O'Leary. The block conjugate gradient algorithm and related methods. Linear Algebra and its Applications, 29:293-322, 1980.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vetito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An imperative style, perfforrffme deep learning library. In Advances in Neural Information Processing Systems, 2019.
Zheng Qu and Peter Richtárik. Coordinate descent with arbitrary sampling i: Algorithms and complexity. Optimization Methods and Software, 31(5):829-88) 857, 2016.
Zheng Qu, Peter Richtarik, Martin Takac, and Olivier Fercoq. SDNA: Stochastic dual Newton ascent for empirical risk minimization. In Proceedings of the 33rd International Conference on Machine Learning, 2016.
13

---

Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, 2007.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian processes for machine learning.
The MIT Press, 11 2005.
Pratik Rathore, Zachary Frangella, Jiaming Yang, Michal Dereziński, and Madeleine Udell. Have ASkotch: A neat solution for large-scale kernel ridge regression. arXiv preprint 2407.10070, 2025.
Peter Richtárik and Martin Takáč. Iteration complexity of randomized coock-coordinate descent methods for minimizing composite function. Mathematical Programming, 144(1):1-38, 2014.
Anton Rodomanov and Dmitry Kropotov. A randomized coordinate descent method with volume sampling. SIAM Journal on Optimization, 30(3):1878-1044 2020.
Alessandro Rudi and Lorenzo Rosasco. Generalization properties of learning with random features. In Advances in Neural Information Processing Systems, 2017.
Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco. Less is more: Nyström computational regularization. Advances in neural information processing systems, 28, 2015.
Ohad Shamir and Tong Zhang. Stochastic gradient descent for non-smooth optimization: Convergence results and optimal averaging schemes. In Proceedings of the 30th International Conference on Machine Learning, 2013.
Stefan Steinerberger. Randomized Kaczmarz converges along small singular vectors. SIAM Journal on Matrix Analysis and Applications, 42(2))608-615, 2021.
Thomas Strohmer and Roman Vershynin. A randomized Kaczmarz algorithm with exponential convergence. Journal of Fourier Analysis and Applications, 15(2):222-278, 2009.
Michalis Titsias Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statistics, 2009.
Lloyd N Trefethen and David Bau. Numerical linear algebra. SIAM, 1997.
Joel A Tropp. Alp Yurtsever, Madeleine Udell, and Volkan Cevher. Fixed-rank approximation of a posidive-semidefiiite matrix from streaming data. In Advances in Neural Information Processing Systems, 2017.
Stephen Tu, Shivaram Venkataraman, Ashia C. Wilson, Alex Gittens, Michael I. Jordan, and Benjamin Recht. Breaking locality accelerates block Gauss-Seide In Proceedings of the 34th International Conference on Machine Learning, 2017.
Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gordon Wilson. Exact Gaussian processes on million data points. In Advances in Neural Information Processing Systems, 2019.
Christopher Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. In Advances in Neural Information Processing Systems, 2000.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In Proceedings of the 32nd International Conference on Machine Learning, 2015.
James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Deisenroth. Efficiently sampling functions from Gaussian process posteriors. In Proceedings of the 37th International Conference on Machine Learning, 2020.
James T. Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth. Pathwise conditioning of Gaussian processes. Journal of Machine Learning Research, 22(105):1-47, 2021.
Kaiwen Wu, Jonathan Wenger, Haydn T Jones, Geoff Pleiss, and Jacob Gardner. Large-scale Gaussian processes via alternating projection. In International Conference on Artificial Intelligence and Statistics, 2024.
14

---

A Related Work
We review the literature on Gaussian process (GP) inference and bketch-and-projeer which are the two key areas that serve as foundations for ADASAP.
A.1 GP inference
Naive GP inference based on solving the linear systems in (1) exactly and generating posterior samples via (3) has cubic cost in number of training and test points, making it prohibitively expensive for datasets whose size exceeds n 10 Given its prominence in scientific computing and machine learning, much work has been to done scale GP inference the big-data setting. Prior work in this area can be (roughly) divided into four approaches: (i) PCG-based inference, (ii) SGD-based inference, (iii) low-dimensional inference, and (iv) variational inference.
A.1.1 PCG-based inference
natural approach to scaling GP inference is to replace the exact linear system solves by direct method with approximate linear system solves by an iterative method. Gardner et al. [2018] showed that by utilizing batched PCG and GPU acceleration, GP inference could be scaled to datasets of size 10 Wang et al. [2019] further showed that by distributing the kernel matrix across multiple GPUs, that inference could be applied to datasets of size n ~ 10 Unfortunately, these approaches still have shortcomings. While Gardner al. [2018]. Wang et al. [2019] show that exact GPs can be scaled to datasets with n ~ 108 by leveraging PCG and GPUs, they do not address the cubic complexity of generating posterior samples in (3). Moreover, when prohibitively expensive from both a memory and computational standpoint, as most preconditioners require (nr) storage [Cutajar et al., 2016, Avron et al., 2017, Díaz al., 2023] and each PCG iteration costs (n²²)
10 PCG can become
A.1.2 SGD-based inference
To address the limitations of PCG- inference, Lin et al. [2023] proposed to replace PCG with SGD. This reduces iteration cost from O(n²) (bb)). where b is the gradient batchsize. addition, they address the challenge of the cubic cost of generating posterior samples by adopting the pathwise conditioning strategy [Wilson et al., 2020, 2021]. Lin et al. [2023] show SGD outperforms PCG and Stochastic Variational Gaussian Processes (SVGP) [Hensman et al., 2013] on various posterior sampling tasks. Lin et al. [2024] further improved upon SGD by introducing Stochastic Dual Descent (SDD). SDD is a coordinate descent algorithm that leverages the dual formulation of the KRR problem to reduce the dependence on the conditioning of from is( (((K + XT)) to x(K + A7). In addition, SDD incorporates momentum and geometric averaging to further enhance performance. Compared to SGD, PCG, and SVGP, Lin et al. [2024] shows that SDD yields the best empirical performance.
While Lin et al. [2023, 2024] address the scaling limitations of PCG and the issue of efficiently sampling from the posterior, there is still room for improvement. In terms of theory, Lin et al. [2023] only establishes that SGD is guaranteed to converge at reasonable rate along the top eigenvectors of K. Convergence along the small eigenvectors is extremely slow, as the rate in Lin et al. [2023] along the t eigenvector depends upon 1/ 1/v For SDD, the convergence rate is unknown, as Lin et al. [2024] provides no convergence analysis. Aside from unsatisfactory convergence analyses, it is known from optimization thcory, that in the worst case, the convergence rate of SGD and SDD is coeeeioooooiii
controlled by the condition number [Nemirovski and Yudin, 1983]. As ill ill-conditioning is defining property of kernel matrices, this means SDD and SGD will converge slowly. Another price incurred by SDD and SGD for their fast performance is the presence of additional hyperparameters. In particular, the stepsize for these methods must be set appropriately to obtain satisfactory convergence. Tuning the stepsize, while possible, can become expensive for large-scale inference datasets like taxi. This is in contrast to PCG, which is stepsize- freee
ADASAP enjoys the scalability benefits of SDD and SGD, but without their limitations. It only operates on a batch of data at each iteration, so it avoids the O(n²) iteration cost of PCG. By incorporating subspace preconditioning, ADASAP addresses the ill-conditioning of the kernel matrix to achieve fast
15

---

convergence. ADASAP also comes with reliable default hyperparameters, which obviates the need for costly hyperparameter tuning.
A.1.3 Low-dimensional inference
When the dimension d of the data matrix X is low (d I 3), GP inference can be made far more efficient by exploiting the structure of the kernel matrix and the low-dimensionality of the data. Common approaches in this vein include rank-structured matrices [Ambikasaran et al., 2015, Minden et al., 2017], sparse Cholesky factorization [Chen et al., 2025], and kernel interpolation [Wilson and Nickisch, 2015, Greengard et al., 2025]. In most cases these methods enable (n) training time and inference, a massive improvement over standard GP inference methods. This linear time complexity has allowed GP training scale datasets with n large 10 [Greengard al., 2025]. Unfortunately, the complexity of these algorithms grows exponentially with the dimension, making them impractical when d >> 3.
A.1.4 Variational inference
To address the challenge of large-scale GP inference, many prior works have focused on developing scalable approximate GP methods. The most popular methods are based on variational inference [Titsias, 2009, Hensman et al., 2013, 2015, 2018). 1J.. The most well-known approach is SVGP [Hensman at., 2013], which selects a set of m inducing points by maximizing an evidence lower bound (ELBO). The use of inducing points leads to only need to work with an m kernel matrix instead of the full n x n kernel matrix. This leads reduces the cost training to O(um² + m²)), a significant improvement over exact inference methods when m can be taken to be much smaller than n. While SVGP and related algorithms significantly reduce the cost of GPs, they generally exhibit worse performance inference than methods which make use of the full data [Wang et al., 2019, Lin et al., 2023, 2024]. Thus, it is advisable to use exact inference whenever possible.
A.2 Sketch-and-project
ADASAP builds off the brrich-and-projeer framework for solving consistent linear systems. Sketch- and-project was first formulated in Gower and Richtárik [2015], who showed that randomized Kaczmarz [Strohmer and Vershynin, 2009], randomized coordinate descent [Leventhal and Lewis, 2010), and randomized Newton [Qu et al., 2016), are all special cases of sketch-and- also established a linear convergence rate for sketch-and-project, which is controlled by the smallest eigenvalue of an expected projection matrix. However, Gower and Richtárik [2015] 5] were unable to provide fine-grained analysis of kketch-and-project in terms of the spectral properties of the linear system, except for few special cases.
They
Follow-up work has combined Aketch-and-projer with techniques such Nesterov acceleration [Tu et al., 2017] and extended sketch-and -ad--andd--jjeeett to quasi-Newtor and Newton type methods Gower et al. [2018, 2019], Hanzely al. [2020]. Despite these extensions, sharp convergence analysis of sketch-and and-jprjecct remained elusive, even in the original setting of consistent linear systems.
Recent work has finally given sharp analysis of kketch-and-projed Using powerful tools from random matrix theory, Dereziński and Rebrova [2024] provided the first sharp analysis of sketch-and- project when the sketching matrix is sub-Gaussiar or sparse. In particular, they were the first to prove that if the matrix exhibits a favorable spectral decay profile, then seecch-and-projett exhibits condition number-free convergence. Dereziñiski and Yang [2024] improve upon these results, showing that by employing a pre-processing step, one can avoid using expensive sketching matrices at each step, and simply sample rows uniformly. They accomplish this by leveraging tools from determinantal point process theory [Mutny et al., 2020, Rodomanov and Kropotov, 2020, Anari et al., 2024]. Dereziński et al. [2024] further refines these results by incorporating Nesterov acceleration and demonstrating an improved computational complexity relative to a fine-grained analysis of PCG-that is, an analysis of PCG that goes beyond the worst-case ((yR) R)) iteration complexity [Trefethen and Bau, 1997].
In addition to improvements in the analysis of exact sketch-and anroo-ee for linear systems, significant strides have been made when approximate subspace preconditioning is employed. Dereziński and Yang [2024], Dereziáski et al. [2024] both use PCG to approximately solve the linear system defined by the subspace Hessian each iteration of SAP, while Dereziński et al. [2025] augment that system with Tikhonov regularization to further improve its conditioning. Rathore et al. [2025] consider an
16

---

approximate bertch-and-projees method for solving kemel ridge regression with one right-hand side. They replace KBS + AI in SAP with KES + pl, where KES is a randomized Nyström approximation of K88- They establish linear convergence for the method, with a rate that is comparable to exact SAP when the kernel matrix exhibits an appropriate rate of spectral decay.
Among the above-mentioned works, ADASAP is closest to Rathore al. [2025], but with significant algorithmic differences such as tail averaging, handling multiple right-hand sides, and distributing bottleneck operations over multiple devices. On the theoretical side, our analysis focuses on the convergence along the top Æ-subspace using exact SAP with -il-averaging. In contrast, Rathore et al. [2025] focus on convergence of approximate the entire space, and as a result, they are only able to show fast convergence guarantee under certain strong conditions on the spectral decay (effectively limiting the condition number of the problem). Consequently, the global convergence analysis of Rathore et al. [2025] does not explain the rapid initial progress that SAP-style algorithms make on test error. We believe that our two-phase analysis consisting of (i) fast condition number-free sublinear convergence, followed by (ii) slower global convergence rate dependent upon the subspace condition number, better captures this phenomenon. In particular, our result is the first step in developing systematic analysis of the generalization properties SAP-style algorithms.
(with and without acceleration) over
To the best of our knowledge, our work is the first to systematically investigate the convergence rate of SAP with tail averaging along the top eigenspace. Epperly et al. [2024a] combines tail averaging with randomized Kaczmarz for solving inconsistent linear systems, but they do not investigate the convergence rate along the top spectral subspaces. Steinerberger [2021], Dereziński and Rebrova [2024] look at the convergence rate of the expected iterates along particular eigenvector. However, as noted in Section 3.2, this does not lead to a meaningful convergence rate for the iterates produced by the algorithm, and this also fails to capture the behavior over an entire subspace. Thus, our analysis provides the first concrete characterization of the convergence rate of tail-averaged SAP along the top eigenvectors.
B Proofs of the Main Results
In this section, we give the detailed proofs for the main theoretical results given Section 3.2, namely Theorem 3.2 and Corollary 3.3. We start with informal overview of the analysis in Section B.1, followed by a formal convergence analysis of SAP for solving positive-definite linear systems in Section B.2, and finally adapting the results to GP posterior mean estimation in Section B.3.
Throughout, we adopt the shorthand A == A + AI, where A square matrix, À R. and denotes the identity matrix of the same size as A.
B.1 Overview of the analysis
Let w. KT²u y E IR" be the solution of the linear system defined by K₂ and y, an let use Wy E TR" to denote the iterates produced by sketch-and-projeet when solving that system. To convert from the posterior mean error the Euclidean norm, we first write the spectral basis functions as (())=
Vxs kk,,, X)u;, where (A,, U;) is the 8th -laggett eigenpair of K. Thus, for any w R and m k(, X)w E H. we have [Iproje(m - m.,) I = 12 w.)ll s [[Qxk)) (w w.)ll- where Q+ Di=1 ".aa is the projection onto the top-l subspace of K.
Having converted to the Euclidean norm, we express the beecch-and-projeee update as a recursive formula for its (scaled) residual vector as follows:
A+t1 K 1/2 (we+1 w.) = [I) (w. -w+) (T -gg)
where ITg K 1/2 T (KB,B X a random projection defined by B. To characterize the expected convergence of sketch-and- projjjee,,,,j,,,,,,,, we must therefore control the average-case properties of Ng. We achieve this by relying on row sampling technique known as determinantal point processes (DPPs) [Kulesza and Taskar, 2012, Dereziński and Mahoney, 2021]. A b-DPP(Kx) is defined as distribution over size index sets C [n]] such that Pr(B] xx det K + A/). DPPs can sampled from efficiently [Calandriello et al., 2020, Anari et al., 2024]: after an initial preprocessing cost of Ö(n²), can produce each DPP sample of size b in Õ(b3) time.
17

---

Prior work [Dereziński and Yang, 2024] has shown that when B is drawn according to 2 then the expectation of 1gg has the same eigenbasis as K. Concretely, we can show that II := E 11g VAVT with A = diag(>i,..., An). where X 2 I+ and V consists of the eigenvectors of K. This leads convergence guarantee for the expected residual vector:
[ A++1 = 11Qu II) E A ll S (1 (1-x))]1Q X
where we used that IT, Xx. and Qe commute with each other. This would suggest that we should obtain a fast linear rate (1 1+olb.z) ) for the sketch-and srrch-and-pren proj-proect iterates in the top-& subspace norm. However, it turns out that the actual residual does not attain the same convergence guarantee as its expectation because, unlike IT, the random projection πg does not commute with K2x or Qe. Indeed,
B|IQr4++1 2 A [(T QT) TB) A
and even though the matrix (1 - Hg QQ Igg) is low-rank, its expectation may not be, since IIs does not commute with Qe. This means that we have no hope showing a linear convergence rate along the top-l subspace. Nevertheless, we are still able to obtain the following bound:
E (1 Hz))Q(( II)) = Qe(1 - (i) + (I - Qc) Iī).
This results in bias-varianes decomposition: the first term is the bias, which exhibits fast convergence in top-l subspace, while the second term is the variance, which accounts for the noise coming from the orthogonal complement subspace. We then tail averaging decay this noise, which allows us to benefit from the fast convergence of expected iterates through a sublinear rate 6l0.2
Finally, to attain the asymptotically faster linear rate for the residual vectors, we observe that every top-f subspace norm is upper-bounded by the full norm: [| S ||||||. Thus, we can effectively repeat the above analysis with = n (since Qn = I), in which case the variance term becomes zero,
and we recover linear rate of the form (1 I+o(b,n)
B.2 Convergence along top-L dimensional subspace for positive-definite matrices
In this section, we consider the general problem of solving a linear system Aw = y for an n n positive-definite matrix A. We provide theoretical results providing the fast convergence guarantees for the iterate sequence generated by SAP along the top-f dimensional subspace of A, with Theo- rem B.3 being the main result of this section. We use Theorem B.3 in the next section to derive the subspace convergence result for posterior mean in GP inference, proving Theorem 3.2.
Preliminaries and notation We start by introducing the notation and some results from existing literature for SAP and DPPs. Let A be an n X n symmetric positive-definite matrix and let VDVT be its eigendecomposition, where D is diagonal with entries \ 2 2 Xn For other matrices we use (M) to denote the eigenvalues of M. For any vector v E RR. define ||e||a == Vo..u. Let V22 E R' consist of the top-& orthonormal eigenvectors of and Q: = V.VV"" be the orthogonal projector corresponding to the top-& eigenvectors of A. We rely on subsampling from fixed- size DPPs [Kulesza and Taskar, 2012 Dereziński and Mahoney, 2021] for our theoretical results. Here, for convenience, we consider the block size be 2b, in contrast to b considered in Algorithm 1. For a fixed and a positive-definite matrix A, we identify 21-DPP(A) A) as a distribution over subsets of [I....,n} of size 26, where any subset has probability proportional the corresponding principal submatrix of A, i.e., det(AAs The subsampling matrix S E R 2bxx corresponding to the set B is defined as a matrix whose rows are the standard basis vectors associated with the indices in B. For notational convenience, we will sometimes say that S is drawn from 20-DPP(A), or S ~ 26-DPP(A), keeping the index set B implicit.
We use the following result from the literature that provides an exact characterization of the eigenvec- tors of the expected projection matrix that arises in the analysis of sketch-and aachh-nn-ppaee and-project methods, when using subsampling from fixed-size DPPs.
Lemma B.1 (Expected projection under 26-DPP(A), adapted from Lemma 4.1 [Dereziński and Yang. 2024]). Let A E Roxn be symmetric positive-defiite matrix with eigenvalues X 2 X2 2 2 X.. and let 1 E b n/2 be fixed. Let 2b-DPP A). Then, we have
5/2 (SAS")) S y
=VD'V
where D' a diagonal matrix with jth diagonal entry is lower bounded by
d Eiss
18

---

As exact sampling from DPPs is often expensive and requires performing operations costly as performing the eigendecomposition of A [Kulesza and Taskar, 2012], numerous works have looked at approximate sampling from these distributions [Calandriello et al., 2020, Anari et al., 2024]. Since we exploit the exact characterization of the eigenvectors of the expected projection matrix, our theory requires a very accurate DPP sampling algorithm. Fortunately, the Markov Chain Monte Carlo tools developed by Anari and Dereziński [2020], Anari et al. [2024] allow near-exact sampling from fixed-size DPPs (in the sense that the samples are with high probability indistinguishable from the exact distribution), while offering significant computational gains over exact sampling.
Lemma B.2 (Sampling from (-DPP(A), adapted from Anari al. [2024]). Given positive-definite matrix A, there exists an algorithm that draws 2 approximate samples from 26-DPP(A) A) in time (nb log" 4 n. + 0b log" n). Furthermore, each drawn sample is indistinguishable from an exact sample from 2b-DPP( A) with probability at least 1 - -n -0(1)
We now provide the main result of this section, incorporating the above results in the analysis of beetch-and-projess (SAP) along the top-l subspace of A. Since throughout this section we focus on a linear system with single right-hand side, we will w, instead of W. denote the SAP iterates, as these are always vectors. The update rule for SAP is now:
W++1 = wi - 5T (SAS") 55 - y)
equivalently / w.) (1 11 (w: w.).
(5)
where II AI=AA' S/ and w. y. Let II denote EII.
In the rest of the section we prove the following result, which is a slight generalization of Theorem 3.2. Theorem B.3 (Fast convergence along top-l subspace). Let b < n/2 be fixed and S RR 2b2 be random subsampling matrix sampled from 2b-DPP(A) A). Furthermore, let 1 S l < be also fixed and (II) > 2An (II). For any t >> 1 define w. E 2i/1/2
1/2 wi where wi generated using (5). Then,
1/2
E/ - -w-)le. S (1- I )
1
|| (wo w.)l3
1/2
8
(1-
1
7 [0(b,)) 1)
||w0 w. 113
20(b,n)
where ((b,p) = =bE
d_ x
The assumption (II) 2 (II) not restrictive buts lets us provide a cleaner analysis by avoiding the edge case of (III = (II). In fact, in this corner case, we can simply rely on the existing SAP analysis from previous works and recover our main result, Theorem 3.2. For completeness, we derive Theorem B.9 for posterior GP mean inference in Section B.3 without the assumption (IT) 2 2A, (IT). The proof of Theorem B.3 appears after the proof of Theorem B.8. We build towards the proof starting with the following result for SAP [Dereziński and Yang, 2024]:
Lemma B.4 (Linear convergence with SAP). Let S E Rb2xx be a random subsampling matrix sampled
from from A6-DPP(A). Then,
El|we+1 111 s (1 X, (T))El w.113.
However, in the following result, we show that along the top-l eigenspace of A, the expected iterates can converge at a much faster rate than 1 - (TI).
Lemma B.5 (Convergence of expected iterates along top-l subspace). Let S E 2bxx be a random subsampling matrix sampled from 2b- Ab). Then
|| [we+1 - w. 113. < (1 >x(TT))*
(((T)))(((/((()))) ||A//( - w. m12.
where E, denotes conditional expectation given wt-
Proof. We have
Q/ - w.) Q 1))11/2 (u) -
19

---

Taking expectation on both sides and squaring we get,
E - mx1ll2, = (w: wx, As A = VDVT we have IT = V /D (SV (VV DV Due to Lemma B.1 we know that IT and A share the same eigenvectors, therefore, I1 and Qr commute. Consequently,
1 - -w.).
| -00t++1-w11112
(w, w.)² ( (1 S (1 - (1) (l,
22 -w.e
We now analyze convergence along the top-l dimensional subspace in L2-norm by considering E/| (wi+1 w.).l We have,
(w - w. lla, = (w: - w.)" w., A - - 11)] (w, - w.).
In partieular, need upper bound E [ Qz) [I]]. We prove the following lemma:
Lemma B.6 (L2-norm enror along top-l subspace). Let S R 2bxn be random subsampling matrix sampled 2b-DPP(A). Then
E/ (Wi+1 w. 1.), S (1 (e(1)) (w, w.)ll. A (1) w.)ll3
Proof. First, rewrite E [[4 11) - 11)] follows:
E - 11)=e(1 - 11)] Q2 nQQr - E[IQu
(1 - 11) - Q + B[IQ
(6)
Crucially, the first two terms terms in (6) live in the top-l subspace, but the third term does not. Nevertheless, are still able to bound it follows:
B E|I(I - (1 - Q.)) E[]] - EQ Q+) 3 π ni( Q.)
where we used Jensen's 's inequality in the last relation. Substituting in (6) we get,
[ - IT)Q,(1 IT) = IT); (T (T)
(7)
Using (7) we now upper bound //I w.
w,)lla 11..
CIIAA w. I12. s (1 (11) (.lu:w.)llle w.)lla. Ae+1 (1) w. w.)lli-rr
which concludes the proof.
Unrolling this recursive bound, and combining it with the convergence in the full norm, we obtain the following convergence guarantee for SAP iterates without tail averaging.
Corollary B.7. Let
(m)) Then, T-A,
E|| (w++1 -ao--w))l w.)113. - (1 (1 - w.). + Aet (
(1 A, (T)))(luo - -w. -.113.
1-a
Proof. Define y 1 - (II). Invoking Lemma B.6 deduce,
A (W++1 w 112. < (1 ((T) w. w + det1 (fi(I M))E[w w.
Now, recursively applying the tower rule and Lemma B.4 yields
B1]44/2(w1+)) 11), - (1 >e (TT) ( w.)1 + (TT(( (T((T-)))>) (wo 1.1.2.
Unfolding the preceding display yields,
E| 2/1A2/2(224 - w.)1lg. I (1 (1 - - 11.)
+ (I((
Σ
711 ((())))-)) ||w0 20o--wll- w=113.
1=0
20

---

Observing that Ei=07'(1 ((t)) -1- y -i= a where = 1-Ar(10 < 1, we have
BB|AA" (w++11 w.)lla, <1 Xe(i)) +1 ( - w.)ll,, + Aet1 (II(( II)
7 f11wo - w. lla
1
which concludes the proof.
We now use the tail averaging idea similar to Epperly et al. [2024a], obtaining fast convergence along the top-l subspace. The proof of Theorem B.3 is then derived from the following result, after combining it with fixed-size DPP sampling guarantees of Lemma B.I.
Theorem B.8 (Fast convergence along subspace for tail-averaged iterate). Let S E 2bxn be random subsampling matrix sampled from 2b-DPP(A). For any t > 2 define wi 2 > t-1 where Wi are generated using update rule (5). Then
2121/2 W)
21 -w))la/ S (1 Ao ((T) | (wo -w,)l
+ Aet1 (TI -I)
t(1 a)² (1 An 1 w. 11á
1-Ax(I0) where a 1- (II)
2/2 wi. Consider / - w. 12a, set la = 1/2, and let Ew,, denote the expectation conditioned on wy. Applying the tower rule yields
Proof. Let úi = Ei2//2
1
(l -w.)lle. (t. - 12)²
El(w, w. A (w,) w.)]
(w, - A. E.,, 1QrAt/2 (0. - w.)]
(t )² 2
s=1a r=ta
For r 5, we have
E QiA w.)] Qe(l E. 1/ (/,2(WW11 w.]]
= QiQe(T - IT) [A2/2( w+)
Qe -1)[- 1 (1 w.)l
Recursing on the above relation yields Ew, 1/2 / (wm
Therefore.
(-1
/||A2//(10x w-)l12 = (t - ta
E - (
w.)
s= r=ta
Σ (1-Ae ((--((1))==))1172 -w-)l1.
(1 12)2
e=1, r=fe
Applying Lemma B.6 upper bound E E//A1/2(w;
w.) 18 obtains
1-1
E|| /|IAA//((0; w.)1l., S (1 - ta)z
1
| (wco w.)1l3. Σ Σ (1- (II))
r=t.
t-1 de+1 ([[[- || ||Al/2 w-)l?
Σ -
+ (t 1 1
1-a
<1 (111))11111//2 w.))lle.
+ det1 ([[(I -mi)] || ///1/2 1000-w,)''' wx)ll [t (t- tat
1
1-a
s=ty r=
21

---

where y = 1 (II)
Using Eist, Li=t, a= < (t - (o) Ei=o a' = i-ts we get,
1-a
Aet1 I1)) Bl|A4/ w1) S (1 ((lI) "M "|A w. ll + (1- (1-1))1 a)² /' 1 w.)ll?
Substituting la = 1/2,, the result immediately follows.
Completing the proof of Theorem B.3 It remains to use the guarantees for the eigenvalues of the matrix II from Lemma B.I.
Proof of Theorem B.3. Recalling the definition of a from Theorem B.8 and our assumption that (TT) 2 22n (IT), we obtain for y = Xn (IT) that (1 - a)- < 27/ (e()) Consequently,
Xe+1 [II(I II)
A< -i))
<
(1 a]=
(1 a]²
1(1 de((T) g2
(III
(e(IT) 1
(1)
Σ
Xe
i>d
(2)
40(b, p)y,
where (1) applies Lemma B.1 and (2) defines ((b, E) = 1/2 Ei>b 2. inequalities: 2Φ(b, n)
Observing the elementary
Xe-
1 Xe(II) S 1 1 1
we immediately deduce from Theorem B.8 that
1/2
E11A4//(00; (w..- - (1 It (+0(6,))
||A w.)l3.
8 + (0b,))
(1-
1
1/2
1|wo w.
n)
B.3 Posterior mean inference along top-l subspace for GPs
We begin by providing some background on GP inference in the Hilbert space setting. This allows for graceful transition from Hilbert space norm over the posterior mean to vector norms over IR... We recall that f is Gaussian process and { w)) J represents the training data. The posterior Gaussian process characterized by N(m,, (+), Knp (s )), where
(. = m(-) kc- X) (K + A k ( k(. k(, X)(K AT)-'A(X..))
Let H be the reproducing kernel Hilbert space (RKHS) associated with the kernel k(.,)) .). Assuming m(-) = 0, the mean function (-)(-) can be identified as an element of the subspace H., defined as
Hn
w;k((,ar;) w E 1
Furthermore, note that In particular, ma Ei (wx) ).k(,, 24), where w. (K A1)- for any element m' E Ha. have |lm''122 nD Kw || The operator C.
22

---

E k(x,;;) @ k(),zz) ) is known as the empirical covariance operator. Let by denote the jth unit eigenvector of K with eigenvalue Aj. is straightforward to show that uj == VA, Li= uj,k(,, x4) unit eigenvector of the unnormalized empirical covariance operator Ei A xz) 8 k(,, x4) with eigenvalue Aj- Let V2 E R consists of top-l orthogonal eigenvectors of K as columns and Q = V be a projection matrix onto the subspace spanned by u1...., Ug- Consider the E-dimensional subspace He defined as
72 He Σ (Qew),k(;
w E R
We claim that He is the subspace formed by top-l eigenvectors of the empirical covariance operator This can be seen clearly by choosing Vj for 1 S S l, we get uj E He. We have the following conclusions.
For any m' Ei= w;k((; have I|m' mn b
|| w. l
The element m/q, = Ei= (Qew)ik(,, mx) is an orthogonal projection of m' onto He. This can be seen as
12
m' Σ 2xk( Σ k(,2,)) Σ ((--e))w, x(-
ma
and finally noting that (m/e m/2g) ( = (I((IQQ))Qaa
For any m' = ne,k(,, 2.), we have [pproj((m)) pm = 11Qc(w w. w.)llk. where proje (mr') denotes orthogonal projection of m onto H..
We now derive the main result of this section by using Theorem B.3. We replace by K2 in the statement of Theorem B.3 and use X to denote the th eigenvalue of Kx, i.c. X where A is the ith eigenvalue of K. Furthermore, Theorem 3.2 can be derived from the following result by noticing that sampling from 26-DPP(Ks) costs time ((nb² log" n) for preprocessing and an additional ( log n) for actual sampling at every iteration (see Lemma B.2). Here is the main result of the section:
Theorem B.9 (Subspace convergence for GP inference). Let b < n/2 and S E 2bxx be random subsampling matrix sampled from (b-D)P((x). For any > 2 define w; = 2 Ex=1/2 w; where wa are generated using the update rule (5). Then, sketch-and -no--oojecct initialized
satisfies
(l)proj(((mn)) pro))(mm)]13 S min 80(b,))
(1 (1- 20(b,
1/2
dr-A where fhe Er werk(, wi) and d(b,p) Eish
AptX
Proof. If (T)) 2 2A, (T), then we rely on Theorem B.3 and replace A by K2. We have the following observation: For any 2, we have
(1-1+100)11
)²=
1/2 ( o(b,8) (11(()))))
) I (6(b,
1/2
(1- 1+
1+0 11100)
=
1 (16
Therefore for all t 2, we have,
8o(b, E) (1 20(b, n²)
1/2
E E 1/2 2 w.)1l S min
||
1/2 w/2w..1..
Let mn Ei= k., where wiss denote tth coordinate of tù, and m. Em Wask(,, 24) where w, = K5 y. Then we have,
23

---

E proje((m)) - proojs(m.)lle BlIQe(wi - w.). 3
< BllQc(ui - w+)liis
min 8p(b, e)
(1-
1
1/2
| 1/2 w.ll
20)
86(b, E)
(1
1
1/2
min
2 (b, n)
On the other hand if (III 2A, (II), then we simply use the SAP analysis and get
1/2
4 E
1-1
El]Qe(wi -..)lis SB|ir
wT 113c
Σ
W
w+)
82
i=t/2
KA
Σ
Bl|ws - w. walk3x
E
Σ
E(w,, w.)Kx(w, - w.)
82
/=1/2
8,r=1/2,02r
1/2
$
| w. 112
20(b, n)
where the last inequality can be obtained using a recursive argument similar to Theorem B.8 by plugging in linear convergence guarantees for SAP (using Lemma B.4 and Lemma B.5 with l = n). We get,
1/2
E (())) - proj.(m..)12 S (1) - w.)l s (1- 20(b, n)
|| w.ll
40(6, n))
(1
1
f/2
min
|1 K
20(b, n
80(b, E)
(1- 2xh(,, n))
1/2
min
|1yll
|
K
t
Combining the results for both scenarios: (ft 2 22, (II) or (e(II < 2X, (II,, we conclude the proof.
Time complexity analysis We now use the above guarantee to provide the time complexity analysis for estimating the GP posterior mean. The following result immediately implies Corollary 3.3.
Corollary B.10. Suppose that the matrix K exhibits polynomial spectral decay, i.e., >;(K) ( O11) and E (0, 1). choosing = 28 we can find
for some 3 > Then for any E {l..... n., in that with probability at least 0.99 satisfies [|proj, (m) - -projy(m.)ll l1? I ((lproj,(mm.))) 2 in
o nL log + (ni² + nl log" n) min log(n/()) (1+ (K) +11
log(mn/te))) n(,,, ()) + +1)
time.
Proof. As y = S(X) + 9. where y (O, AT) and ((X) ~ NN K), we have y ~ NO, K2). This implies E[yv"] Kx. Therefore, EI|K 1/2 10.11 and |||Q,K] w.ll2 l. So, using standard Gaussian concentration of measure, it follows that with probability 0.999, K 1/2 w. || s OL1) || w. 1 Furthermore.
K 1.1.. (1+ 2n) |K |k w.... (1+ 2 [m))))
We get the following relation:
24

---

n11 + A/Dc)
6(6,1) (1 (1-
1
1/2
El (f)))) - pmoj.(m.)l - O(1)
min
pproj
E
_ ni)
Now let 6 > f Eiss Are + A). We have ((6,6) S t EisD N+tx +
d
n2
S 1
nd
b(x++
b(r+tA
We consider following two cases:
Case 1: B(r+t] S 1. In this case we get @(b, () < 2. After t = ( iterations we get
nd
(mpproj((m,)]13 (l)) - pm)))((m)))l S
Case 2: b(++t)) 1. We get A bλ ++ x 2.. Therefore, after 1 = O proj, (m,,) 11 S clIprojs(mm,) (m, On the other hand, either case
nd
In this have m) :/ Eisb An+
n
log nil m
iterations, we get EE] El[proje
A +
n1
(h + ()n + Aj
A de
n b
de +
(bb, m)
b
1
+
<
A, + +
$>6
1>b
as we assumed b >
Eisb
d,
Furthermore, for given spectral decay for any 3 > 1, This implies for b 26 obtain e accuracy result log(n/Ee))). where we used = O(1) to get
Aetx
have l + Eisb de
2c.
in I = o 0 min { lo)(n//)) (1 log (7(1+ (7(1+1/2e)) a/Xe))
Aet X.+
O (log(?) 7)). Then, after applying Markov's inequality, m. = m, obtains the desired e accuracy with probability 0.999. The total time complexity follows by combining the number of iterations with cost of approximate sampling from 26-DPP(A) (Lemma B.2) and additional per iteration cost of Onl + e log" n). Note that while the sampling algorithm is not exact, taking the union bound with respect the high probability guarantees in Lemma B.2 ensures that all of the samples are indistinguishable from true DPP with probability 0.999. Finally, union bounding over the three 0.999 probability events we have invoked concludes the proof.
Improved guarantee for very ill-conditioned problems Here, we show that when À is very small (i.e., Ky is highly ill-conditioned), then we can obtain an even better time complexity in the first phase of the convergence, addressing the claim in Remark 3.4.
Corollary B.11. Suppose that the matrix K exhibits polynomial spectral decay, i.e., {x() = @ for some B 1 and X < Ib then we can find in that with probability at least 0.99 satisfies [[proj (m.)) S (l(projj((m,)) 4 in
(nb² log n + (n² 2 + nb² log ()(B/))=-1/e)
Proof. We reconsider case 1 from the proof of Corollary B.10. Using the given spectral decay profile
for K and that À < NT we get o(b. f) = O (b) Therefore after t = n 6-8 approximation guarantee. The total runtime complexity becomes:
Ee t we obtain the e
(nb² log n + (n² nb² log m)(e/b)8-1/c)
C Additional Algorithmic Details for ADASAP
We provide additional implementation details for ADASAP. Appendix C.I describes how we dis- tribute matrix-matrix products across rows and columns in the algorithms ColDistMatMat and RowDistMatMat. Appendix C.2 describes the practical implementation of the randomized Nyström approximation and provides pscudocode for RandNysAppx. Appendix C.3 describes how we compute
25

---

preconditioned smoothness constants and provides pseuodocode for GetStepsize, Appendix C.4 provides pseudocode for Nesterov acceleration in NestAcc. Appendix C.5 provides scaling plot illustrating the speedups achieved by using multiple GPUs in ADASAP. Appendix C.6 investigates the impact of tail averaging on the performance of ADASAP.
All operations involving kernel matrices are performed using pykeops [Charlier et al., 2021], which allows us avoid instantiating kernel matrices explicitly in memory. To see the full details of our implementation, we recommend the reader to view our codebase.
C.1 Distributed matrix-matrix products
Here, provide details for how implement the distributed matrix- products in ADASAP. ColDistMatMat (Algorithm 3) shows how we distribute the matrix-matrix product KBn Wi in ADASAP and RowDistMatMat (Algorithm 4) shows how distribute the calculation of the Nyström sketch in ADASAP. Our implementation of these algorithms uses torch forch.multiproeessing to spawn CUDA context on each device (i.e., worker) and uses pykeops generate the column and row block oracles.
Algorithm 3 ColDistMatMat
Require:
Right-hand side matrix W E TR row indices B, workers (Wi,... W Partition rows of W as (Wi.... W,,,,WW.. WN Send W; to W;
WN Waokkrr
Generate column block oracle Xyy using B (K.,.W): Kw [W.]
Compute column block products in parallel
Aggregate KB. W < E Nuorkore (KBn W): return KBn W
Algorithm 4 RowDistMatMat
Require:
Right-hand side matrix 52 e R"xr, row indices B, workers (Wi.... WN workems
Send Ω to each W;
Generate row block oracle x.W. using
(KBBl)), t K [St]
Compute row block products in parallel
Aggregate K88Ω + [(KBBSt) (KBBl) Morlinee.
return KES
C.2 Randomized Nyström approximation
Here, we present a practical implementation of the Nyström approximation used in ADASAP (Algo- rithm 2) in RandNysAppx (Algorithm 5). The Randomized Nyström approximation of KgS with test matrix Ω E Rxxr [TTropp et al., 2017] given by:
K (KB KBS (B)l)) (KBB)))
The preceding formula is numerically unstable, so ADASAP uses RandHysAppx, which is based on Tropp et al. [2017, Algorithm 3]. eps(æ:) is defined as the positive distance between and the next largest floating point number of the same precision as The resulting Nyström approximation M given by UST where U E [gr is an orthogonal matrix that contains the approximate top-r eigenvectors of M, and S E R" contains the top-r eigenvalues of M. The Nyström approxima- tion positive-semideffnni but may have eigenvalues that are equal to 0. In our algorithms, this approximation is always used in conjunction with regularizer ensure positive definiteness.
The dominant cost in Algorithm computing the SVD of B at a cost O(pr²). This the source of the 3pr cost in Phase III of ADASAP.
26

---

Algorithm 5 RandNysAppx
Require: sketch MΩ E Roxr, sketching matrix Q E Rexy, approximation rank S p A < eps(12 MS2.dtype) Tr( MΩ)
D Compute shift for stability Cholesky decomposition: CTC S2" MΩ + AΩ Ω
C chol( MΩ + 2
B < YC-1
b Triangular solve Thin SVD
[U, <- svd(B) (,
S + max diag(S² A
D Compute eigs, and remove shift with element-wise max
return U,S
C.2.1 Applying the Nyström approximation vector
In our algorithms, often perform computations of the form (M + 9 (USU + pl)- This computation can be performed in ((rp) time using the Woodbury formula [Higham, 2002]:
(USUT + - (S UT g 9- UUT
(8)
We also use the randomized Nyström approximation to compute preconditioned smoothness constants in GetStepsize (Algorithm 6). This computation requires the calculation (P + - for some w E IR,, which can also he performed in (pr)) time using the Woodbury formula:
(USUT + p11-- = u + pl) 1/2 U U - (v- (U-UUT) v
(9)
In single precision, (8) is unreliable for computing (P pf)- This instability arises due to roundoff error: the derivation of (8) assumes that Ü = I, but we have empirically observed that orthogonality does not hold in single precision. To improve stability, we compute Cholesky decomposition LL, of ps + U, which takes O(pr²) time since we form U U. Using the Woodbury formula and Cholesky factors,
1
(USUT + p
U + +0-
p
P
1 1 UL. UT.
9
p
This computation can be performed in Opr)) time, since the O(r²) cost of triangular solves with LT and L is negligible compared to the Oppr) cost multiplication with U and U.
Unlike Eq. (8), we find using (9) in GetStepsize yields excellent performance, i.e.. we do not need to perform any additional stabilization.
C.3 Computing the stepsize
Here, we provide the details of the GetStepsize procedure in ADASAP (Algorithm 2) for automat- ically computing the stepsize. Our procedure is inspired by Rathore al. [2025], who show that approximate SAP with the Nyström approximation converges when
1/A1 (PPs 1/) AI)(PB pl)-1/2)
TIS
That is, the correct stepsize use is the reciprocal of the "preconditioned subspace smoothness constant"..
To compute \1 ((Ps + 1l) + (PB + 1/2), GetStepsize uses randomized pow- ering [Kuczyński and Wožniakowski. 1992]. This technique has been used in several previous works on preconditioned optimization to great effect [Frangella al., 2024a,b, Rathore et al., 2025]. Given a symmetric matrix H, preconditioner P, and damping P. randomized powering computes
( H (P pT) 11)-1
using only matrix-vector products with the matrices H and p11-/) When P is calculated using RandNysAppx, GetStepsize can efficiently compute matrix-vector product with (P + -1/
27

---

taxi (n 3.31 10P)
4.0
3.5
Speedap
2.3
2.0
1.3
1.0
Namber GPUs
Idval Spredhap
Messured Spredup
Figure 4: Multi-GPU scaling of ADASAP the taxi dataset. ADASAP obtains near-linear scaling with the number of GPUs.
using (9). In practice, we find that 10 iterations of randomized powering are sufficient for estimating the preconditioned smoothness constant. The pseudocode for randomized powering, based the presentation in Martinsson and Tropp [2020], is shown in Algorithm 6. Since ADASAP only runs GetStepsize for 10 iterations. the total cost of the procedure is O(b?), which makes the source of the (b²) cost of Phase III in Algorithm 2.
Algorithm 6 GetStepsize
Require: symmetric matrix H, preconditioner P, damping maximum iterations N + 10 0 + randn( (
to + ub/ /o/l/oo|22
> Normalize
for i = 0, I..... N do
U1+1 + (P+ p
U/41 +- Hu1+1
U:11 - (P+ pP+ -1/2
Ui41 + 0241/1lus411/
> Normalize
end for
X < (UN- (UN-1)=oo) JNNN
return 1/A
C.4 Nesterov acceleration
We present pseudocode for Nesterov acceleration in NestAcc (Algorithm 7).
Algorithm 7 NestAcc
Require: iterates W,, V,, Z4, search direction D,, stepsize f/is, acceleration parameters 8, a
W Z, ng D,
V-11 < 3V, + (1 B) 74
miB Dr
741 <- aV,
a) +1
return W441, Ve+1, Ze+1
C.5 Parallel scaling of ADASAP
Here present Fig. 4, which shows the parallel scaling of ADASAP on the taxi dataset.
28

---

b
100
109
109
10-1
16-1
Sabaparre Groot
Sulequee léror
Salequce Ecror
10m2
10
10
10-1
10
10-1
10
10-1
10-6
10-²
10
Data Pasees
Data Pasere
Data Paeses
100
10-1
10-1
10-1.
Salepooo Etree
10-
Satepace Elnee
Salepero Enue
10-
10
10"
10-4
10m3.
10-3
10-5.
10-7.
10
70
10-7
10
0
10
0
10
Data Pasees
Data Passes
Data Pieses
(noo-nvemagd))
(non-avemged))
(pom-aoeraapd)
= 10 (tail averaged)
190 (tail averaged)
1000 (tail averaged)
Figure 5: Performance of ADASAP with and without tail averaging. One "data pass" corresponds to pass through the kernel matrix. Tail averaging does not improve convergence substantial margin.
C.6 Impact tail averaging on performance
The theoretical results in Section 3 require tail averaging, Here we assess whether tail averaging leads to practical improvements in the convergence of sketch-and aesch-and-project algorithms for kernel ridge regression. To do so, we run synthetic experiment using an RBF kernel with 1000 samples.
Fig. 5 displays the relative errors the top-E subspace, which computed using the expression
|proje((it) (mo..(m...) 112// ///oo/ (m..))..
When / = 10 and the blocksize is small, tail averaging slightly improves the convergence rate. However, when E {100, 1000}, tail averaging does not improve the convergence rate. In fact, as the blocksize increases, tail averaging results in slower convergence!
D Additional Details for Experiments
Here we provide additional details for the experiments that are not provided in the main paper.
D.I Determining hyperparameters for regression
We use zero-mean prior for all datasets. We train the kernel variance, likelihood variance, and lengthscale (we use separate lengthscale for each dimension X) using the procedure Lin et [2023], which we restate for completeness:
29

---

1. Select centroid point from the training data X uniformly at random.
2. Select the 10,000 points in the training data that are closest to the centroid in Euclidean norm.
3. Find hyperparameters by maximizing the exact GP likelihood over this subset of training points.
4. Repeat the previous three steps for 10 centroids and average the resulting hyperparameters.
D.2 Optimizer hyperparameters
We present the hyperparameters for ADASAP, ADASAP-I, PCG, and SDD that were not described in the main paper.
For GP inference on large-scale datasets, we use blocksize b = nz/100 ADASAP, ADASAP-I, and SDD; blocksize b n n220000 for transporation data analysis, and b = n/5 for Bayesian optimization,
We set the rank 100 for both ADASAP and PCG.
Similar to Lin et al. [2024], we set the stepsize SDD to be one of (1/n, 10/n, 100/n) (this grid corresponds to SDD-1, SDD-10, and 100, the momentum 0.9, and the averaging parameter to 100/Tmax-
D.3 Additional details for GP inference experiments
song and houseelec are from the UCI repository, yolanda and acsincome from OpenML, and benzene and malonaledehyde are from sGDML [Chmicla et al., 2017]. We select the kernel function k for each dataset based on previous work [Lin et al., 2023, Epperly et al., 2024b, Rathore et al., 2025].
We standardize both the features and targets for each dataset. For fairness, we run all methods for an equal amount of passes through each dataset: we use 50 passes for yolanda, song, benzene, and malonaldehyde and 20 passes for acsincome and houseelec.
We use pathwise conditioning with 2,048 random features (approximately) sample from the GP posterior.
D.4 Additional details for transporation data analysis
We standardize both the features and targets and RBF kernel. Due computational constraints, we use a single 99%-1% train-test split, and run each method for a single pass through the dataset.
D.5 Additional details for Bayesian optimization
Our implementation of Bayesian optimization largely mirrors that of Lin et al. [2023]. We only present the high-level details here, and refer the reader to Lin et al. [2023] for the fine details of the implementation.
We draw the target functions f [0, 1] - R from a zero-mean GP prior with Matérn-3// 3/2 kernel using 5,000 random features. At each iteration, we choose the acquisition points using parallel Thompson sampling [Hernández-Lobato al., 2017]. As part of this process, we use the multi-start gradient optimization maximization strategy given in Lin et al. [2023]. At each step, we acquire 1,000 new points, which we use to evaluate the objective function. Concretely, if Knew is an acquired point, compute Vuew f(znew) training data for the next step of optimization. In our experiments, we initialize all methods with a dataset consisting of 250,000 observations sampled uniformly at random from [0, 1]s.
where (00 10-5). We then add (ænew Vnew) to the
D.6 Additional timing plots for Section 5.1
Here we present timing plots for the datasets used in Section 5.1 that were not shown in the main paper. Fig. 6 shows that all the methods (except PCG) perform similarly on both test RMSE and test mean NL... However, Fig. 7 shows that ADASAP and PCG attain a much lower train RMSE than the competitors. This suggests that the similar performance of the methods test RMSE and test
30

---

vilverda (a 340. 109)
-4 1.0-1000 109
1.5 (F)
215
1.00
1.75
0.23
0.50
110
vilordda (n 3.00 10P)
waap (n 68 169)
aroineeme (n 1.50
24
2011
3.70
1.2
Tho
50033
Figure Performance of ADASAP and competitors on RMSE and mean NLL, as a function of time, for benzene, malonaldehyde, and houseelec. The solid curve indicates mean performance over random splits of the data; the shaded regions indicate the range between the worst and best performance over random splits of the data. ADASAP performs similar to the competition.
yolanda (n 3.60 109)
song (n 4.64 109)
2.0
2.0
1.5
1.5
Toan RMSE
RMSE
1.0
1.0
Than
0.5
0.5
0.0
0.0
600
Tame
Time
PCG
SDD-10
ADASAP
SDD-1
ADASAP-I
Figure 7: Performance of ADASAP and competitors on train RMSE, as a function of time, for yolanda and song. The solid curve indicates performance random splits of the data; the shaded regions indicate the range between the worst and best performance over random splits of the data.
mean NLL is not due to differences in optimization, but rather, because the datasets are not well-modeled by Gaussian processes.
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The theoretical and empirical claims made in the abstract and introduction are justified in the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: We discuss the limitations of our work. In particular, we describe the key limitation of our work, which is that we do not show a strong theoretical guarantee for ADASAP due its differences with SAP.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
32

---

Answer: [Yes]
Justification: We clearly state all assumptions in our theorem statements, and a proof of each claim is provided within the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Yes, we include all the code to reproduce our experiments, and provide the details for the experiments within the paper itself.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
33

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We have provided all the code for the experiments GitHub, and linked to the code in the main paper. We also provide instructions in our GitHub repository on how to download the data.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.cm co/
guuiees/Coeeubbln
for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https:
/ /nips. wwtppcccccc/ccccccccccccccccctct gwwpc//gg//////iiioiiooicc /tt/ggddiC///oeuSmmiiiiii
(CodeSubmissi wwwiisssssioiooooooooo
for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: We provide all of the essential details, such as data splits, hyperparameters, and optimizers within the main paper and appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Our experiments for GP inference and Bayesian optimization five random splits to ensure that the improvements attained by ADASAP are not just due to randomness. We also show the worst and best performance of each method in the timing plots, which demonstrates that ADASAP is more robust to the random seed compared to the competitors.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
34

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπτог rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We clearly state the hardware and number of iterations that we use to run each algorithm.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https:/ /ttpr///////u..eeeeeee cctppblic//etteeeeceeeeee petublcc/thiccccideeineeee
Answer: [Yes]
Justification: After reviewing the NeurIPS Code of Ethics, we believe that our research conforms with this code.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., - if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
The work performed in this paper does not have clear societal impact.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
35

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), faimess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We do not believe the datasets used in this paper require safeguards for use.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but encourage authors to take this into account and make best. faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Yes, properly credit the authors of data and code that use in experiments.
Justification:
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
36

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. www.coo//aaooo/o datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: This paper involves neither crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: This paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
37

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
We do use LLMs for any of the core methods in this research.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: //twerris...../eeee cc/Conferences//2222 2002LLL for what should or should be described.
38


Reviews:
[{"review_id": "RK8qZ8wJ4t", "reviewer": "Reviewer_P6N1", "summary": "The authors propose a new algorithm for Gaussian Process (GP) inference, with a focus on addressing large-scale problems. They identify two key issues in previous methods (PCG and SDD): robustness to ill-conditioning and scalability. Their method addresses both, offering good default hyperparameters and distributed execution capability. The paper also provides convergence guarantees for the posterior approximation under ill-conditioning (applicable more generally to sketch-and-project methods), and demonstrates superior performance on a series of benchmark tasks.", "questions": "1) When the dataset becomes very large (as in the main scenario considered in the paper), in my experience, it often makes sense to move away from GPs and use alternative models like Bayesian Neural Networks or Deep Kernel Methods, which don’t suffer from the same scaling issues. Could the authors elaborate on why GPs remain a preferable choice in this context? Some comparative discussion would strengthen the case for GPs.\n\n2) Existing state-of-the-art GP methods already handle datasets on the order of 10^6. So, could authors provide any details on how much additional value is gained by pushing beyond this scale? This ties into the previous question about the motivation for continuing with GPs at such scales.\n\n3) I didn’t fully grasp two points that might be critical. Could the authors clarify what they mean by “condition number-free” and what exactly is considered “ill-conditioning”? I assume this refers to challenges in inverting the matrix K + \\lambda I, but I’m not entirely certain that’s the intended meaning.\n\n4) In Algorithm 1, the variable ‘b’ appears without an immediate explanation. While it becomes clearer later in Section 3.2, it would help to define it properly when it first appears, especially since it plays a role in justifying the SAP algorithm’s reduced cost. Also, practical guidance on how to choose the block size ‘b’ would be very helpful.\n\n5) Could the authors provide references or direct links to the proofs of Theorem 3.2 and Corollary 3.3 in the Appendix when they are first introduced? Currently, it’s a bit inconvenient to search the paper and then the Appendix for these proofs.", "limitations": "The authors acknowledge certain limitations, such as not covering approximate preconditioning, acceleration, or uniform sampling in their theoretical analysis, which is appreciated.\n\nHowever, the discussion of limitations and potential future work could be expanded. For example, what kinds of problems are particularly challenging for the method? How do the authors recommend practitioners use the method in different contexts? It would be helpful to draw comparisons to non-GP alternatives.\n\nWhile the authors emphasize that the algorithm performs well with default hyperparameters (which is great for adoption), more discussion around what might go wrong and how one should adjust the hyperparameters accordingly would be useful for practitioners.\n\nLastly, should the term “large-scale” be explicitly mentioned in the paper’s title? It could help clarify the paper’s focus and improve visibility beyond the GP research community.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 2, "originality": 4, "strengths_and_weaknesses": "Strengths:\n\n1) The proposed algorithm performs well, and its ability to scale to datasets with more than 3 x 10^8 observations is particularly impressive.\n\n2) The paper establishes convergence of the posterior mean, and the theoretical contributions appear sound (though I have some follow-up questions in the respective section).\n\n3) The distributed nature of the algorithm is a practical strength that could boost its adoption.\n\n4) The paper is well-written and, to the best of my knowledge, presents original work.\n\nWeaknesses:\n\n1) The justification for the algorithm primarily focuses on the posterior mean (which is indeed crucial), but GPs are often valued for their uncertainty estimates. The quality of uncertainty quantification isn't clearly demonstrated here.\n\n2) There are some motivational gaps, which I address in the “Questions” section. I assigned a significance score of 2 because the problem being addressed feels quite specific to the GP community. I'm open to being convinced otherwise.\n\n3) Some aspects of the solution are not explained in sufficient detail (e.g., the notion of \"condition number-free\"), which affected my clarity and confidence scores.", "was_revised": true, "final_justification": "The paper presents a well-motivated and technically solid contribution to scalable GP inference, with strong empirical results and theoretical backing. While further evaluation of the uncertainty estimates would strengthen the work, the proposed method addresses important limitations of prior approaches and represents a meaningful advancement.", "rebuttal": "> **Comment 1.1** — The justification for the algorithm primarily focuses on the posterior mean\n(which is indeed crucial), but GPs are often valued for their uncertainty estimates. The quality of\nuncertainty quantification isn’t clearly demonstrated here.\n\nOur method can be used for both mean and uncertainty estimates, as demonstrated in our\nBayesian optimization experiments, which require both. Our exposition focuses on the mean just for\nclarity. We'll clarify in the paper that our method can handle both the mean and uncertainty estimates.\n\n> **Comment 1.2** — When the dataset becomes very large (as in the main scenario considered in\nthe paper), in my experience, it often makes sense to move away from GPs and use alternative\nmodels like Bayesian Neural Networks or Deep Kernel Methods, which don’t suffer from the same\nscaling issues. Could the authors elaborate on why GPs remain a preferable choice in this context?\nSome comparative discussion would strengthen the case for GPs.\n\nWhile Bayesian neural networks (BNNs) and deep kernel methods (DKMs) are attractive\nalternatives for Bayesian modeling, we believe that GPs still retain advantages in certain scenarios.\n\nFirst, regarding scalability, current state-of-the-art implementations of BNNs and DKMs rarely scale\nto large datasets (e.g., $n \\gg 10^6$ samples) (Ober et al., 2021; Arbel et al., 2023). In contrast, our work\ndemonstrates full GP inference on datasets with $n \\sim 10^8$ samples, which is two orders of magnitude\nlarger than what is typically feasible for BNNs or DKMs.\n\nSecond, in terms of prior knowledge and interpretability, GPs allow practitioners to specify kernels\ninformed by domain expertise. If a suitable kernel is already known, using GPs would be advantageous\nover learning the kernel via DKMs. Furthermore, DKMs are prone to overfitting: Ober et al. (2021)\nshows that DKMs can overfit even more than neural networks trained using the standard maximum\nlikelihood approach.\n\n> **Comment 1.3** — Existing state-of-the-art GP methods already handle datasets on the order\nof $10^6$. So, could authors provide any details on how much additional value is gained by pushing\nbeyond this scale? This ties into the previous question about the motivation for continuing with\nGPs at such scales.\n\nWe live in an era of big data, where datasets often have size  $n > 10^{6}$.\nFor datasets that are this large, methods like PCG are unsuitable for GP inference (e.g., taxi, for which $n \\geq 10^{8}$).\nMethods like SDD can scale to datasets of this size, but they have to be carefully tuned to achieve optimal performance, as Figure 3 shows. \nHowever, such tuning is prohibitively expensive (see the runtimes in Figure 3).\nMoreover, even for datasets where $n \\sim 10^{6}$, such as houseelec in Figure 1, existing SOTA solvers cannot match the performance of ADASAP.\nFor PCG, this arises because of the time required to complete each iteration, while SDD performs poorly (even when tuned) due to ill-conditioning of $K + \\lambda I$. \nThus, in addition to being essential for tackling problems of larger scale, ADASAP solves problems at previous limit scales much more efficiently than existing SOTA solvers when $K + \\lambda I$ is ill-conditioned.\nThis is crucial as kernel matrices are often ill-conditioned in practice.\n\n> **Comment 1.4** — I didn’t fully grasp two points that might be critical. Could the authors clarify what they mean by “condition number-free” and what exactly is considered “ill-conditioning”? I assume this refers to challenges in inverting the matrix $K + \\lambda I$, but I’m not entirely certain that’s the intended meaning.\n\nWe are happy to clarify our terminology.\nThe reviewer is correct that ill-conditioning refers to the conditioning of the kernel matrix $K + \\lambda I$. \n$K + \\lambda I$ is ill-conditioned when the condition number $(\\lambda_1(K)+\\lambda)/(\\lambda_n(K)+\\lambda)$ is large.\nAs discussed in Section 2.2 (Line 84), the convergence rate of SDD, a SOTA exact inference method, is controlled by the condition number of $K + \\lambda I$.\nThus, SDD converges slowly on problems where $K +\\lambda I$ is ill-conditioned, as can be seen in our experiments. \nCorollary 3.3 (Line 169) shows that for certain levels of precision, the convergence rate of the posterior mean along the top $\\ell$-basis functions is completely independent of the condition number of $K + \\lambda I$ and the value of $\\lambda$. \nSee Remark 3.4 (Line 172) for details. \nThis is what we mean when we say ``condition number-free'' in the paper. \nWe will make the meaning of this terminology more explicit in the revision.\n\n> **Comment 1.5** — In Algorithm 1, the variable ‘b’ appears without an immediate explanation. While it becomes clearer later in Section 3.2, it would help to define it properly when it first appears, especially since it plays a role in justifying the SAP algorithm’s reduced cost. Also, practical guidance on how to choose the block size ‘b’ would be very helpful.\n\nThanks for pointing this out.\nWe agree that the blocksize $b$ should be defined earlier in the paper---we'll fix this in the revision.\n\n> **Comment 1.6** — Could the authors provide references or direct links to the proofs of Theorem 3.2 and Corollary 3.3 in the Appendix when they are first introduced? Currently, it’s a bit inconvenient to search the paper and then the Appendix for these proofs.\n\nWe will add the references to the proofs of Theorem 3.2 and Corollary 3.3 to the main paper.\n\n> **Comment 1.7** — However, the discussion of limitations and potential future work could be expanded. For example, what kinds of problems are particularly challenging for the method? How do the authors recommend practitioners use the method in different contexts? It would be helpful to draw comparisons to non-GP alternatives.\n\nThe most challenging type of problem for ADASAP is one in which the posterior mean is uniformly spread across the spectral basis functions of the RKHS $\\mathcal H$.\nIn this case, learning the posterior mean well along the top spectral basis functions is insufficient for achieving good predictive performance.\nWe note this challenge is not unique to ADASAP, since methods like PCG and approximate kernel methods also suffer in this regime.\nFortunately, this setting is atypical for kernel and GP problems. \nAs we mention in the beginning of Section 3, many prior works (Dicker et al., 2017; Lin et al., 2020;\n2023) have shown that the top $\\ell$-spectral basis functions (for some $\\ell \\ll n$) are most responsible for achieving good generalization both in theory and practice. \nNevertheless, we will add discussion of this problem class to the revision.\n\nWe believe ADASAP will be valuable to practitioners whenever they are using GP and Bayesian methods that require the solution of large kernel linear systems.\nWe would recommend replacing solvers like Choleksy and PCG (which are unfriendly for large-scale problems) with ADASAP. \nNote that linear systems of this size arise in models that go beyond the basic GP setting.\nFor instance, we believe ADASAP could be valuable for scaling deep kernel earning to larger datasets, as certain derivative calculations in the log-likelihood require solving a linear system with a large kernel matrix (see Eq. 7 in Wilson et al. (2016)).\nWe will include this discussion in our updated revision.\n\n> **Comment 1.8** — While the authors emphasize that the algorithm performs well with default hyperparameters (which is great for adoption), more discussion around what might go wrong and how one should adjust the hyperparameters accordingly would be useful for practitioners.\n\nWe agree that the submission could benefit from including more details about how to tune the hyperparameters for ADASAP. Here are some examples, for intuition:\n* Blocksize $b$: Increasing $b$ should lead to faster per-iteration convergence, at the cost of more expensive iterations (both in terms of wall-clock time and memory). \nOn the other hand, decreasing $b$ should lead to slower per-iteration convergence, while leading to cheaper iterations.\nIn general, we recommend selecting $b$ to maximize GPU utilization without running out of memory.\n\n* Rank $r$: The behavior for adjusting $r$ is similar to that of adjusting $b$.\nLarger values of $r$ lead to faster convergence, but more expensive iterations; smaller values of $r$ lead to slower convergence, but cheaper iterations.\n\n* Acceleration parameters $\\mu$ and $\\nu$: Increasing $\\mu$ (or decreasing $\\nu$) could improve the convergence rate in practice.\nHowever, practitioners should be careful to ensure that $\\mu < \\nu$, otherwise this will make some of the Nesterov acceleration parameters negative, which could lead to unstable behavior/divergence.\n\n**References**\n\nJulyan Arbel, Konstantinos Pitas, Mariia Vladimirova, and Vincent Fortuin. A primer on bayesian\nneural networks: Review and debates. arXiv preprint arXiv:2309.16314, 2023.\n\nLee H Dicker, Dean P Foster, and Daniel Hsu. Kernel ridge vs. principal component regression:\nMinimax bounds and the qualification of regularization operators. Electronic Journal of Statistics,\n11:1022–1047, 2017.\n\nJihao Andreas Lin, Javier Antoran, Shreyas Padhy, David Janz, José Miguel Hernández-Lobato,\nand Alexander Terenin. Sampling from Gaussian process posteriors using stochastic gradient\ndescent. In Advances in Neural Information Processing Systems, 2023.\n\nJunhong Lin, Alessandro Rudi, Lorenzo Rosasco, and Volkan Cevher. Optimal rates for spectral algorithms with least-squares regression over Hilbert spaces. Applied and Computational Harmonic\nAnalysis, 48(3):868–890, 2020.\n\nSebastian W Ober, Carl E Rasmussen, and Mark van der Wilk. The promises and pitfalls of deep\nkernel learning. In Uncertainty in Artificial Intelligence, pages 1206–1216, 2021."}, {"review_id": "EPnhNtfhXP", "reviewer": "Reviewer_UCmX", "summary": "The paper provides a new optimization method for Gaussian Processes that is very scalable.  It builds on the Sketch-and-Project (SAP) framework that uses a Nystrom sampling sketch at each step of the optimization.  The advantage of the new approach is better dependence on batch size b, by using a numerically stable method to update the estimation of the sampled gram matrix K.  The procedure also uses Nesterov's acceleration and can be parallelized.  \n\nThe paper provides some analysis which shows the convergence depends on the stable rank of the gram matrix, and with sub-quadratic update time if K exhibits polynomial spectral decay.  \nThe experiments use cross-validation on large datasets to validate the method.  They show comparison to a couple leading methods from the past couple years, and show better fit.  They also scale to a NYC Taxi dataset with over 30 million data points -- where other methods fail.", "questions": "I would like to see more ablation studies to explain which of the improvements on any baselines and the basic SAP make the biggest contribution in the scaling.  Is this possible to do?  If not, please try to give insight into which components are most essential in the improvements.", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 2, "significance": 4, "originality": 3, "strengths_and_weaknesses": "There is room for improvement in the writing.  The algorithm description came after the section describing the convergence results.  So it was fairly difficult to understand the paper without going back and forth several times.  \n\nAlso, there seems to be a few extensions to the basic SAP framework.  The faster update, the acceleration, and the parallelization.  While it was useful to see the ADASAP-I that did not do preconditioning, it would be useful to see the effect of each of those three steps (the A, D, and A) in separate runs of the algorithm so we can see which is responsible for the increased improvement.  \n\nGaussian Processes are a core technique in ML, and scaling them has been a challenge.  This paper while not centered around a single grand idea, does provide a very scalable solution for GP, and I think that make it a strong contribution.", "was_revised": true, "final_justification": "This is a real solid contribution in that it provides what is now the state of the art way to scale up a widely used approach: Gaussian Processes.  I think the experimental section could be better, but the paper explain how it works well, and the results are clear important advancements, so I think it should be accepted.", "rebuttal": "> **Comment 2.1** — There is room for improvement in the writing. The algorithm description came after the section describing the convergence results. So it was fairly difficult to understand the paper without going back and forth several times.\n\nWe appreciate hearing about your experience as a reader, and will add more guideposts and references in the revision to enable smoother transitions between sections.\nWe structured the paper so that results on each algorithm (SAP, ADASAP) appear in the section introducing that algorithm.\nThe attractions and limitations of SAP are what motivate the development of ADASAP, which we introduce in Section 4.\n\n> **Comment 2.2** — Also, there seems to be a few extensions to the basic SAP framework. The faster update, the acceleration, and the parallelization. While it was useful to see the ADASAP-I that did not do preconditioning, it would be useful to see the effect of each of those three steps (the A, D, and A) in separate runs of the algorithm so we can see which is responsible for the increased improvement.\n\nHere is how we view the impacts of the extensions to the basic SAP framework:\n\n* Faster update: Incorporating the Nyström approximation allows ADASAP to scale to larger blocksizes.\nWe did not run ADASAP without the Nyström approximation, since it would not scale to our largest datasets.\nFor example, on taxi $(n = 3.31 \\cdot 10^8)$, we run ADASAP with blocksize $b = 1.65 \\cdot 10^5$.\nComputing the Cholesky decomposition to solve the $b \\times b$ linear system with $K_{\\mathcal B \\mathcal B} + \\lambda I$ would require $\\mathcal O(b^3)$ compute and $\\mathcal O(b^2)$ memory, which is not practical for this blocksize. In terms of wall-clock time, computing the Cholesky decomposition in this setting would probably correspond to over a day on an A100 GPU.\n        \n* Acceleration: Several works have shown that incorporating acceleration into KRR solution methods leads to consistent practical improvements  (Tu et al., 2017; Rathore et al., 2025).\nTherefore, in the interest of time, we did not run experiments without acceleration.\n        \n* Parallelization: We show the impact of parallelization on the taxi dataset in Figure 4 of the Appendix. \nADASAP achieves a speedup of 3.4$\\times$ when using 4 GPUs, demonstrating that our method benefits from parallelization.\n\n> **Comment 2.3** — I would like to see more ablation studies to explain which of the improvements on any baselines and the basic SAP make the biggest contribution in the scaling. Is this possible to do? If not, please try to give insight into which components are most essential in the improvements.\n\nPlease see the response to comment 2.2.\nWe believe that the biggest contribution in scaling is due to the faster update (which allows for larger blocksizes), followed by acceleration (helps for ill-conditioned problems like KRR), followed by parallelization.\nFor smaller problems, we would not expect parallelization to be as helpful, since small problems may not fully utilize all of the GPUs.\n\n**References**\n\nPratik Rathore, Zachary Frangella, Jiaming Yang, Michal Dereziński, and Madeleine Udell. Have\nASkotch: A neat solution for large-scale kernel ridge regression. arXiv preprint 2407.10070, 2025.\n\nStephen Tu, Shivaram Venkataraman, Ashia C. Wilson, Alex Gittens, Michael I. Jordan, and\nBenjamin Recht. Breaking locality accelerates block Gauss-Seidel. In Proceedings of the 34th\nInternational Conference on Machine Learning, 2017."}, {"review_id": "fNTpEO4CBq", "reviewer": "Reviewer_WCDC", "summary": "The paper proposes a GP inference algorithm based on sketch and project (a sort of sketched Newton method) which is then combined with approximate preconditioning, distributed operations, and Nesterov acceleration. It derives guarantees (specifically, for posterior mean estimation) for the sketch and project algorithm, characterizing the convergence rate (Theorem 3.2) and time complexity (Corollary 3.3) as a function of the spectrum of the GP kernel matrix. These results show transitions from sublinear to linear rates and show dependence (or lack thereof) on a reduced (smoothed) condition number.", "questions": "See weaknesses", "limitations": "Yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths\n\nThe proposed algorithm lies between (preconditioned) coordinate descent and stochastic gradient descent methods, trading off stability and the handling of ill-conditioning. The theoretical results give an informative characterization of the convergence and time complexity of SAP for the mean path of GP posteriors. Results parallel those in the setting of kernel ridge regression developed by Rathore et al. The theory is supported by numerical evaluations comparing two exact solver baselines (PCG and SDD) on 3 tasks. The paper is accompanied by a well-structured and clean codebase.\n\n\n## Weaknesses :\n\n- The related work in Section A.2 does not substantially compare the results in Theorem 3.2 and Corollary 3.3 to existing ones. In particular, the result in Theorem 3.2 appears to closely parallel those for ASkotch from Rathore et al. but with a smoothed conditioned number. Considering the close relation between KRR and the GP mean, a proper comparison is warranted.\n\n- The \"condition number-free\" results do not appear to be particularly relevant, given that it only applies for large enough $\\epsilon$. In other words, the algorithm is conditioned number-free for large error compared to the kernel matrix conditioning. These do not seem to be situations that are meaningful in practice: if $K$ is severely ill-conditioned, we could really only expect to recover $m$ along very few dimension $\\ell \\ll n$ using a larger $\\lambda$ to regularize the problem. Otherwise, the second term will end up dominating the first one for any reasonable choice of $\\epsilon$. Hence, the \"condition number-free\" regime does not appear to be of practical interest (except perhaps in a very narrow setting). Claims of \"condition number-free\" convergence throughout the paper should be put in significantly more context or be toned down (particularly in the abstract and introduction, e.g., lines 10 and 49).\n\n- On line 39, the manuscript states that ADASAP \"obtains the robustness to ill-conditioning of PCG\". I imagine that the paper alludes here to the dependence on the smoothed condition number $\\phi$, which remains a condition number. What is more, even in the benign setting of Corrolary 3.3, any robustness gains compared to PCG appear to be due to the restriction of the estimation to the principal $\\ell$ components of the kernel matrix. As such, there is an additional \"approximation error\" to take into account, namely, $\\|m_n - proj_\\ell(m_n)\\|$ which is on the order of the tail of the spectrum of $K$. As such, the robustness gains also comes at the cost of a larger approximation error. This point (which should be made significantly clearer in the manuscript) suggest the importance of including an experimental investigation to further support the robustness claims (or additional theoretical results showing that equilibrating approximation and estimation error is not harmful).\n\n- While the paper is well-written and clear, addressing a few minor issues could improve the presentation:\n   - Some notation is not defined: $n_{rhs}$ (line 110), $K_{\\lambda}$ (Theorem 3.2), $\\rho$ (Line 199), $\\mathcal{Y}_{\\mathcal{B}}$ (Line 106)\n   - The suffix for SDD-10, SDD-100 etc. is only explained in the appendix (Line 927). This should have been explained before or in the caption of Figure 1.\n   - The significance of tail averaging is only discussed in Sec 4.2.1.", "was_revised": true, "final_justification": "I do not dissent from the consensus of accepting the paper, but I do feel that there are substantial changes in the claims that are required for the paper to fairly present its contributions and place them in the existing literature. I have kept my score at 4 because this is the stage at which I believe the current manuscript is. But if this entails rejecting the paper, then I'm happy to raise my score to 5.", "rebuttal": "> **Comment 3.1** — The related work in Section A.2 does not substantially compare the results in Theorem 3.2 and Corollary 3.3 to existing ones. In particular, the result in Theorem 3.2 appears to closely parallel those for ASkotch from Rathore et al. but with a smoothed conditioned number. Considering the close relation between KRR and the GP mean, a proper comparison is warranted.\n\nWe are happy to discuss how the theoretical results we have developed compare to those for ASkotch in Rathore et al.\nOur analysis in Section 3 focuses on the convergence of the posterior mean along the top $\\ell$-subspace using exact SAP with tail-averaging.\nOur motivation for considering this is prior works such as \n Dicker et al. (2017); Lin et al. (2020; 2023), which show theoretically and empirically that a low-dimensional top eigenspace of $K + \\lambda I$ is most relevant for achieving good generalization. \n\nIn contrast, Rathore et al. focuses on convergence of approximate sketch-and-project (with and without acceleration) *over the entire space*, and as a result, they are only able to show a fast convergence guarantee for ASkotch under certain strong conditions on the spectral decay (effectively limiting the condition number of the problem).\nConsequently, the global convergence analysis of Rathore et al. does not explain the rapid initial progress that SAP-style algorithms make on test error.\nWe believe our two-phase analysis consisting of (i) fast condition number-free sublinear convergence, followed by (ii) a slower global convergence rate dependent upon the subspace condition number, better captures this phenomenon. \nIn particular, our result is the first step in developing a systematic analysis of the generalization properties of SAP-style algorithms.\nIn the future, we would like to provide a systematic analysis of the convergence and generalization properties of approximate SAP methods by combining our approach in this work with that of Rathore at al.\n\nIn the revision, we shall add the preceding discussion on how our results compare to those of Rathore et al. to the related work section.\n\n> **Comment 3.2** — The \"condition number-free\" results do not appear to be particularly relevant, given that it only applies for large enough $\\epsilon$. In other words, the algorithm is conditioned number-free for large error compared to the kernel matrix conditioning. These do not seem to be situations that are meaningful in practice: if $K$ is severely ill-conditioned, we could really only expect to recover along very few dimension $\\ell \\ll n$ using a larger $\\lambda$ to regularize the problem. Otherwise, the second term will end up dominating the first one for any reasonable choice of $\\epsilon$. Hence, the \"condition number-free\" regime does not appear to be of practical interest (except perhaps in a very narrow setting). Claims of \"condition number-free\" convergence throughout the paper should be put in significantly more context or be toned down (particularly in the abstract and introduction, e.g., lines 10 and 49).\n\nWe agree that the condition number-free sublinear convergence in Theorem 3.2 and Corollary 3.3 only applies up to a certain moderate precision (large enough $\\epsilon$), whereas for high precision (smaller $\\epsilon$), the condition number-dependent linear rate takes over. \nHowever, arguably, the moderate precision regime (large enough $\\epsilon$) is particularly relevant in the context of GP as well as other machine learning settings, where the ultimate goal is to optimize the test error, as opposed to the training error. \nOnce the algorithm reaches a certain moderate precision in the training error, it has typically achieved close to optimal test error and further optimization yields little further benefit. \nIn fact, the phase-one sublinear convergence rate appears to align more closely with the behavior of the test error in our empirical results than the second-phase linear rate does. \n\nThat being said, we accept that our condition number-free sublinear rate is relevant only for a range of $\\epsilon$ values. \nWe will provide additional context to make this clear in the paper, particularly in the introduction.\n\n> **Comment 3.3** — On line 39, the manuscript states that ADASAP \"obtains the robustness to ill-conditioning of PCG\". I imagine that the paper alludes here to the dependence on the smoothed condition number $\\phi$, which remains a condition number. What is more, even in the benign setting of Corrolary 3.3, any robustness gains compared to PCG appear to be due to the restriction of the estimation to the principal $\\ell$ components of the kernel matrix. As such, there is an additional \"approximation error\" to take into account, namely, $|m_n - proj_{\\ell}(m_n)|$ which is on the order of the tail of the spectrum of $K$. \nAs such, the robustness gains also comes at the cost of a larger approximation error. This point (which should be made significantly clearer in the manuscript) suggest the importance of including an experimental investigation to further support the robustness claims (or additional theoretical results showing that equilibrating approximation and estimation error is not harmful)\n\nAs mentioned above, prior works such as Dicker et al. (2017); Lin et al. (2020; 2023) provide both theoretical and empirical evidence that a low-dimensional top eigenspace of $K + \\lambda I$ is most relevant for achieving good generalization. \nIn other words, the error restricted to the top principal components tends to align better with the test error, whereas the remaining \"tail\" approximation error tends to be associated with overfitting. \nWhile these statements are informal (and also, depend on the choice of $\\ell$), the robustness to ill-conditioning is observed empirically in our results, given the strong test error convergence of ADASAP. \nWe also note that one can freely set $\\ell=n$ both in Theorem 3.2 and in Corollary 3.3, which brings the tail approximation error to 0 and still leads to a convergence guarantee.\n\nNevertheless, we fully agree with the reviewer's point about clarifying the meaning of robustness. \nIn particular, we will discuss the role of the approximation error along the tail principal components, and how this affects the robustness of our algorithms.\n\n> **Comment 3.4** — While the paper is well-written and clear, addressing a few minor issues could improve the presentation:\n> * Some notation is not defined: $n_{\\mathrm{rhs}}$ (line 110), $K_\\lambda$ (Theorem 3.2), $\\rho$ (Line 199),  $\\mathcal Y_B$ (Line 106)\n> * The suffix for SDD-10, SDD-100 etc. is only explained in the appendix (Line 927). This should have been explained before or in the caption of Figure 1.\n> * The significance of tail averaging is only discussed in Sec 4.2.1.\n\n* Thanks for pointing out the undefined notation---we will fix this in the revision.\n\n* We also explain the suffix around line 240 in the main submission.\nHowever, we agree that explaining the suffix in the caption of Figure 1 would improve clarity, and we will update this in the revision.\n\n* We also provide a high-level overview of why we use tail averaging in section 3.1 (see lines 108-109).\nDue to limited space, we could not elaborate further in the main submission---the details of the analysis are in Appendix B.\nIf there are particular details you would like us to add to the main submission, please let us know.\n\n**References**\n\nLee H Dicker, Dean P Foster, and Daniel Hsu. Kernel ridge vs. principal component regression: Minimax bounds and the qualification of regularization operators. Electronic Journal of Statistics, 11:1022–1047, 2017.\n\nJihao Andreas Lin, Javier Antoran, Shreyas Padhy, David Janz, José Miguel Hernández-Lobato, and Alexander Terenin. Sampling from Gaussian process posteriors using stochastic gradient descent. In Advances in Neural Information Processing Systems, 2023.\n\nJunhong Lin, Alessandro Rudi, Lorenzo Rosasco, and Volkan Cevher. Optimal rates for spectral algorithms with least-squares regression over Hilbert spaces. Applied and Computational Harmonic Analysis, 48(3):868–890, 2020."}, {"review_id": "VOuwBWeKtf", "reviewer": "Reviewer_5L3E", "summary": "This paper proposes a scalable linear solver for Gaussian process (GP) posterior sampling, building on the sketch-and-project (SAP) framework of Gower and Richtárik (2015). The proposed method ADASAP integrates approximate subspace preconditioning via the randomized Nyström approximation (Tropp et al., 2017), distributed matrix-matrix products, and Nesterov acceleration. The authors also provide a convergence analysis of the SAP algorithm and evaluate ADASAP empirically on several large-scale regression datasets. Their results demonstrate the method's ability to scale to datasets with up to $10^8$ points, while showing competitive test RMSE and negative log-likelihood (NLL), and low variance with respect to different random data partitions.", "questions": "* How are RMSE and NLL computed? \n* How are the RMSE and NLL of ADASAP influenced by input data dimensionality, block size, and tail averaging?", "limitations": "The authors are do not explicitly discuss limitations of their method. I do not see any obvious potential negative social impact.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### **Quality**\n\nThe submission is technically sound and builds upon the recent GP inference literature. It combines existing components (SAP, randomized Nystrom approximations, distributed computation, and acceleration) in a coherent way, and the relation to prior work for GPs sampling is clearly explained. Although I am not familiar with the more theoretical results in the literature, the analysis of the error of the SAP algorithm seems like a meaningful contribution.\n\nOn the empirical side, the method performs well overall, particularly on the benzene, malonaldehyde, and houseelec datasets (Figure 2). However, the method’s relative performance is less convincing on the yolanda, song, and acsincome datasets (Appendix Figure 5). It would improve the presentation to move these results to the main text. Overall, the ADASAP the mean test RMSE and NLL seem to have low variance with respect to different random partitions of the data.\n\nThe presentation of experimental results in Table 1 is problematic. The use of mean values over 5 random partitions without standard errors makes it difficult to assess statistical significance even though Figures 2 and 5 suggest that this is not an issue. Furthermore, the policy of bolding results if they are within 0.01 of the best RMSE or NLL lacks rigor and should be changed with a more appropriate statistical comparison (e.g., paired sample t-tests) should be included.\n\nFurthermore, some ablations that would enhance the empirical evaluation and improve the practical relevance of the work are currently missing, including the influence on RMSE and NLL of:\n* Input feature dimensions (most of the datasets in the evaluation have $< 15$ input dimensions);\n* Block size;\n* Tail averaging.\nThese would help clarify when the method is most effective.\n\n### **Clarity**\n\nThe paper is well written and enjoyable to read. The structure and the scope of the work are clear. Connections to related work are also well presented. Code is provided though I did not executed it. Nevertheless, there are areas where clarity could be improved:\n\n* A more formal introduction to the SAP framework would benefit readers unfamiliar with this class of methods.\n* Important implementation details (e.g., computation of RMSE and mean NLL) are not reported.\n* More algorithms should be moved to the main text to make the paper more complete.\n\n### **Significance**\n\nThe work addresses the important challenge of scaling GP posterior sampling. While the algorithmic contribution is somewhat incremental, the proposed method scales to large datasets and shows good performance. The convergence analysis adds some theoretical depth.\n\nHowever, the claim that this is \"the first time that full GP inference has been scaled to a dataset of this size\" (Line 92) is overstated. As far as I understand, the large-scale experiment on the TAXI dataset ($10^8$) only computes the posterior mean (which is equivalent to the kernel ridge regression estimator) for computational budget reasons. Similar scaling has been achieved previously for the kernel ridge regression estimator, e.g., in [1], on an even larger partitioning of the TAXI dataset which exceeds $10^9$ examples. The authors should acknowledge this work, consider comparing with this method, and clarify the novelty of their approach more carefully.\n\nFurthermore, the paper does not include an ablation study of the effect of the approximate preconditioner and acceleration which would improve the paper's significance.\n\n[1] Kernel methods through the roof: handling billions of points eﬃciently, Meanti et al., 2020\n\n### **Originality**\n\nThe paper combines well-established components in a practical and effective way. The improved SAP for posterior sampling in GPs, along with a convergence analysis in RKHS norm, constitutes an incremental but meaningful contribution. \n\n### **Recommendations**\n\nOverall, I recommend accepting this work. To strengthen the paper, I suggest the following:\n\n* Improve statistical reporting in Table 1 (add standard errors, t-tests);\n* Add more algorithms to the main text;\n* Include the results in Appendix Figure 5 to the main text;\n* Clarify computation for RMSE and NLL;\n* Introduce the SAP method more formally;\n* Consider adding ablations on input dimensionality, block size, and tail averaging;\n* Moderate the claims about scalability and acknowledge relevant prior work more precisely.", "was_revised": true, "final_justification": "The authors have clarified my minor points during the discussion period. I maintain my score of 5 and recommend accepting the paper.", "rebuttal": "> **Comment 4.1** — The presentation of experimental results in Table 1 is problematic. The use of\nmean values over 5 random partitions without standard errors makes it difficult to assess statistical\nsignificance even though Figures 2 and 5 suggest that this is not an issue. Furthermore, the policy of\nbolding results if they are within 0.01 of the best RMSE or NLL lacks rigor and should be changed\nwith a more appropriate statistical comparison (e.g., paired sample t-tests) should be included.\n\nWe appreciate your suggestions for improvement.\nIn the revision we will provide standard errors, and we will use a paired t-test for comparing the methods in Table 1. \n\n> **Comment 4.2** — Furthermore, some ablations that would enhance the empirical evaluation and improve the practical relevance of the work are currently missing, including the influence on RMSE and NLL of:\n> * Input feature dimensions (most of the datasets in the evaluation have $< 15$ input dimensions);\n> * Block size;\n> * Tail averaging. These would help clarify when the method is most effective.\n\n* The dimensions of the datasets are: yolanda $(d = 100)$, song $(d=90)$ benzene $(d = 66)$, malonaldehyde $(d=36)$, acsincome $(d=9)$, houseelec $(d=9)$, and taxi $(d=9)$, so the majority of the inputs do not have feature dimensions $< 15$. \nWe will be sure to include this information in the revision.\n\n* The theory developed in our paper shows that the convergence rate improves as the blocksize increases.\nHowever, the improvement in convergence rate must be balanced with increased computational and memory costs.\nThis phenomenon has been verified in many prior SAP papers. \nHence, for brevity we have not studied this here, but we will point to the references that show this in the revision.\nIn practice, we recommend using the blocksize that maximizes GPU utilization while not running out of memory.\n\n* As discussed in Section 4.2.1, we believe the value of tail averaging is primarily theoretical, as has been observed in the SGD literature (Shamir and Zhang, 2013; Johnson and Zhang, 2013), where the last iterate often delivers comparable performance to the averaged iterate.\nIn the revision, we will add experiments to show the impact of tail averaging on convergence.\n\n> **Comment 4.3** — Nevertheless, there are areas where clarity could be improved:\n> * A more formal introduction to the SAP framework would benefit readers unfamiliar with this class of methods.\n> * Important implementation details (e.g., computation of RMSE and mean NLL) are not reported.\n> * More algorithms should be moved to the main text to make the paper more complete.\n\n* Section 2 provides a brief introduction to SAP and its properties most relevant for our setting.\nWe would like to have provided more intuition for SAP, but could not owing to length constraints.\nFinally, section A.2 provides a detailed literature review of SAP methods, from which unfamilar readers can learn more.\n\n* Our computation of these quantities follows Lin et al. (2023; 2024).\nThe provided codebase also provides all the code used to compute these quantities in our experiments.\nThe precise mathematical formulae we used are provided in response to Comment 4.6.\n\n* Unfortunately, due to length constraints were forced to push the description of many algorithmic subroutines to the appendix. \nWe agree this is not necessarily ideal.\nIf accepted, we get an extra page, and will try and bring more of these algorithms forward to the main paper.\n\n> **Comment 4.4** — However, the claim that this is \"the first time that full GP inference has been scaled to a dataset of this size\" (Line 92) is overstated. As far as I understand, the large-scale experiment on the TAXI dataset ($10^8$) only computes the posterior mean (which is equivalent to the kernel ridge regression estimator) for computational budget reasons. Similar scaling has been achieved previously for the kernel ridge regression estimator, e.g., in [1], on an even larger partitioning of the TAXI dataset which exceeds $10^9$ examples. The authors should acknowledge this work, consider comparing with this method, and clarify the novelty of their approach more carefully.\n\nWhen we refer to full GP inference, we mean performing GP inference on the entire dataset, without using inducing points.\nWhile Falkon in [1] runs on $10^{9}$ samples, they use $10^{5}$ inducing points, which means they are processing a kernel matrix of size $10^5 \\times 10^9$.\nOn the other hand, the kernel matrix that we process has size $3.31 \\cdot 10^8 \\times 3.31 \\cdot 10^8$---which is three orders of magnitude larger than the kernel matrix processed by Falkon.\nFurthermore, inducing points methods for GPs and kernel ridge regression often perform worse than methods that solve the full problem (Gardner et al., 2018; Lin et al., 2023; 2024; Rathore et al., 2025).\nTherefore, we do not think it is necessary to perform a comparison with Falkon at this time.\n\n> **Comment 4.5** — Furthermore, the paper does not include an ablation study of the effect of the approximate preconditioner and acceleration which would improve the paper's significance.\n\n* Approximate preconditioner: Incorporating the Nyström approximation allows ADASAP to scale to larger blocksizes.\nWe did not run ADASAP without the Nyström approximation, since it would not scale to our largest datasets.\nFor example, on taxi $(n = 3.31 \\cdot 10^8)$, we run ADASAP with blocksize $b = 1.65 \\cdot 10^5$.\nComputing the Cholesky decomposition to solve the $b \\times b$ linear system with $K_{\\mathcal B \\mathcal B} + \\lambda I$ would require $\\mathcal O(b^3)$ compute and $\\mathcal O(b^2)$ memory, which is not practical for this blocksize.\nSince the \"full\" preconditioner does not scale to this blocksize, we omitted it from our experiments.\n\n* Acceleration: Several works have shown that incorporating acceleration into KRR solution methods leads to consistent practical improvements (Tu et al., 2017; Rathore et al., 2025).\nTherefore, in the interest of time, we did not run experiments without acceleration.\n\n> **Comment 4.6** — How are RMSE and NLL computed?\n\nLet the test set be $\\lbrace (x_i, y_i) \\rbrace_{i = 1}^{n_{\\mathrm{tst}}}$.\nFurthermore, let $m_n(\\cdot)$ be the mean of the GP posterior.\n\nWe compute the RMSE using the expression $\\sqrt{\\frac{1}{n_{\\mathrm{tst}}} \\sum_{i = 1}^{n_{\\mathrm{tst}}}(m_n(x_i) - y_i)^2}$.\n    \nFor the mean NLL, we compute 64 samples from the posterior GP using pathwise conditioning, and we use those samples to compute empirical variance estimates $\\sigma_i^2$ for each point $i$ in the test set.\nWe then compute the negative log-likelihood for each point $i$ in the test set using $\\mathcal N(m_n(x_i), \\sigma_i^2)$ as the \"reference\" distribution.\nFinally, we take the mean of these negative log-likelihoods to get the mean NLL.\n\nOur approach for computing RMSE and mean NLL is identical to that of Lin et al. (2023; 2024).\nPlease see our code for more details.\n\n> **Comment 4.7** — How are the RMSE and NLL of ADASAP influenced by input data dimensionality, block size, and tail averaging?\n\nPlease see our response to Comment 4.2.\n\n> **Comment 4.8** — Contributions 3 and 4 are not significantly different. Perhaps they could be merged into a single \"Experimental Results\" point.\n\nWe originally separated these contributions because they were addressing different experimental settings (GPs and Bayesian optimization) which have different use cases.\nHowever, we will consider merging these in the revision to improve clarity.\n\n> **Comment 4.9** — Line 72: The claim is inaccurate—conditioning requires cubic time in the number of data points.\n\nWe are confused by the point the reviewer is making---in general, computing the square root of a $n \\times n$ matrix requires $\\mathcal O(n^3)$ time.\nPlease let us know if we have misinterpreted your point.\n\n> **Comment 4.10** — Line 135: The formulation “averaging the second half” is unclear; please clarify.\n\nWhen we say ``averaging the second half'' we are referring to the tail averaging in Algorithm 1 (see the if statement in the algorithm block).\nWe will make sure that this is clarified in the revision.\n\n**References**\n\nJacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPyTorch:\nBlackbox matrix-matrix Gaussian process inference with GPU acceleration. In Advances in\nNeural Information Processing Systems, 2018.\n\nRie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance\nreduction. In Advances in Neural Information Processing Systems, 2013.\n\nJihao Andreas Lin, Javier Antoran, Shreyas Padhy, David Janz, José Miguel Hernández-Lobato,\nand Alexander Terenin. Sampling from Gaussian process posteriors using stochastic gradient\ndescent. In Advances in Neural Information Processing Systems, 2023.\n\nJihao Andreas Lin, Shreyas Padhy, Javier Antoran, Austin Tripp, Alexander Terenin, Csaba Szepesvari, José Miguel Hernández-Lobato, and David Janz. Stochastic gradient descent for Gaussian\nprocesses done right. In International Conference on Learning Representations, 2024.\n\nPratik Rathore, Zachary Frangella, Jiaming Yang, Michal Dereziński, and Madeleine Udell. Have\nASkotch: A neat solution for large-scale kernel ridge regression. arXiv preprint 2407.10070, 2025.\n\nOhad Shamir and Tong Zhang. Stochastic gradient descent for non-smooth optimization: Convergence results and optimal averaging schemes. In Proceedings of the 30th International Conference on Machine Learning, 2013.\n\nStephen Tu, Shivaram Venkataraman, Ashia C. Wilson, Alex Gittens, Michael I. Jordan, and\nBenjamin Recht. Breaking locality accelerates block Gauss-Seidel. In Proceedings of the 34th\nInternational Conference on Machine Learning, 2017."}]

Author rebuttal:

Meta-review: The reviewer consensus is that this is a good paper. Reviewers praise high performance, theoretical grounding in terms of convergence results, and the ability to implement in a distributed manner as strengths. Reviewers also praise clarity, and describe the work as well-written. Most of the weaknesses have to do with minor typos, and some presentation details with how ill-conditioning is discussed (Reviewer WCDC, who gave the most critical score). From the text of the reviews, I infer that the reviewers see the strengths as significantly outweighing the weaknesses.

The reviewer consensus is that it should be accepted, and I concur.
Decision: Accept (poster)