Improving Linear System Solvers for Hyperparameter Optimisation in Iterative Gaussian Processes
Authors: Jihao Andreas Lin
Shreyas Padhy
Bruno Kacper Mlodozeniec
Javier Antoran
José Miguel Hernández-Lobato
Venue: NeurIPS 2024 poster (2024)
Keywords: Gaussian process
marginal likelihood
iterative linear system solver
pathwise
warm start
conjugate gradient
alternating projection
stochastic gradient descent
Primary area: probabilistic_methods
TL;DR: We accelerate iterative linear system solvers in the context of marginal likelihood optimisation for Gaussian processes.

Abstract: Scaling hyperparameter optimisation to very large datasets remains an open problem in the Gaussian process community. This paper focuses on iterative methods, which use linear system solvers, like conjugate gradients, alternating projections or stochastic gradient descent, to construct an estimate of the marginal likelihood gradient. We discuss three key improvements which are applicable across solvers: (i) a pathwise gradient estimator, which reduces the required number of solver iterations and amortises the computational cost of making predictions, (ii) warm starting linear system solvers with the solution from the previous step, which leads to faster solver convergence at the cost of negligible bias, (iii) early stopping linear system solvers after a limited computational budget, which synergises with warm starting, allowing solver progress to accumulate over multiple marginal likelihood steps. These techniques provide speed-ups of up to $72\times$ when solving to tolerance, and decrease the average residual norm by up to $7\times$ when stopping early.

Improving Linear System Solvers for Hyperparameter Optimisation in Iterative Gaussian Processes
Jihao Andreas Lin 1.2
Shreyas Padhy!
Bruno Mlodozeniec
Javier Antorán!:3 1.3
José Miguel Hernández-Lobato-
University of Cambridge
MPI for Intelligent Systems, Tübingen SAngström Al
Abstract
Scaling hyperparameter optimisation to very large datasets remains an open prob- lem in the Gaussian process community. This paper focuses iterative methods, which use linear system solvers, like conjugate gradients, alternating projections stochastic gradient descent, to construct estimate of the marginal likelihood gradient. We discuss three key improvements which are applicable across solvers: pathwise gradient estimator, which reduces the required number of solver iterations and amortises the computational cost of making predictions, (ii) warm starting linear system solvers with the solution from the previous step, which leads to faster solver convergence at the cost of negligible bias, (iii) early stopping linear system solvers after a limited computational budget, which synergises with warm starting, allowing solver progress to accumulate over multiple marginal likelihood steps. These techniques provide speed-ups -pp of up 72x when solving to tolerance, and decrease the average residual norm by up to 7x when stopping early.
Introduction
Gaussian processes 22 (GPs) are a versatile class of probabilistic machine learning models which are used widely for Bayesian optimisation black-box functions (24), climate and earth sciences (6). and data-efficient learning in robotics and control (6]. However, their effectiveness depends on good estimates of hyperparameters, such kernel length scales and observation noise. These quantities are typically learned by maximising the marginal likelihood, which balances model complexity with training data fit. In general, the marginal likelihood is non-convex function of the hyperparameters and evaluating its gradient requires inverting the kernel matrix. Using direct methods, this requires compute and memory resources which are respectively cubic and quadratic in the number of training examples. This is intractable when dealing with large datasets of modern interest.
Methods improve the scalability of Gaussian processes can roughly be grouped into categories. Sparse methods [20] 27 12] approximate the kernel matrix with low-rank surrogate, which is cheaper invert. This reduced flexibility may result in failure to properly fit increasingly large or sufficiently complex data [5). On the other hand, iterative methods [1I) express GP computations in terms of systems of linear equations. The solution these linear systems is approximated up specified numerical precision with linear system solvers, such as conjugate gradients (CG) [I]
(30]. alternating projections (AP) [23 28 33.. or stochastic gradient descent (SGD) [15 16). These
methods allow for a trade-off between compute time and accuracy. However, convergence can be slow the large data regime, where system conditioning is often poor.
In this paper, we focus on iterative GPs and identify techniques, which were important to the success of previously proposed methods, but did not receive special attention in the literature. Many of these amount to amortisations which leverage previous computations to accelerate subsequent ones. We analyse and adapt these techniques, and show that they can be applied accelerate different linear solvers, obtaining speed-ups of up to 72x without sacrificing predictive performance (see Figure
38th Conference on Neural Information Processing Systems (NeurIPS 2024).

---

Conjugate Gradients
Alternating Projections
Stochastic Gradient Descent
Relative Time
100%
100%
100%
80%
80%
80%
60%
60%
60%
40%
40%
20%
20%
20%
0%
0%
0%
POL
BIKE
ELEV
POL
BIKE
Datasct Name
Dataset Name
Dataset Name
Standard
Pathwise
Standard + start
Pathwise
start
Hatched areas represent the fraction of the total runtime which used by the linear system solver
Figure 1: Comparison of relative runtimes for different methods, linear system solvers, and datasets. The linear system solver (hatched areas) dominates the total training time (coloured patches). The pathwise gradient estimator requires less time than the standard estimator. Initialising at the previous solution (warm start) further reduces the runtime of the linear system solver for both estimators.
In the following, we summarise our contributions:
We introduce a pathwise estimator of the marginal likelihood gradient and show that, under real-world conditions. the solutions the linear systems required by this estimator are closer to the origin than those of the standard estimator, allowing our solvers to converge faster. Additionally, these solutions transform into samples from the GP posterior without further matrix inversions, amortising the computational costs of predictive posterior inference.
We propose to warm start linear system solvers throughout marginal likelihood optimisation by reusing linear system solutions to initialise the solver in the subsequent step. This results in faster convergence. Although this technically introduces bias into the optimisation, we show that, theoretically and empirically, the optimisation quality does not suffer.
We investigate the behaviour of linear system solvers on limited compute budget, such that reaching the specified tolerance is not guaranteed. Here, warm starting allows the linear system solver to accumulate solver progress across marginal likelihood steps, progressively improving the solution quality of the linear system solver despite early stopping.
We demonstrate empirically that the methods above either reduce the required number of iterations until convergence without sacrificing performance or improve the performance a limited compute budget hinders convergence. Across different UCI regression datasets and linear system solvers, we observe average speed-ups of up 72x when solving until the tolerance is reached, and increased performance when the compute budget is limited.
Source code available at:
2 Gaussian Process Regression and Marginal Likelihood Optimisation
Formally, a GP is stochastic process f R. such that, for any finite subset { C X. the set of random variables (f(x;;) is jointly Gaussian. In particular, f is uniquely identified by mean function uu.)) = B[](-)] and positive-defintte kernel function (;) = Co)(f()), f()) with kernel hyperparamcters 0. We use a Matémn3/2 kernel with length scales per dimension and scalar signal scale and write ~ GP k) to express that f is a GP with mean l and kernel
For GP regression, let the training data consist of n inputs x C X and targets y E R". We consider the Bayesian model y = fz,) + e, where each & (O, o²) i..d. and f GP k). We assume 0 without loss of generality. The posterior of this model is fly GP( k(1)): with
( (-) k(,, æ; (k( x;;
y.
( /:
k(, æ; (k(x, x; 0) kk)
(2)
where k(,, x; -), +) and k(æ, x; +) refer to pairwise evaluations, resulting in a n row vector, n X 1 column vector and x matrix respectively.
Pathwise Conditioning Wilson al. 31 32 express GP posterior sample as random function
(flv)()) = f(.) + k(, x 1 (16 x. xz 0) - (f(x) + +c),
(3)

---

where e ~ ~N(, is random vector. f ~ GP(O, k) is a zero-mean prior function sample, and
fx) is its evaluation at the training data. Following previous work [31 32 15 16). we efficiently approximate the prior function sample using random features [21 25 (see Appendix B for details).
Using pathwise conditioning. single linear solve suffices to evaluate a posterior function sample we amortise the cost of this single
at arbitrary locations without further linear solves. In Section 3
linear solve during marginal likelihood optimisation to obtain posterior samples efficiently.
The Marginal Likelihood and Its Gradient With hyperparameters O = { o) and regularised
kernel matrix Hg = k(x,, æ: 2) + o²I E R=*" the marginal likelihood C as function of 0 and its gradient Vo, & with respect to 0% can be expressed as
1 log det Hg
n2
((0) 2y 9 y
log 2T,
(4)
2
2
1 2 (Ha²y) y) HH
0Hg
Vo, (())
H
y
2 ir
H 0
(5)
00k
00
where the partial derivative of Hg with respect to 02 is x matrix. We assume n is too large to compute the inverse or log-determinant of Hg and iterative methods are used instead.
2.1 Hierarchical View of Marginal Likelihood Optimisation for Iterative Gaussian Processes
Marginal likelihood optimisation for iterative GPs consists of bi- level optimisation, where the outer loop maximises the marginal likelihood 44) using stochastic estimates of its gradient (5). Comput- ing these gradient estimates requires the solution systems of linear equations. These solutions are obtained using an iterative solver in the inner loop. Figure Fillustrates this three-level hierarchy.
Outer-Loop Optimiser L-BPGS, Adam
Gradient Estimator
Hutchinson trace estimator
Outer-Loop Optimiser
The outer-loop optimiser maximises the
Linear System Solver e.g. conjugate gradients
marginal likelihood L (4) using its gradient (5). Common choices
are L-BFGS (3), when exact gradients are available, and Adam (14) in the large-data setting, when stochastic approximation is required. We consider the case where gradients are stochastic and use Adam.
Figure Marginal likelihood optimisation for iterative GPs.
Gradient Estimator The gradient (5) involves two computationally expensive components: linear solves against the targets Hg y and the trace term ir (H, g (H,/00s), An unbiased estimate of the latter can be obtained using probe vectors and Hutchinson's 's trace estimator (13),
HHo 00k
HH, 00
HHg 09 2;-
H
Ez
z H
H
(6)
where the probe vectors 2j E R" satisfy Vj : =[ = I. and = "Ho is obtained using a linear solve. We refer to this as the standard estimator and set 64, unless otherwise specified.
Linear System Solver Substituting the trace estimator (6) back into the gradient ( we obtain an unbiased gradient estimate in terms of the solution batch of systems of linear equations,
Ho (Dy: V1, U..
y. 21. 2.
which share the same coefficient matrix Hg. Since Ho is dositive-definite, the solution v = H to the system He b can be obtained by finding the unique minimiser of the quadratic objective
1 =u H2 ubb.
arg min
(8)
u 2
facilitating the use of iterative solvers. Most popular in the GP literature are conjugate gradients (CG) 28
[9] (2), alternating projections (AP) [23]
2 33) and stochastic gradient descent (SGD) [5] 16].
We consider these in our study and provide detailed descriptions of them in Appendix B Sol vers are
often run until the relative residual norm b Hau|
|/||b] reaches a certain tolerance 30 18 33
We set 0.01, following Maddox et al. [18 The linear system solver in the inner loop dominates the computational costs of marginal likelihood optimisation for iterative GPs, as shown in Figure] Therefore, improving linear system solvers is the main focus of our work.
3

---

Initial Distance
Solver Iterations
Top Eigenvalue Hg"
Noisc Precision -2
6M
240k
1500
1500
POL
4M
160k
1000
1000
2M
80k
500
500
0
0
24k
18
18
ELEV
110k
16k
12
12
55k
8k
0
0
25
50
75
100
25
50
75
100
25
50
75
100
0
25
50
75
100
Marginal Likelihood Steps
Standard
Pathwise
Standard (thcoretical)
Pathwise (thcoretical)
Figure 3: On the POL and ELEVATORS datasets, the pathwise estimator results in lower RKHS dis- lance (12) between solver initialisation and solution, as predicted by theory [14 (left). This results in fewer AP iterations until reaching the tolerance (left middle). When using the standard estimator, the initial distance follows the top eigenvalue of H (right middle), which is strongly related to the noise precision (right). The latter tends to increase during marginal likelihood optimisation when fitting the data. The effects are greater on POL due the higher noise precision.
Pathwise Estimation of Marginal Likelihood Gradients
We introduce the pathwise estimator, an alternative to the standard estimator (6) which reduces the required number of linear system solver iterations until convergence (see Figure 3 Additionally, the estimator simultaneously provides us with posterior function samples via pathwise conditioning, hence the name pathwise estimator. This facilitates predictions without further linear solves.
We modify the standard estimator (6) absorb Hgy
into the distribution of the probe vectors
HHo
AHe 08
HHg
HHe
H
H
H
Eg
ET
Ejy
00
08
08
j=1
where Vj 2[ = H Probe vectors & with the desired second moment can be obtained as f(x) ~ NO, k(æ, 0);
= g (0, Ho) = 2 = (/ H,*)
(10)
c~ C~N(00.7T)
where g = f(x) + e. Akin to the standard estimator in Section 2.1 we obtain By and ž, by solving He Ug. 21...., 2, y. E.. E.
(11)
Initial Distance to the Linear System Solution Under realistic conditions, the pathwise estimator moves the solution of the linear system closer to the origin. To show this, we consider the generic linear system Hgu = and measure the RKHS distance between the initialisation Unitt and the solution
H₂¹ b as || u uIlio With Unnt
O, which is standard 9 30 15
33
16
Uimit u Ili ||u||| li, = u Hou = 6 H2 H 6 H bb.
(12)
Since b is random vector (= in 7 and in [1
analyse the expected squared distance
H b] [tr (b HH b)] = [tr (bb (bbTH )] = tr [bb]] Hg').
(13)
For the standard estimator (6). we substitute b = z with E 2z T = I, yielding
E []] Vinit
Ho. tr (B [zz H = tr (IHH')) tr (H²²).
(14)
For the pathwise estimator
9 we substitute = c with [Ec He. yielding
E [[] - u I = tr (E [EE H
tr (H, H₂²) tr (I) =
(15)
The initial distance for the standard estimator is equal to the trace of H whereas it is constant for the pathwise estimator. Figure Figurr@iilussratte illustrates that this trace follows the top eigenvalue, which roughly matches the noise precision. As the model fits the data better, the noise precision increases, increasing the initial distance for the standard but not for the pathwise estimator. practice, the latter leads to faster solver convergence, especially for problems with high noise precision (see Table

---

Test Log-Likelihood
Time (min): CG
Time (min): AP
Time (min): SGD
0.20
4.8
48
80
0.15
3.6
36
60
0.10
2.4
24
40
0.05
1.2
12
20
0.00
0.0
0
16
64
1024
25
50
75
100
0
25
50
75
100
0
25
50
75 100
Number Samples
Marginal Likelihood Steps
Absolute difference between pathwise
Number of probe vectors for pathwise estimator
64 (default)
conditioning and predictions
16
Figure 4: On the POL dataset, increasing the number of posterior samples improves the performance of pathwise conditioning until diminishing returns start to manifest with more than 64 samples (left). Furthermore, with 4x as many probe vectors, the total cumulative runtime only increases by around 10% because the computational costs are dominated by shared kernel function evaluations (right).
Amortising Linear Solves for Optimisation and Prediction The name of the pathwise estimator comes from the fact that solving the linear systems (TT) provides us with all of the terms we need to construct set of posterior samples via pathwise conditioning (3. Each of these is given by
(+) = ʃ()) k(,, x; 0) Hp' (y - <) = s() + k(, ax (), 2).
(16)
We can use these to make predictions without requiring any additional linear system solves.
How Many Probe Vectors and Posterior Samples Do We Need? In the literature 9
18 2 33
it is common to use s S 16 probe vectors for marginal likelihood optimisation. However, larger number of posterior samples, around s = 64, is necessary to make accurate predictions 2 (see 4. Thus, to amortise linear system solves across marginal likelihood optimisation and prediction, we must use the same number of probes for both. Interestingly, shown in Figure 4 using 64 instead of probe vectors only increases the runtime by around 10% because the computational costs are dominated by kernel function evaluations, which are shared among probe vectors.
15
Estimator Variance The standard estimator with Gaussian probe vectors and the pathwise estimator have the same variance if Hg and 0He 000; commute with each other (see Appendix A.IIp. There has been work developing trace estimators with lower variance 19 88. however, we did not pursue these as we find variance to be sufficiently low, even when relying on only 16 probe vectors.
Approximate Prior Function Samples Using Random Features In practice, the pathwise estima- tor requires samples from the prior f (0, k). which are intractable for large datasets without the of random features (3 32] In Figure 5 we show that, despite using random features, most of the time the marginal likelihood optimisation trajectory of the pathwise estimator matches the trajectory of exact optimisation using Cholesky factorisation and backpropagation. Further, we confirm that deviations of the pathwise estimator are indeed due to the use of random features by demonstrating that we can remove these deviations using exact samples from the prior instead.
Pathwise Trajectory
Length Scale 9
Length Scale 12
Length Scale 17
40%
2.6
2.6
2.6
30%
2.2
2.2
2.2
20%
1.8
1.8
1.8
10%
1.4
1.4
0%
1.0
1.0
0.5
0.5
25
75
100
0
Difference to Exact
Exact
Marginal Likelihood Steps
Pathwise estimator with random features
Pathwise estimator with exact prior
Figure 5: Across all datasets and marginal likelihood steps, most hyperparameter trajectories of the pathwise estimator rarely differ from exact optimisation, as shown by the histogram illustrating the differences between hyperparameters (left). On selected length scales of the ELEVATORS dataset, the pathwise estimator deviates due to the use of random features to approximate prior function samples. With exact samples from the prior, the pathwise estimator matches exact optimisation again (right).
5

---

Mean
Probe Vectors
Initial Distance to Solution
Second Eigendircctoon
1.0
1.2
Mran
0.6
0.G
20
0.2
Probe
-0.2
-1
-0.2
0.2
0.6
0
25
50
75
100
First Eigendirection
First. Eigendirection
Marginal Likelibood Steps
Initialisc at zero (default)
Initialisc at previous solution (warm start)
Solution
Figure 6: Two-dimensional cross secsiseccionn along top eigendirections of the inner-loop quadratic objective after 20 marginal likelihood steps on the POL dataset. The current solution is placed at the origin of coordinates (left and middle). Warm starting significantly reduces the initial root-mean- square RKHS distance to the solution throughout marginal likelihood optimisation (right).
CG
AP
SGD
Test. L.H.H
225k
210k
Solver Iterations
POL
150
150k
140k
1.27
75
75k
70k
1.26
0
1.25
48
24k
6k
0.37
ELEV
32
0.38
16
8k
2k
0.39
0
-0.40
25
50
75
100
25
50
75
100
0
25
50
75
100
CG
AP
SGD
Marginal Likelihood Steps
Standard estimator
Standard estimator + warm start
Figure 7: Required number of solver iterations to the tolerance 0.01 during marginal likelihood optimisation on the POL and ELEVATORS datasets. Warm starting with the previous solution reduces the required number of iterations to reach the tolerance without sacrificing predictive performance.
4 Warm Starting Linear System Solvers
Linear system solvers are typically initialised at zero 30
15
16 However, because the 33
outer-loop marginal likelihood optimisation does not change the hyperparameters much between consecutive steps, we expect that the solution to inner-loop linear systems also does not change much between consecutive steps (see Appendix A.2 for more formal argument). Therefore, we suggest to warm start linear system solvers by initialising them at the solution of the previous [17]. This requires that the targets of the linear systems, z; or E,, are resampled throughout optimisation, which can introduce bias [5]. However, we find that warm starting consistently provides gains across all linear system solvers for both the standard and the pathwise estimator, and that the bias is negligible.
Visualising Warm Starts Segurelvssuualseee the two top eigendirections of the inner-loop quadratic objective on POL. Throughout traiming, warm starting solvers at the solution to the previous linear system results substantially smaller initial distance the current solution.
Efffects on Linear System Solver Convergence Reducing the initial RKHS distance the solution reduces the required number of solver iterations until the tolerance - 0.01 is reached for all solvers and all five datasets, as shown in Figure 7 Table and Appendix C However, the effectiveness depends on the solver type. CG is more sensitive to the direction of descent rather than the distance the solution because uses line searches to take big steps. only obtains 2.1x speed-up on average. AP and SGD benefit more, with average speed-ups of 18.9 and 5. respectively.
Notable exceptions are Artemev et al. 3 who warm start vn in sparse lower bound on L. and Antorán et al. 2 who warm start stochastic gradient descent solver for miite-dimensionl linear models.

---

Table Test log-likelshoods, total training times, and average speed-up among datasets for CG, AP, and SGD after 100 outer-loop marginal likelihood steps with learning rate of 0.1. We consider five datasets with < 50k; which allows us to solve to tolerance, and report the mean over 10 data splits.
path
warm
Test Log-Likelihood
Total Time (min)
Average Speed-Up
wise
POL
ELEV
BIKE
PROT
KEGG
POL
ELEV
BIKE
PROT
KEGG
1.27
-0.39
2.15
-0.59
1.08
4.83
1.58
5.08
29.9
28.0
1.27
-0.39
2.07
-0.62
1.08
3.96
1.49
4.41
20.0
26.4
1.2
CC
1.27
-0.39
2.15
-0.59
1.08
2.28
1.03
2.74
11.5
12.8
2.1
1.27
-0.39
2.06
-0.62
1.08
2.47
1.00
3.07
13.7
13.0
1.9
1.27
-0.39
2.15
-0.59
493.
77.8
302.
131.
24
AP
1.27
-0.39
2.07
-0.62
1.08
27.9
1.67
19.9
16.4
211.
5,4 x
1.27
-0.39
2.15
-0.59
1.08
44.0
36.4
35.1
55.8
491.
> 18.9
1.27
-0.39
2.06
-0.62
1.08
3.90
1.21
5.40
12.3
14.0
72.1
1.27
-0.39
2.15
-0.59
1.08
139.
5.54
412.
75.2
620.
SGD
1.27
-0.39
2.07
-0.63
1.08
73.6
4.58
156.
24.0
412.
2.1
1.27
-0.39
2.15
-0.59
1.08
26.5
1.22
74.3
11.2
168.
1.27
-0.39
2.06
-0.62
1.07
17.9
1.14
64.2
11.9
58.7
7.2 X
Warm Start Trajectory
Length Scale
Length Scale 12
Length Scale
60%
2.6
2.6
2.6
45%
2.2
2.2
30%
1.8
1.8
15%
14
1.4
0%
1.0
1.0
1.0
0.2
0.2
25
50
75
100
50
75
100
0
25
50
75
100
Difference to Exact
Marginal Likelihood Steps
Standard estimator + warm start
Exact
Figure 8: Across all marginal likelihood steps and datasets, warm starting results in hyperparameter trajectories which barely differ from exact optimisation, as shown by the histogram (left). On the same selected length scales from Figure warm starting matches exact optimisation (right).
Does Warm Starting Introduce Bias? A potential concern when warm starting is that the latter introduces bias into the optimisation trajectory because the linear system targets not resampled throughout optimisation. Although individual gradient estimates are unbiased, estimates are correlated along the optimisation trajectory. In fact, after fixing the targets, gradients become deterministic and it is unclear whether the induced optimum converges to the true optimum 2 Fortunately, one can show that the marginal likelihood the optimum implied by these gradients will converge in probability to the marginal likelihood of the true optimum.
Theorem (informal) Under reasonable assumptions, the marginal likelihood L of the hyperparam- eters obtained by maximising the objective implied by the warm-started gradients O will converge in probability to the marginal likelihood of a true maximum 0*: L((*) 2, L())) + cc.
See Appendices] A.3 and for details. In practice, small number of samples seems to be sufficient. In Appendix C. we illustrate that optimisation trajectories of warm-started solvers are almost identical trajectories obtained by non-wand non-warm-started solvers across solver types and datasets.
Warm Starting the Pathwise Estimator One advantage of the pathwise estimator from Section the reduced RKHS distance between the origin and the solution. When warm starting, the inner-loop solver no longer initialises at the origin, and thus, one may be concerned that we lose this advantage. However, empirically, this is not the case. As shown in Table combining both techniques further accelerates AP and SGD, reaching 72.1x 1x and 7.2x average speed-ups across our datasets relative to the standard estimator without warm starting. Furthermore, since we run solvers until reaching the tolerance, the predictive performance is almost identical among all methods and solvers.
"The concern might be likened to how pointwise convergence of integrable functions does not always imply convergence of the integrals of those functions, potentially biasing the optima of the limit of the integrals.
7

---

CG
AP
SGD
Test L.H.
6
0.6
1.2
3.1
Residual Norm
4
0.4
0.8
1.3
2
0.2
0.5
0
0.0
0.0
2.3
1.2
0.12
0.6
1.5
0.8
0.08
1.3
0.4
0.04
0.2
1.1
0.0
0.00
0.0
0.9
0
25
50
75
100
25
50
75
100
25
50
75
100
CG
AP
SGD
Marginal Likelihood Steps
Standard
Pathwise
Standard + warm start
Pathwise + warm start
Maximum solver epochs per marginal likelihood step:
10
20
30
40
50
Figure 9: Relative residual norms of the probe vector linear systems each marginal likelihood step on the POL dataset when solving until the tolerance or a maximum number of solver epochs reached. Increasing the compute budget generally reduces the residual norm. Given the same compute budget, the pathwise estimator reaches lower residual norms than the standard estimator. Adding warm starts further reduces the residual norm for both estimators. However, the final test log-likel legllike does not always match the residual norm. Surprisingly, good predictive performance can be obtained even if the residual norm is much higher than the tolerance T = 0.01.
5 Solving Linear Systems on a Limited Compute Budget
Our experiments far have only considered relatively small datasets with n < 50k, such that inner- loop solvers can reach the tolerance in a reasonable amount of time. However, on large datasets, where linear system conditioning may be poor, reaching a low relative residual norm can become computationally infeasible. Instead, linear system solvers are commonly given a limited compute budget. Gardner et al. [ limit the number of CG iterations
20, Wu et al. [33 use 11 epochs of
AP, Antorán et al. (2) run SGD for 50k iterations, and Lin et al. [5] 16] run SGD for 100k iterations. While effective for managing computational costs, is not well understood how early stopping before reaching the tolerance affects different solvers and marginal likelihood optimisation. Furthermore, it unclear whether certain tolerance is required to obtain good downstream predictive performance.
The Effects of Early Stopping We repeat the experiments from Table but introduce limited compute budgets: 10, 20, 30, 40 or 50 solver epochs, where one epoch refers to computing each value in He once (see Appendix B for details), In this setting, solvers terminate upon either reaching the relative residual norm tolerance or when the compute budget exhausted, whichever occurs first.
In Figure 9 we illustrate the relative residual norms reached for each compute budget on the POL dataset (see Figures| 14 17 in Appendix for other datasets). In general, the residual norms increase as Hg becomes more ill-conditioned during optimisation, and as the compute budget is decreased. The increase in residual norms is much larger for CG than the other solvers, which is consistent with previous reports of CG not being amenable to early stopping [15). AP seems to behave slightly better than SGD under limited compute budget. Both the pathwise estimator and warm starting combine well with early stopping, reaching lower residual norms when using budget of 10 solver epochs than the standard estimator without warm starting using a budget of 50 solver epochs.
In terms of predictive performance, we see that CG with the standard estimator and no warm starting suffers the most from early stopping. Changing the pathwise estimator and warm starting recovers good performance most of the time. SGD also shows some sensitivity to early stopping, but there to be stronger correlation between invested compute and final performance. Surprisingly, AP generally achieves good predictive performance even on the smallest compute budget, despite not reaching the tolerance of 7 0.01. Overall, the relationship between reaching low residual norm and obtaining good predictive performance seems to be weak. This is an unexpected yet interesting observation, and future research should investigate the suitability of the relative residual norm as metric to determine solver convergence.
Because kernel function evaluations dominate the computational costs of linear system solvers, this results similar time budgets methods while preventing compute wastage due to mesed stopping.
8

---

Conjugate Gradients Test LLH
Alternating Projections Test LLH
Stochastic Gradient Descent
Res. Norm
Res. Norm
Res. Norm
Test LLH
BUZZ HOUSE
1.8
0.09
2.6
0.20
1.8
144
WWW
0.0G
1.3
0.15
1.0
1.0
0.03
0.10
0.2
0
5 10:15
5 15
10
0
5 10
0
5 10
10
15
4.4
0.11
-0.03
0.6
0.0
2.7
-7
0.00
0.08
0.3
-0.1
1.0
144
0.01
0.13
-0.2
0 10.00 30
10 20 30
10
20 30
0
10 20 30
0 10 20 30
10 20 30
SONG
2.40
1..
0.08
11.1
0.10
1.11
1.25
1.6
0.04
1.12
0.05
1.12
0.10
2.2
0.00
1.13
0.00
1.13
0 10
20
30
10 20 30
10
20 30
0 10 20 30
0 10 20
30
10 20 30
Marginal Likelihood Steps
Pathwisc estimator
Puthwisc estimator
+ warm start
Figure 10: Relative residual norms and test log-likelihoods during marginal likelihood optimisation on large datasets using the pathwise estimator. Warm starting allows solver progress to accumulate over multiple marginal likelihood steps, leading to decreasing residual norms. Without warm starting, residual norms tend to remain similar or increase during optimisation. Despite reaching significantly lower residual norms, the predictive performance does not always improve, akin Figure
Demonstration on Large Datasets After analysing early stopping on small datasets, we now turn to larger UCI datasets 391k < n 1.8M, where solving until reaching the tolerance becomes computationally infeasable. Thus, we introduce compute budget of 10 solver epochs per marginal
likelihood step. Hyperparameters are initialised with the heuristic of Lin et al. [15] and optimised using a learning rate of 0.03 for 30 Adam steps (15 for HOUSEELECTRIS due to high computational costs). We use the pathwise estimator because it accelerates solver convergence (see Section 3 and it enables efficient tracking of predictive performance during optimisation. See Appendix for details.
Figure 10 visualises the evolution of the relative residual norm of the probe vector linear systems and the predictive test log-likelihood during marginal likelihood optimisation. A full of results is in Appendix
C For all solvers, warm starting leads to lower residual norms throughout outer-loop steps. This suggests a synergistic behaviour between early stopping and warm starting: the latter allows solver progress to accumulate across marginal likelihood steps. This can be interpreted as amortising the inner-loop linear system solve over multiple outer-loop steps. Despite the lower residual norm, CG is brittle under early stopping, obtaining significantly worse performance than AP and SGD on BU77 and HOUSFFLECTRIC AP and SGD seem to be more robust to early stopping. However, lower residual norms do not always translate to improved predictive performance. Furthermore, we find that SGD can suffer due the optimal learning rate changing as the hyperparameters change.
6 Conclusion
Building upon a hierarchical view of marginal likelihood optimisation, this paper consolidates several iterative GP techniques into a common framework, analysing them and showing their applicability across different linear system solvers. Overall, these provide speed-ups of up to 72x when solving until a specified tolerance is reached, and decrease the average relative residual norm by up to 7x under a limited compute budget. Additionally, our analyses lead to the following findings: Firstly, the pathwise gradient estimater accelcrates linear system solvers by moving solutions closer to the origin, and also provides amortised predictions as an added benefit by turning probe vectors into posterior samples via pathwise conditioning. Secondly, warm starting solvers at previous solutions during marginal likelihood optimisation reduces the number of solver iterations to tolerance at the cost of introducing negligible bias into the optimisation trajectory. Furthermore, warm starting combines well with pathwise gradient estimation. Finally, stopping linear system solvers after exhausting limited compute budget generally increases the relative residual norm. However, when paired with warm starting, solver progress accumulates, amortising inner-loop linear system solves over nultiple outer-loop steps. Nonetheless, we observe that low relative residual norms are not always necessary to obtain good predictive performance, which presents an interesting avenue for future research.

---

Acknowledgments
Jihao Andreas Lin and Shreyas Padhy were supported by the University of Cambridge Harding Distinguished Postgraduate Scholars Programme. José Miguel Hemándezz-ooato ledgess support from a Turing Al Fellowship under grant EP/V023756/1. We thank Runa Eschenhagen for helpful discussions. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.esd3.cg www.csd3.cam.ac.uuk ca).uukkkk provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/T022159) and DiRAC funding from the Science and Technology Facilities Council www.dirac.cd ak..). This work was also supported with Cloud TPUs from Google's 's TPU Research Cloud (TRC).
References
[1] J. Antorán, R. Barbano, Leuschner, J. M. Hernández-Loobtto, and B. Jin. Uncertainty Estimation for Computed Tomography with a Linearised Deep Image Prior. arXiv:2203.00479
2022. Cited on pages 6
[2] J. Antorán, S. Padhy, R. Barbano, E. T. Nalisnick, D. Janz, and M. HerándeezLobard
Sampling- based inference for large linear models, with application to linearised Laplace. In International Conference on Learning Representations, 2023. Cited on pages 4
[3] A. Artemev, D. R. Burt, and M. van der Wilk. Tighter Bounds on the Log Marginal Likelihood of Gaussian Process Regression using Conjugate Gradients. In International Conference on Machine Learning. 2021. Cited on pages 3
[4] J. Bradbury, R. Frostig. P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman- Minnl,, and Q. Zhang. JAX: composable transformations of Python+NumPy programs, 2018. Cited Page[7
[5] H. Chen, L. Zheng, R. Al Kontar, and G. Raskutti. Stochastic Gradient Descent in Correlated Settings: A Study on Gaussian Processes. In Advances in Neural Information Processing
Systems, 2020. Cited on page [
[6] M. P. Deisenroth, D. Fox, and C. E. Rasmussen. Gaussian Processes for Data-Efficient Learning
in Robotics and Control. IEEE Trans. Puttern Anal. Mach. Intell., 2015. Cited on page
[7]
D. Dua and C. Graff. UCI Machine Learning Repository, 2017. Cited on pages 17
[8] E. N. Epperly, A. Tropp, and R. J. Webber. XTrace: Making the Most of Every Sample in Stochastic Trace Estimation. Matrix Analysis and Applications, 45(1), 2024. Cited on page 5
[9] J. Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson. GPyTorch: Blackbox Matrix-matrix Gaussian Process Inference with GPU Acceleration. In Advances in Neural
Information Processing Systems, 2018. Cited pages 1 3 6 g
18 19 35
[10] P. Ghasemi, M. Karbasi, A. 7. Nouri, M. S. Tabrizi, and H. M. Azamathulla. Application
of Gaussian process regression to forecast multi-step ahead SPEI drought index. Alexandria
Engineering Journal, 60(6), 2021. Cited on page [
[11] M. N. Gibbs and J. C. MacKay. Efficient Implementation of Gaussian Processes. Technical report, University of Cambridge, 1997. Cited on page [
[12] J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian Processes for Big Data. In Uncertainty in Artificial Intelligence, 2013. Cited on page [
[13] M. Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in Statistics Simulation and Computation, 19(2). 1990. Cited on page
[14] D. P. Kingma and J. Ba. Adam: A Method for Stochastic Optimization. In International Conference on Learning Representations, 2015. Cited on pages 18
[15] J. A. Lin, J. Antorán, S. Padhy, D. Janz, J. M. Hermández-Lobato, and A. Terenin. Sampling
from Gaussian Process Posteriors using Stochastic Gradient Descent. In Advances in Neural
Information Processing Systems, 2023. Cited on pages I 3 9 18 20 35
[16) J. A. Lin, S. Padhy, J. Antorán, A. Tripp, A. Terenin, C. Szepesvári, J. M. Hernández-LLobato, and D. Janz. Stochastic Gradient Descent for Gaussian Processes Done Right. International Conference on Learning Representations, 2024. Cited on pages
18
20
35
10

---

[17] J. A. Lin, S. Padhy, B. Mlodozeniec, and J. M. Hemándezz-obato Warm Start Marginal Like- lihood Optimisation for Iterative Gaussian Processes. In Advances in Approximate Bayesian Inference, 2024. Cited on page 6
[18] W. J. Maddox, S. Kapoor, and A. G. Wilson. When Iterative Gaussian Processes Reliably
Accurate? In ICML OPTML Workshop, 2021. Cited on pages 3
[19] R. Meyer, C. Musco, C. Musco, and D. Woodruff. Hutch++: Optimal Stochastic Trace Estima- tion. Symposium on Simplicity in Algorithms, 2021, 2021. Cited on page
[20] J. Quifionro--Cadern and C. E. Rasmussen, A Unifying View of Sparse Approximate Gaussian
Process Regression. Journal of Machine Learning Research, 6. 1, 2005. Cited page
[21] A. Rahimi and B. Recht. Random Features for Large-scale Kernel Machines. In Advances in
Neural Information Processing Systems, 2008. Cited on pages B 18
[22] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press,
2006. Cited on page I
[23] S. Shalev-Shwartz and TT. Zhang. Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization. Journal of Machine Learning Research, 14, 2012. Cited pages I 3
[24] J. Snoek, H. Larochelle, and R. P. Adams. Practical Bayesian Optimization of Machine Learning Algorithms. In Advances in Neural Information Processing Systems, 2012. Cited on page
[25] D. J. Sutherland and J. Schneider. On the error of random Fourier features. Uncertainty in Artificial Intelligence, 2015. Cited on pages 18
[26] K. Tazi, J. A. Lin, R. Viljoen, A. Gardner, T. John, H. Ge. and R. E. Turner. Beyond Intuition, Framework for Applying GPs to Real- World Data. In ICML Structured Probabilistic Inference
& Generative Modeling Workshop, 2023. Cited on page [
[27] M. K. Titsias. Variational learning of inducing variables sparse Gaussian processes. In
Artificial Intelligence and Statistics, 2009. Cited on page
[28] S. Tu, R. Roelofs, S. Venkataraman, and B. Recht. Large Scale Kernel Learning using Block
Coordinate Descent, arXiv: 1602.05310, 2016. Cited on pages 3
[29] R. Vershynin. Highh-iiensional Probability: An Introduction with Applications in Data Science. Cambridge University Press, 2018. Cited on pages 14
[30] K. A. Wang, G. Pleiss, R. Gardner, S. Tyree, K. Q. Weinberger, and G. Wilson. Exact Gaussian Processes on a Million Data Points. In Advances in Neural Information Processing
Systems, 2019. Cited on pages [
19 35
[31] J. T. Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Efficiently Sampling Functions from Gaussian Process Posteriors. In International Conference on Machine Learning, 2020. Cited on pages
18
[32] J. T. Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Pathwise Conditioning of Gaussian Processes. Journal of Machine Learning Research, 22, 1, 2021. Cited on pages
[33] K. Wu, J. Wenger, H. Jones, G. Pleiss, and J. R. Gardner. Large-Scale Gaussian Processes via Alternating Projection. In International Conference on Artificial Intelligence and Statistics, 2024. Cited pages]
19 35

---

A Mathematical Derivations
In this appendix, we provide mathematical derivations for claims in the main paper.
A.1 Variance of Standard and Pathwise Gradient Estimator
To compare the variances of the standard estimator (6) and the pathwise estimator (
we calculate
HHo
with (O. and Var [x eTHHo
Var
z H
&
with ~ N(". (17)
a
08g
The variance of the standard estimator is given by
HHg
HHe
HHg
HHg
Var
= H
= ir
H
H
H
(18)
00k-
00;
00k
00,
HHg
HHg
HHg
HHg
tr
0
202
H
0
d0n
+ tr
H
80
00
H
(19)
The variance of the pathwise estimator is given by
HH
(HH
HHg
Var
#
2
2r
H
H
(20)
00p
0
002
HHg
HHo
HHe
HHg
H
H
H
H
(21)
00
00
00g
00,-
Therefore, the variances of both estimators share the first trace term and only differ in the second trace term. Hence, their variances will be identical if
Lr H HHe HHg
1 Ir H
HH.
H
HHe 00g-
(22)
00
00
a0,
which is the case if H and OH,/00, commute with each other. For example, consider the derivative with respect to the noise scale
HH. a do do
(k(2, x; 1 + I 2I1.
(23)
In this case, H and 0Hg/80 'o commute with each other, such that both estimators have the same variance. In general, sufficient condition for matrix multiplication be commutative is simultaneous diagonalisability of two matrices.
A.2 Taylor Approximation View of Warm Start
At iterations / and ( 1 of the outer-loop marginal likelihood optimiser, associated with (()) and g(++1) the linear system solver must solve two batches of linear systems, namely
H
(t)
y. (4)
(a)
and
(24)
H
((+1)
(1+1)
1j (++1)
((+1)
y 21 (++1)
(1+1)
(25)
where ( and ((-")) are related through the change from (()) to (++1) and vy
(t) and vy
(I+1)
further related through sharing the same right-hand side y in the linear system. In such setting, where the coefficient matrix only changes slightly and the right-hand side remains fixed, we can approximate ((++) using a first-order Taylor expansion (((1))
H (1((t)) -1
(H(() (1(5))' (H!
H (l+1)
H²0) (H!)) (H!')²
(26)
(++1) s (t) (H!))
H
H(!)) (1))
(27)
H is small then (t) will be close to (1+1), such that we can reuse (1) to initialise the linear system solver when solving for (f++1) To satisfy the condition of fixed right-hand sides, we must set ( 2j at the cost of introducing some bias throughout optimisation, which turns out to be negligible in practice (see Section 4 and Appendices A.3 and A.4 for details).
If A = H (++1)
12

---

A.3
Convergence of Warm Starting Marginal Likelihood Optimisation
Recall the gradient of the marginal likelihood objective:
aC(0) 1 (Hg'y)) 0
HH Hg' y 00
1 br H 2
HHo 00
k E (l,....,o).
where Ho E RR is positive semi-definite symmetric matrix, E I a real vector, n is the number of "data" examples, and de is the number of hyperparameters. Also, recall the the warm start estimator Ou (0) to the gradient /000)
OHg
On (0)
2 H
Zj.
00k
==1
where the probe vectors 2j random variables with identity second moments: [[zzz]] I, and is the number of probe vectors in the trace estimator.
We will write Sn-1 {x E R" ||ح|| 1} for sphere in R" For real matrix A E IR mx we will denote the operator (spectral) norm supæesa-L yupyesm-1 y Ax with |
Notation
Definition 2 (Sub-gand Suub-gaussian norm). The sub-gaussian norm of sub- gaussian random variable X is defined as:
X 11
inf 0: D: 2)
Definition 3 (Sub-exponential norm). The sub-exponential norm of a sub-exponential random variable X is defined as:
X lv = inf(t 0: [[ 2}
We will denote the optimisation domain for the hyperparameters θ where we assume C IR
Theorem Assume that the probe vectors (Zi....; 2.) are zero meun, coordinate-wise independent, and that elements of 2j are sub-gaussian with norm وz = o'v {I,.... nn.j (I,..., j (1...., 8).
Assume that the of the singular values of H d0,
HHe for all k {1,..., do upper-bounded
on the domain of 0 by max Then, for all 8 0:
max (V log
P
OC g(0)
00)
9dg 25
n C1 log
9dg 25
C2or X maa
1 1-65 8,
00
The crux of the proof of Theorem comes from bounding the spectral norm of the difference (Ej=12 - To do so, it is useful to introduce the following definitions and lemmas.
-1)
Lemma 5 (Computing the operator norm on a net [29 Exercise 4.4.3)). Let A be an n x matrix and e [0, 1). Then, for any e-net E,, of the unil sphere 1 we have:
sup x S |p S .EΣ,
sup x Ax. #<Σ.
1 - 25
Lemma 6 (Size of E-net on - 29 Corollary 4.2.13]). There exists an E-net on S with cardinality at most (2 1))
Lemma 7. Let M E
be an n random matris, where zÿ are independent and identically distributed sub- sabb-gassiann random variables with sub- gaussian norm ||=illyy a, and Sn-1 any unit vector in Raa Then:
C2
P 2 2 3 2e
Cgo
where C1 C2 absolute constants.
13

---

Proof. We can rewrite:
2
x Ma x
Σ [2)2]
1)
Σ
2x)
x.
Fizej
-1
3=1
Zijij weighted sum of independent sub-gaussian random variables, it is also sub- gaussian with squared norm:
Since E
2ij
C
2 I = C
1 = Co?
==
where C' is an absolute constant 29 Proposition 2.6.1]. I. Hence, since (E z;2iy) ) is the square of sub-gaussing gagssiasss random variable, must be sub-exponential with the sub-exponential norm (see [29
Lemma 2.7.6]):
Σ
Tizij
I
S CCo,
Vi
Lastly, since EE T; is sub-exponential, (E miza))
so
will the mean-centered counterpart
E
(2 1= )3 )
(E miziy) )²
29 Exercise 2.7.10) with sub-exponential
norm:
-1
S C20
where C2 is another absolute constant. Hence, since (Ei= i= x4248) sub-exponential, zero-mean and independent, we can apply Bernstein's 's inequality bound the tail probability of x: Ma:
for j E {.....,s} are
Pl =
Σ IzZt
2
1
>
2e
C2 min
.2 Cg.
Ezo
holds for any B2 O, where C1, C2 are absolute constants.
Lemma 8. Ler M = zzz] I be an n n random matrix, where 2 are independent and identically distributed sub- Gaussian random variables with aab-gaussian norm || a. Then;
|M 2 max
(V log
(
9
n C1 log
9 28
C4
<
Vg 0,
28
where C1, C4 are absolute constants.
Proof. Pick an s-net Σg on S of size at most (3 + 1)" (Lemma 6).
Then, we can bound the tail
probability of the operator norm as:
P | 2
8
< P
sup
TM 2 (1- - 25)
A By Lemma
LEE
P Mc 2 (1 - 2=)
A Union bound
me
A p(1-20)
< Σ 2e
Cc min
clo²
A By Lemma
2 + +1)
-C min
(1 Bl1-20)
<
2e
oz
8
Setting g gives:
M 2/8
5 2e
-C min से. وه
2
14

---

where C1 is an absolute constant.
Lastly, to get the bound into the form P
M op
2 f(o, n, s)
S 6, we note that:
C2 min
8
B
9
2e
c40
2
min
C1 log
0f
C40
28
max
"
Ci log
C1 log
9
Cạ
28
26
and so by setting 3 = max (V log
(2))
_ C1 log
(3))
C4a we get the desired result:
P M 2 max
" log
n C1 log
9
C40
S 6
25
26
Proof of Theurem Let Eisi pp be a singular value decomposition (SVD) of def
H
HH
a0,
where {qs ) and (p:) are two sets of orthonormal vectors. First, note that we can rewrite:
OC(0)
2 Az
Ez x=Az)
00,
=
j=1
Σ x]
Aaqipk
2j -
Ez
xT Σ Aappp
(28)
j=1
Σ
A
2) q:p! 2j
Σ
Ne Ez 1 s.q.p.a. x]
j=1
i=1
A
2 P.P. Zj
Ez |
c.P. 2|
j=1
Σ
N
Σ
Zjz
T
E2 zz
p;
1
j=1
pbx
Σ 2
Di.
(29)
M
Therefore, can bound the norm of the difference as
Th (0) 9% (0)
< Adl l Mp:
(30)
| max [ op
(31)
By Lemma 8 with probability at least 1 we can bound the operator nomn of M as:
Thx (0) gz (0) < X max [M] < max (V2C. log
9 n C1 log
C20 max
25
25
with C1, C2 absolute constants.
15

---

We can apply a union bound over all & E {1....,,d) do to bound the probability of the &o -oom of the gradient deviating by a certain amount:
ac 00)
max (v c log
9
9
I
g(0)
C1 log
C20 max
> 1 dØd,
00
28
26
or:
Oc g(0) 00
(0)
max
V Ci log
9dg
n
C1 log
9dg
)
C20
>1 d.
28
25
Now, if (0) is conservative field, and so is implicitly gradient of some (approximate) objective
E: ④ + IR. the above result allows to bound the enror the solution found when optimising using the approximate gradient 9 instead of the actual gradient g = V.. However, in general, g()) need not be strictly conservative. In practice, since (0) converges to conservative field the more samples we take, we may assume that it is close enough to being conservative for the purposes of optimisation hardware with finite numerical precision. Assuming that g(0) is conservative allows us to show the following bound on the optimum found when optimising using g(0), which is a restatement of T
Theorem 9. Let 9 and C be defined as in Theorem 4 Assume the optimisation domain O is convex, closed and bounded. Then, with probubility at least 1 8:
Assume ğ O - R is a conservative field.
((0')> max (1
" log
9do 25
"dg 26
C1 log
C20 max /
where AO def supo.ore ||0' - 0|| is the maximum distance between two elements in e.
Proof. Let L: θ - RR be approximate objective implied by the gradient field g, namely scalar field such that VL g. Such scalar field exists if is conservative field, and is unique up to constant (which does not affect the optimum).
For any points 0, O, with A0 def 0' 0, have that
(C(0') - L(0)) (((0') [(0)
A Replace difference in values with integral along path from 0'
/'' Ōt L (0 + AOt) dt [-" Ot
[ Alt) dt
AO VC + AGt) di f A0 vā A0t) dt
AO (VC AOt) vā stt)) dt
< AO (VC(0+ 200) V4 + 201) dt
[ ||0||| | (VL 200) vÂ A01) )
|| ( (VL (0 + 201) - vē + 201)
dddt
A Bound fo-m by the loc norm
< || |A0 max V
C₁ log
9dg
a C₁ log
9do
))
C2aA max dodt
28
s
25
16

---

Difference of gradients bounded with probability at least (1 - 6) by Theorem
Th (v
9dg
9dg
S max
log
C1 log
C2a max A0 do.
28
28
The above inequality holds with probability at least (1 6). Hence,
Negative because 0" is maximum 2
((0') - 2(0)) ((0)) ((O') (E(07) 2(00))
$ [(0') C(0") - (Z((*) [(0))
max (12Ci log
9da
9dg
- -i log
C20r max AOV do.
20
26
Remark 10. Theorem 9 above implies that the objective of the optimum 0" obtained with the approximate gradients converges to the objective of the Irue optimum 0" in probability:
Ya > 0: P - > a) + +0
0c
(32)
which trivially follows from the implication of Theorem that for every a, 8 O. we can find an s E B such that L(O')= 2 C(0") - a with probability at least 1 1 8.
Remark 11 (Convexity of θ). We also note that the convexity of the hyperparameter domain O is fairly mild assumption which is satisfied in the majority of practical settings. For example, optimising kernel length scales and noise scale bounded intervals (10-19, 10 1100) falls within the assumptions, but introducing a "hole" " into the domain (e.g. introducing a constraint like || - 1| 2 0.5) would break the assumption. In particular, this assumption is NOT a statement about the convexity of the objective C- in the proof. we allow for the objective to be arbitrarily non-convex, and only assume its differentiability.
A.4 Convergence of Warm Starting the Pathwise Estimator
The result in Appendix A.3 can be trivially extended for the pathwise estimator in P for any pairwise independent probe vectors Ê, (with second moment Hg²) that upon rescaling by H will be zero-mean with independent coordinates. This is true for probe vectors 2, that are either i.d.d.
(O, Hp'
"1'')ddisribbute )-distributed or obtained by transforming Radamacher random variables by H
B Implementation and Experiment Details
In this appendix, provide details about our implementation and experiments.
General Our implementation uses the JAX library 4 All reported experiments were conducted on internal NVIDIA A A100SXXM4-80GB GPUs using double floating point precision. Some additional experiments and ablations were performed on Google Cloud TPUs (v4). The total compute time, including preliminary and failed experiments, and evaluation is around 4500 hours. The compute The source code is available here
time of individual runs is reported in Tables 2 to I0
Datasets Our experiments are conducted using the datasets and data splits from the popular UCI regression benchmark [7]. They consist of various high-dimensionas multivariate regression tasks and are available under the Creative Commons Attribution 4.0 International (CC BY 4.0) license. In particular, we used the POL (n = 13500, d 26), ELEVATORS (n 14940, d = 18), BIKE (n 15642, d = 17), PROTEIN (n. = 41157, d = 9), KEGGDIRECTED (n 43945, d = 20). 3DROAD (n = 391387, d = 3), SONG
463811, d = 90), BUZZ (n 5524925, d = 77), and
HOUSEELECTRIC (n 1844352, d = 11) datasets.
17

---

Kernel Function and Random Features In all experiments. we used the Matérn-Y2 kernel, param- eterised by scalar signal scale and length scale per input dimension. For pathwise conditioning (3) and the pathwise gradient estimator (9). we used random Fourier features 21 25 (1000 sin/cos pairs, 2000 features in total) to draw approximate samples from the Gaussian process prior. For an explanation about how to efficiently sample prior functions from a Gaussian process using random features, we refer to existing literature [31] 32, 15]. However, we want to discuss some details in terms of using this technique for the pathwise estimator from Section B
For pathwise gradient estimation, the linear system solver must solve linear systems of the form
Ha (yy & ,,.....,,..))) 1 EI,...,,..]]
(33)
with ş = f(xr) e, where [(x) ~ (V, k(x, æ;; 2) is a prior function I sample evaluated at the training data at and E <~ _V( is Gaussian random vector. Both quantities are resampled in each outer-loop marginal likelihood step if the pathwise estimator is used without warm starting. With warm starting enabled, the right-hand sides of the linear system must not be resampled. In this case, f and are sampled once and fixed afterwards. However, f depends on & and depends on ø, and both & and are hyperparameters which change in each outer-loop step. Therefore, what does mean to keep f and fixed?
For c, this amounts to the reparameterisation fixed afterwards, such that e becomes deterministic. For f. this refers to fixing the parameters of the random features, for example the frequencies in the case of random Fourier features. Intuitively, this corresponds to selecting a particular instance of prior sample, although the distribution of the sample can change due changes in the hyperparameters. In each outer-loop step, the random features evaluated using the fixed random feature parameters and the updated kernel hyperparameters, and the prior function sample is then evaluated the training data using the updated random features. Both of these operations are (n) and efficient as long as the number of random features is reasonable.
w. where ~ (VOO.)) I) is sampled once and
Iterative Optimiser To optimise hyperparameters @ given estimate of V., used the Adam optimiser (14) with default settings except for the learning rate. For all small datasets (n < 50k), initialised the hyperparameters at 1.0 and used a learning rate of 0.1 perform 100 steps of Adam. For all large datasets (n > 50k), we initialised the hyperparameters using a heuristic and used a learning rate of 0.03 perform 30 steps of Adam (15 for HOusEELECTRIC due to high computational costs). The heuristic to obtain initial hyperparameters for the large datasets consists of:
1. Select centroid data example uniformly at random from the training data.
2. Find the 10k data examples with the smallest Euclidean distance to the centroid.
3. Obtain hyperparameters by maximising the exact marginal likelihood using this subset.
4. Repeat the procedure with 10 different centroids and average the hyperparameters.
This heuristic has previously been used by Lin et al. [15
[6 to avoid aliasing bias.
To enforce positive value constraints during hyperparameter optimisation, we used the softplus function. In particular, we reparameterise each hyperparameter 0₂ R> as 0 = log(1 + exp(v*)) and apply optimiser steps to 142 IR instead, to facilitate unconstrained optimisation.
Gradient Estimator For all experiments, unless otherwise specified, the number of probe vectors was set to s = 64 for both the standard and the pathwise estimator. The distributions of the probe vectors are z ~ NN I) for the standard estimator and & := Hg 5 ~ NO H for the pathwise estimator. See Section for details about how to generate samples from (0, H)²).
The probe vectors used by Gardner et al. (9 have conceptual similarities but the motivation is different. They used probe vectors = ~ (O, P). where P is constructed using a low-rank pivoted Cholesky decomposition to implement the preconditioner. In contrast, our pathwise probe vectors are sampled using random features, and, for CG, we used the pivoted Cholesky preconditioner in addition to the pathwise probe vectors. Furthermore, used the solution of the pathwise probe vector systems to construct posterior samples via pathwise conditioning, which has not been done by Gardner et al. (9).
The name of the puthwise estimator can be related to the reparameterisation trick by viewing a sample from the GP posterior as a deterministic transformation of sample from the GP prior, and observing that the latter itself is an affine transformation of standard normal random variable.
18

---

Linear System Solver We conducted two sets of experiments which only differ in the termination criterion of the linear system solver. In the first set, we stop linear system solvers once they reach a relative residual norm tolerance of 0.01. In the second set, we also restrict the maximum number of solver epochs to 10, 20, 30, 40 or 50, such that most of the time the residual norm does not reach т.
For a generic system of linear equations Hg u = b, the residual is defined as r = b - Hg u and the relative residual norm is defined as ||r||/||||| In practice, improve numerical stability, the relative residual norm tolerance is implemented by solving the system Hg u b. where b = b/ (l]bb] + e), until ||=|| := b Hg أ||l S T and then returning u = (|]bi| + c) čž, where epsilon is set to a small constant value prevent division by zero. Since we solving batches of systems of linear equations of the form Hg Uyy vi....., Us = y, 21..... 2, ]. we track the residuals of each individual system and calculate separate residual norms for the mean and for the probe vectors, where the residual norm for the mean ||r,y || corresponds to the system Hg By y and the residual norm for the probe vectors ||=|| is defined the arithmetic average over residual norms corresponding the systems He [Mi...., V, == [>1,-.., ]. Both relative residual norms must reach the tolerance 7 to satisfy the termination criterion. We use separate residual norms because ||r typically converges faster than ||=]], such that an average other all systems tends to dilute the latter (see Figures 14 to 17
Conjugate Gradients The conjugate gradients algorithm 9 30 computes necessary residuals as part of the algorithm. In terms of counting solver epochs, every conjugate gradient iteration counts as one solver epoch because in every iteration each value of Hg is computed once. Following previous work, we used a pivoted Cholesky preconditioner of rank 100 for all experiments [30]. We initialised conjugate gradients either at zero (no warm start) or the previous solution (warm start). Otherwise, conjugate gradients does not have any other parameters. Pseudocode is provided in Algorithm [
Algorithm 1
Conjugate gradients for solving He (Vy. .1.... Us. y, 21...., 3.
Require: Linear operator (.),), targets [y, 21, Fi,..., #s tolerance т, maximum epochs T
Require:
Preconditioner P..)
1: Let (.), denote parallel execution over (.)y (),...., (-)
2: < 0 (or previous solution if warm start)
3: T. + b. (H(
4: p. + P(r.)
5: d. + Pz
6: % rr.p.
7: t+0
8: while
< T and ||r, > r and Ei
||rl|| = |rz|| > do
9:
a. s 7./d. Ho (d.)
10:
u. u. + a.dd
11:
T. + r. - a. (.H)
12:
P. P(r.)
13:
B. r_p-//
14:
7. + r.P. p.
d, r p. + 3.d.
16:
1+1+1
17:
end while
18: return By. 01 035
Alternating Projections The alternating projections algorithm 33 also keeps track of the residuals as part of the algorithm. In terms of counting solver epochs, we convert the number of maximum solver epochs
maximum number of solver iterations by multiplying with n/b, where 6 is the block size, because every iteration of alternating projections computes b/ of all entries of Hg. We used a block size of b 1000 for all datasets, except PROTEIN and KEGGDIRECTED, where we used b = 2000 instead. We initialised alternating projections either at zero (no warm start) or at the previous solution (warm start). During each marginal likelihood step, the Cholesky factorisation of every block is computed once and cached afterwards (although, in practice, the Cholesky factorisation does not dominate the computational costs). In each iteration of alternating projections, the block with largest residual norm is selected to be processed. Pseudocode is provided in Algorithm 2
19

---

Algorithm 2 Alternating projections for solving Hg [Vy: 01. y, ,,,...,
Require: Linear operator H:):), targets
3/ #1...... 2, tolerance 7, maximum epochs T
Require: Block size b, block partitions [1], [2]..... [[8]]
Let (-), denote parallel execution over (.) (.),......
2: < 0 (or previous solution if warm start)
3: r. + b. Ho(v.)
4: t+0
5: while 7 and |ryl and - by |r||| > T do
6:
[]] f arg arg_max([Irg
E j= 7, 1,
| Ej= ?
7:
v.[i] v. + chhol_solee_____ Ho i i...).)))
8:
T. + r. (chol_ Hp
9:
1+- 1+1
10:
end while
11: return By. U1 U/
16 does not compute Stochastic Gradient Descent The stochastic gradient descent algorithm 15 residuals as part of the algorithm. Therefore, we estimate the current residual by keeping a residual vector in memory and updating it sparsely whenever leveraging the property that the negative gradient is equal to the residual. practice, we find that this estimates approximate upper bound the true residual, which becomes fairly accurate after few iterations. In terms of counting solver epochs, we apply the same procedure as for alternating projections. The number of maximum solver epochs is converted to maximum number of solver iterations by multiplying with n//b, where b is the batch size. We used a batch size of 500, momentum of p 0.9, and Polyak averaging, because averaging would interfere with our residual estimation heuristic. We use learning rates of 30, 20, 30, 20, and 20 respectively for the POL, ELEVATORS, BIKE, KEGGDIRECTED and PROTEIN datasets, picking the largest learning rate from grid [5, 10, 20, 30, 50, 60. 70, 80, 90, 100] that does not cause the inner linear system solver diverge on the very first outer marginal likelihood loop. For the larger datasets, we use learning rates of 10, 10, 50, and 50 for 3DROAD, BUZZ, SONG and HOUSEELECTRIC, picking half of the largest learning rate above. We find that the larger datasets are more sensitive to diverging when the hyperparameters change, and therefore we choose half of the largest learning rate possible at initialisation. Pseudocode is provided Algorithm
compute the gradient on batch of data,
strictly necessary (16) and
Algorithm 3 Stochastic gradient descent for solving Ho [uy U1...., U. y. 21.....
Require: Linear operator (.).). targets = [y, 21,...,, 2,),, tolerance т, maximum epochs T Require: Batch size b, learning rate " momentum, p
1: Let (-), denote parallel execution over (-)y, (.).....,() (-),
2: v. < 0 (or previous solution warm start)
3: r. + b.
4: m. +0 0
5: t+0
6: while
T and ||rl|| and ! E |ry||
|r>/| > 7 do
7:
[ + uniform_he wwtp__aahhh__aal(((((e
8:
g. 0
9:
[ + Hol :l (v.) b
10:
m. + pm. 9-
11:
v. = =V,- m.
12:
1 [
13:
1 1+111
14:
end while
15:
return By. U1, 0,,..,,, 00
20

---

C Additional Empirical Results
In this appendix, we provide additional result from our experiments.
For our first experiment (solving until reaching the relative residual norm tolerance), we present the predictive performance, time taken, and speed-ups in Tables to Hyperparameter trajectories illustrated in Figures to 13 Required number of solver iterations are shown in Figure 21
For our second experiment (solving until reaching the tolerance exhausting the compute budget), we tabulate the predictive performance, time taken, and average residual norms in Tables TO The behaviour of residual norms is visualised in Figures to 17 The evolution of residual norms and predictive performance via pathwise conditioning is depicted in Figures 18 20
Table 2: Results on POL when solving until convergence (mean + standard error over 10 splits).
palh
warm
POL (n
13 500, d 26)
wise
start
Test RMSE
LLH
Total Time (min)
Solver Time (min)
Speed-Up
0.0750 0.0010
1.2682 0.0084
4.8263 + 0.0356
4.6138 0.0355
CG
0.0754
0.0010
1.2716 = 0.0077
3,9567 = 0.0204
3,7148 = 0.0205
1.2 x
0.0750 0.0010
1.2681 + 0.0084
2.2844 + 0.0131
2.0755 + 0.0131
2.1 x
0.0758 0.0010
1.2666 + 0.0074
2.4652 + 0.0133
2.2513 + 0.0133
2.0 x
0,0750 0.0010
1.2682 + 0.0084
493.05 + 3.0738
492.84 # 3.0737
AP
0.0754 0.0010
1.2715 + 0.0077
44.006 + 0.0722
43.765 = 0.0725
11.2 x
0.0750 0.0010
1.2681 + 0.0084
27.915 0.2707
27.706 0,2708
17.7 x
0.0758 0.0010
1.2666 + 0.0074
3.8962 + 0.0174
3.6831 = 0.0173
126.6 x
0.0750 0.0010
1.2681
0.0084
138.63 0.8169
138.60 0.8170
SGD
0.0754 0.0010
1.2708 1 0.0074
73.546 1 0.2661
73.434 1 0.2661
1.9 X
0,0750 0.0010
1.2682 1 0.0084
26.483 1 0.1002
26.461 1 0.1001
S.2 X
0.0757 0.0010
1.2678
= 0.0076
17.938 = 0.1083
17.854 = 0.1083
7.7 x
Table 3: Results on ELEV when solving until convergence (mean + standard error over 10 splits).
path
ELEVATORS (n
14 940, d = 18)
start
Test RMSE
Test LLII
Total Time (min)
Solver Time (min)
Speed-Up
0.3550 = 0.0034
-0.3856 = 0.0065
1.5811 + 0.0063
1.3661 = 0.0062
CG
0,3562 = 0.0033
-0.3868
0.0065
1.4940
0.0039
1.2484 0.0039
1.1
0.3550 + 0.0034
-0.3856
0.0065
1.0308
0.0037
0.8191 0.0036
1.5
0.3558 + 0.0034
-0.3856 + 0.0066
0.9977 + 0.0029
0.7835 + 0.0028
1.6
0,3550 = 0.0034
-0.3856 = 0.0065
77.787 = 0.2694
77.572 = 0.2693
0.3562 = 0,0033
-0.3668 = 0.0065
36.346 + 0.0821
36.100 = 0.0820
2.1
AP
0.3550 - 0.0034
-0.3856 + 0.0065
1.6658 + 0.0081
1,4541 + 0.0081
46.7
0.3558 + 0.0034
0.0066 -0.3856
0.0032 1.2065
0.9928 0.0031
64.5
0.3550 = 0.0034
-0.3855
0.0066
5.5408
0.0101
5.4091 = 0.0101
SGD
0.3562 = 0.0033
-0.3668
0.0065
4.5758
0.0086
4.4146 0.0086
1.2
0.3550 - 0.0034
-0.3855 1 0.0066
1,2221
1 0.0042
1.0927 1 0.0041
4.5
0.3558 = 0.0034
-0.3554 + 0.0068
1.1437 + 0.0033
1.0105 - 0.0032
4.8 x
21

---

Table 4: Results on BIKE when solving until convergence (mean + standard error over 10 splits).
path
BIKE (n
15 642, d = 17)
Test RMSE
Test LLH
Total Time (min)
Solver Time (min)
Speed-Up
0.0326 0.0031
2.1500 + 0.0180
5.0797 + 0.0280
4.8412 = 0.0280
CG
0.0326 0.0030
2.0674 0.0167
4.4147 = 0.0225
4.1410 0.0221
1.2
0.0327 = 0.0031
2.1508 = 0.0181
2.7392 = 0.0137
2,5039 = 0.0135
1.9 x
0.0329 0.0030
2.0615 + 0.0151
3.0658 1 0.0168
2.8268 1 0.0167
1.7 x
0.0326 0.0031
2.1504 = 0.0180
302.26 + 1.7735
302.03 = 1.7735
0,0325 0.0030
2.0668 = 0.0167
35.081 = 0.0773
34,811 = 0.0772
8.6 x
AP
0.0326 0.0031
2.1503 = 0.0181
19.892 = 0.2146
19.657 = 0.2146
15.2 x
0.0330 0.0030
2.0616 + 0.0150
5.4041 + 0.0394
5.1653 + 0.0393
56.0 X
0.0326 0.0031
2.1535 = 0.0181
412.17 10.460
412.03 = 10.450
SGD
0.0324 0.0030
2.0692 = 0.0174
156.24 2.2113
156.05 2.2109
2.6 x
0.0327 0.0031
2.1524 + 0.0179
74.341 = 1.2532
74.184 = 1.2524
5.5 x
0.0332 0.0030
2.0562 + 0.0144
64.145 + 1.1086
63.983 + 1.1083
6.4 x
Table 5: Results on PROT when solving until convergence (mean + standard error over 10 splits).
path
warm start
PROTEIN (n 41 157, 9)
wise
Test RMSE
Test L...
Total Time (min)
Solver Time (min)
Speed-Up
0.5024 + 0.0036
-0.5871 + 0.0096
29.849
0.2463
28.208 + 0,2463
0.4909 = 0.0035
-0.6210 = 0.0082
19.984 + 0.1429
18.317 = 0.1430
1.5
0.5026 = 0.0036
-0.5871 + 0.0096
11.542 0.1058
9.8816 = 0.1059
2.6
0.4912 = 0.0034
-0.6214
0.0079
0.0959 13.744
12.085 = 0.0956
2.2
0.5024 = 0.0036
-0.5871
0.0096
130.93
1.0450
129.29 1.0451
SGD AP
0.4907 1 0.0035
-0.6214
0.0082
55.790
0.2054
54.125 0.2053
0.5027 + 0.0036
-0.5871 1 0.0097
16.425 1 0.1816
14.765 1 0.1818
8.0
0.4912 1 0.0034
-0.6213 1 0.0079
12.341 + 0.0347
10.682 + 0.0346
10.6
0,5026 = 0.0036
-0.5871 = 0.0096
75.205 + 1.4358
72.932 = 1.4357
0.4894 = 0.0035
-0.6268 = 0.0077
24.030 0.1364
21.726 = 0.1364
3.1
0.5027 0.0037
-0.5878
0.0096
11.226 0.1430
8.9541 0.1433
6.7
0.4911 + 0.0034
-0.6217
0.0078
11.939 0.0700
9.6617 + 0.0700
Table 6: Results on KEGG when solving until convergence (mean + standard error over 10 splits).
path
KEGGDIRECTED (n Test LLH
43 945, d 20)
wise
start
Test RMSE
Total Time (min)
Solver Time (min)
Speed-Up
0.0837 0.0016
1.0818 0.0170
27.974 + 0.3172
25.543 + 0.3120
0.0837 + 0.0016
1.0818 + 0.0169
26.362 + 0.2851
23.897 = 0,2804
1.1 x
0.0837 + 0.0016
1.0816 0.0171
12.754
+ 0.1314
10.326 + 0.1266
2.2 x
0.0836 + 0.0016
1.0819 + 0.0166
12.998 + 0.1383
10.559 + 0.1330
2.2 x
> 24 h
> 24 h
SGD AP
0.0837 + 0.0016
1.8200 + 0.0166
491.41 + 0.4624
488.95 +'t 0.4657
> 2,9
0.0837 0.0016
1,0818
0.0172
211.28 + 1,9504
208.85 + 1.9556
6.8
0.0836 + 0.0016
1.0817 - 0.0166
14.013 + 0.0768
11.574 + 0.0711
> 102.8 x
0.0837 0.0016
1.0816 116--1000 0.0173
620.07 1 6.3224
617.35 1 6.3194
0.0837 1 0.0016
1.0822
1 0.0164
411.58 1 6.1055
408.86 1 6.0874
1.5 x
0.0837 1 0.0016
1.8211
1 0.0170
168.38 1 1.6669
165.66 1 1.6636
3.7 x
0.0839 + 0.0016
1.0725
58.679 + 0.5457
55.952 + 0.5418
10.6 x
22

---

Conjugate Gradients
Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale
1.0
Selected Length Stales
2
0.5
2
POL,
0.0
0
0.5
3
3
3
0.0
Selocted Length Scales
0.5
3
2
ELEVATORS
0.0
0.5
2
0.0
Selected Length Scales
BIKE
0.0
0
0.0
1.0
Selected Length Scales
2
PROTEIN
0.0
1.5
2
0.5
1.0
Selected Length Scales
5
KEGGDIRECTED
1.0
5
5
50
100
100
50
100
50
100
50
100
Marginal Likelihood Steps
Exact
Standard
Pathwise
Standard warm start
Pathwise start.
Figure Evolution of hyperparameters during marginal likelihood optimisation on different datasets using conjugate gradients as linear system solver. Most of the time, the behaviour of exact gradient computation using Cholesky factorisation is resembled.
23

---

Alternating Projections
Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale
1.0
Selected Lenuth Scales
2
2
3
0.5
POL.
1.0
5
5
0.5
3
3
3
1.0
Selected Lenggh Scales
3
0.5
3
2
2
BLEVATORS
0.0
1.0
3
0.5
2
2
2
0.0
1.0
Selected Length Seales
10
10
0.5
2
5
5
BIKE
0.0
0
0
1.0
10
10
0.5
5
5
0.0
Selected Length Scales
0.5
PROTEIN
0.0
0
0
1.0
0.5
Selected Length Seales
5
KEGGDIRECTED
0.5
3
0.0
1.0
0.5
3
0.0
50
100
50
100
0
50
100
0
50
100
50
100
Marginal Likelihood Steps
Exact
Standard
Pathwise
Staudard
warm start
Pathwise + warm start
Figure 12: Evolution of hyperparameters during marginal likelihood optimisation different datasets using alternating projections as linear system solver. Most of the time, the behaviour of exact gradient computation using Cholesky factorisation is resembled.
24

---

Stochastic Gradient Descent
Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale Signal Scale Noise Scale
1.0
Selected Length Scales
2
2
3
0.5
POL.
1.0
5
5
0.5
3
3
3
1.0
Selected Lenggh Scales
3
0.5
3
2
2
BLEVATORS
0.0
1.0
3
0.5
2
2
2
0.0
1.0
Selected Length Seales
10
10
0.5
2
5
5
BIKE
0.0
0
0
1.0
10
10
0.5
5
5
0.0
Selected Length Scales
0.5
PROTEIN
0.0
0
0
1.0
0.5
Selected Length Seales
5
KEGGDIRECTED
0.5
3
0.0
1.0
0.5
3
0.0
50
100
50
100
0
50
100
0
50
100
50
100
Marginal Likelihood Steps
Exact
Standard
Pathwise
Staudard
warm start
Pathwise warm start
Figure 13: Evolution of hyperparameters during marginal likelihood optimisation different datasets using stochastic gradient descent as linear system solver. Most of the time, the behaviour of exact gradient computation using Cholesky factorisation is resembled.
25

---

Table 7: Results on 3DROAD with 10 maximum solver epochs (mean standard error over 10 splits).
3DROAD (n 391 387, d = 3)
Average Residual Norm of Probe Vectors
Test RMSE
Test LLH
Total Time (h)
of Mean
CG
1.2177 = 0.0296
-1.5463 = 0.0187
1.9308 = 0.0073
0,9815 = 0.0217
0,9174 + 0.0106 0.9003 0.0077
0.5371 0.0598
-0.9143 = 0.0336
1.9411 = 0.0047
0.6530 : 0.0152
AP
0.1042 0.0017
0.8237 + 0.0172
2.1782 = 0.0059
0.0950 + 0.0017
0.0541 + 0.0010
0.0563 0.0009
0.9309 + 0.0186
2.1805 = 0.0047
0.0469 + 0.0010
0.0651
0.0010
SGD
0.1430 0.0014
-0.4662 = 0.0399
1.2717 = 0.0014
0.1130 + 0.0017
0.0633
0.0020
0.0654 0.0010
0.9276 0.0126
1.2772 = 0.0031
0.0561 + 0.0010
0.0797 0.0022
Table 8: Results on SONG with 10 maximum solver epochs (mean + standard error over 10 splits).
SONG (n = 463 811, d 90)
Average Residual Norm of Probe Vectors
start
Test RMSE
Test LLH
Total Time (h)
of Mean
CG
2.1573 0.0672
-2.0688 + 0.0303
18. 18.107 + 0.0602
1.6528 + 0.0773
1.6356
0.0716
0.8698 0.0091
-1,3025 + 0.0140
17.879 + 0.1979
0.3793 0.0171
0.4147 0.0186
AP
0.7428 0.0019
-1.1197 = 0.0024
18.256 = 0.0272
0.0421
0.0016
0.0499
0.0011
0.7420 = 0.0019
-1.1184 + 0.0023
15.114 = 0.2373
0.0085 + 0.0004
0.0125 + 0.0002
SGD
0.7426 0.0019
-1.1205 + 0.0024
17.160 = 0.0669
0.0688
0.0017
0.0834
0.0014
0,7419 0.0019
-1.1184 + 0.0023
16.756 = 0.0842
0.0100 + 0.0003
0.0117 + 0.0001
Table 9: Results BUZZ with 10 maximum solver epochs (mean = standard error 10 splits).
warm start
B177 (n = 524 925, d 77)
Average Residual Norm of Probe Vectors
Test RMSE
Test LLH
Total Time (h)
of Mean
CG
2.0042
0.0382
-1,8249 + 0.0393
22.012 + 0.0047
1.9248 1 0.0340 3.0380 + 0.1826
4.5317
0.4186
-11.973 + 1.4876
21.725 = 0.2588
2.1218 1 0.1021
1.6968 + 0.0424
0,2770 0.0030
-0.0483 + 0.0039
22.466 = 0.0048
0.0516 + 0.0017
0.0851 + 0.0009
AP
0.2743 0.0030
-0.0366 = 0.0034
22.470 = 0.0056
0.0143 = 0.0004
0.0263 = 0.0001
SGD
0.2851
0.0035
-0.1029 1 0.0152
16.722 + 0.0012
0.2252 1 0.0004
0.3906 + 0.0115
0.2735 0.0030
-0.0457
= 0.0045
16.702 = 0.0016
0.0767 = 0.0022
0.1544
0.0039
Table 10: Results on HOUSE with 10 maximum solver epochs (mean + standard error over 10 splits).
HOUSEELECTRIC (n = 844 352, d 11)
Average Residual Norm of Probe Vectors
Test RMSE
Test LLH
Total Time (h)
of Mean
CG
0.8328 0.0215
-4.0695 = 0.2445
32.368 = 0.3881
0.7495 + 0.0224
1.4922 + 0.0826 1.1958 0.0454
0.4519 0.0220
-1.2631 = 0.2168
32.770 = 0.0029
0.6063 0.0252
AP
0.0320 0.0005
1,9051 +: 0.1403
29.090 + 0.0075
0.0215 1 0.0005
0.0598
0.0027 0.0409
0.0292 0.0007
2.3509 +: 0.0654
29.058 + 0.0064
0.0104 4: 0.0002
0.0014
SGD
0.0449
0.0003
0.9635
0.2071
23.646 = 0.0135
0.0409
0.0005
0.1645
0.0085
0.0002 0.0334
1.6021 + 0.0652
23.445 + 0.1611
0.0321 + 0.0006
0.1230 0.0081
26

---

CG
AP
SGD
Test RMSE
1.2
0.03
0.12
0.09
POL
0.8
0.02
0.08
0.08
0.4
0.01
0.04
0.07
0.0
0.00
0.00
0.06
0.39
0.39
ELEV
0.2
0.06
0.26
0.37
Residual Norm Mean
0.35
0.0
0.00
0.00
0.33
0.07
BRKE
0.8
0.02
0.08
0.05
0.1
0.01
0.04
0.03
0.0
0.00
0.00
0.01
2.1
0.15
0.36
0.54
PROT
1.4
0.10
0.24
0.52
0.7
0.05
0.12
0.50
0.0
0.00
0.00
0.48
2.1
0.03
0.12
0.39
KBGG
0.7
0.01
0.04
0.13
0.00
100
25
25
CG
AP
SGD
Marginal Likelihood Steps
Standard estimator
Standard estimator + warm start
Maximum solver epochs marginal likelihood step:
Figure 14: Relative residual norms of the mean at each marginal likelihood step and final test root- mean-square errors using the standard estimator on different datasets. The linear system solver is terminated upon either reaching the tolerance or exhausting maximum number of solver epochs.
CG
AP
SGD
Test RMSE
0.9
0.03
0.12
0.09
POL
0.0
0.02
0.08
0.08
0.3
0.01
0.04
0.07
0.0
0.00
0.00
0.06
0.3
0.06
0.36
0.37
BLEV
0.2
0.04
0.24
0.36
Residual Norm of Mean
0.1
0.02
0.32
0.35
0.0
0.00
0.00
0.31
BIKE
0.6
0.02
0.08
0.06
0.03
0.0
0.00
0.00
0.00
1.8
0.12
0.36
0.51
PROT
1.2
0.08
0.24
0.50
0.6
0.04
0.12
0.49
0.0
0.00
0.00
0.48
0.6
0.03
0.09
0.12
KEGG
0.4
0.02
0.06
0.09
0.2
0.01
0.03
0.06
0.00
0.03
100
CG
AP
SGD
Marginal Likelihood Steps
Pathwise
Pathwise + warm start
Maximum solver epochs per marginal likelihood step:
Figure 15: Relative residual norms of the mean at each marginal likelihood step and final test root- mean-square errors using the pathwise estimator different datasets. The linear system solver is terminated upon either reaching the tolerance or exhausting maximum number of solver epochs.
27

---

CG
AP
SGD
Test LLH
0.6
1.2
1.6
POL
0.4
0.8
0.3
0.2
0.4
1.0
0.0
2.3
Residual Norm of Prube Vectors
BIKE BLEV
0.2
0.18
0.8
-166
0.1
0.09
0.4
3.2
0.0
0.00
0.0
-.8
4.5
0.45
1.2
2.5
3.0
0.30
0.8
2.0
1.5
0.15
0.4
1.5
0.0
0.00
0.0
1.0
0.45
0.9
-0.4
PROT
2
0.30
0.6
-0.6
0.15
0.3
-.88
0.00
0.0
-1.0
0.33
1.2
1.6
KEGG
0.11
0.4
-1.8
0.00
75
25
100
50
100
CG
AP
SGD
Marginal Likelihood Steps
Standard estimator
Standard estimator + warm start
Maximum solver epochs per marginal likelihood step:
Figure 16: Relative residual norms of the probe vectors at each marginal likelihood step and final test ligeliillsooods using the standard estimator on different datasets. The linear system solver terminated upon either reaching the tolerance or exhausting maximum number of solver epochs.
CG
AP
SGD
Test LLH
1.2
0.12
0.6
1.35
POL.
0.8
0.08
1.20
0.4
0.04
0.2
1.05
Residual Norm of Probe Vectors
0.12
-0.38
ELEV
0.1
0.08
0.4
00.99
0.2
0.04
-0.10
0.0
0.00
0.0
00.11
BIKE
0.8
0.08
0.4
1.38
0.0
0.00
0.0
-0.74
1.8
0.12
-0.58
PROT
1.2
0.08
0.2
-0.2
0.6
0.04
0.66
0.0
0.00
0.0
0.70
1.8
0.09
0.6
1.58
KBGG
0.6
0.03
0.2
1.01
0.0
0.00
2.35
CG
Marginal Likelihood Steps
Pathwise
Pathwise + warm start
Maximum solver epochs per marginal likelihoad step:
30
Figure 17: Relative residual norms of the probe vectors at each marginal likelihood step and final test log-likell ligliiklisoods using the pathwise estimator on different datasets. The linear system solver is terminated upon either reaching the tolerance or exhausting maximum number of solver epochs.
28

---

Conjugate Gradients
Res. Norm Mean
Res. Norm Probes
Test RMSE
Test LLH
HOUSE
0.7
1.4
VVW
0.7
-4
0.4
1.0
0.1
-8
10
15
0
5
10
15
BU72
2.7
2.6
-7
www
wwW
0.6
1.0
0.6
14
10
20
30
10
20
30
10
20
30
0
10
20
30
2.2
2.40
2.4
1.0
SONG
0.10
0.6
2.2
10
20
30
10
20
30
10
20
30
0
10
20
30
3DROAD
1.4
1.1
W
1.4
0.8
WWW
0.9
0.9
0.9
1.2
w mm
ww Ww
1.6
20
10
20
0
10
20
30
Marginal Likelihood Steps
Pathwise estimator
Pathwise estimator
warm start
Figure 18: Relative residual norms, test root-mean-squee errors and test log-likelihoods during marginal likelihood optimisation large datasets using the pathwise gradient estimator and conjugate gradients as linear system solver.
Alternating Projections
Res. Norm Mean
Res.
Norm Probes
Test RMSE
Test. LLH
0.04
0.09
0.04
2.6
HOUSE
0.02
0.06
0.03
1.3
0.03
0.02
10
10
0
15
0.06
0.11
0.31
0.03
BUZZ
0.03
0.06
0.29
0.08
0.00
0.01
0.27
0.13
0.06
0.08
0.76
1.11
SONG
0.03
0.01
0.75
1112
0.00
0.00
U.74
-1.33
20
10
20
10
10
30
3DROAD
0.11
0.09
0.12
1.4
0.06
0.06
0.08
0.5
0.01
0.03
0.04
0.4
10
20
30
10
20
30
10
20
30
10
30
Marginal Likelihood Steps
Pathwise estimator
Pathwise estimator
start
Figure 19: Relative residual norms, test root-mean-squee enrors and test log-likelihoods during marginal likelihood optimisation on large datasets using the pathwise gradient estimator and alternat- ing projections as linear system solver.
29

---

Stochastic Gradient Descent
Res. Norm Mean
Norm Probes
Test RMSE
LLH
0.05
0.20
0.05
1.8
HOUSE
0.04
0.15
0.04
1.0
0.03
0.10
0.03
0.2
BUZZ
0.15
0.3
0.29
-0.1
0.00
0.0
0.26
-0.2
10
20
30
10
20
30
10
20
30
10
30
0.10
0.10
0.76
1.11
SONG
0.05
1.12
0.00
0.00
0.74
1.13
10
20
30
10
20
30
10
20
30
10
20
30
3DROAD
0.15
0.10
0.15
1.5
0.09
0.07
0.10
0.03
0.04
0.05
-2.5
10
20
30
10
20
30
10
20
30
10
20
30
Marginal Likelihood Steps
Pathwise estimator
Pathwise estimator warmn
Figure 20: Relative residual norms, test root-mean-squee meaan-squre errors and test log-likelshoods during marginal likelihood optimisation large datasets using the pathwise gradient estimator and stochastic gradient descent as linear system solver.
CG
AP
SGD
Test LLH
225
225k
210k
1.28
POL
150
150k
140k
1.27
75
75k
70k
1.26
0
1.25
ELEV
32
16k
4k
-0.88
Solver Itorations
0
0
0.40
225
180k
600k
2.20
MKE
150
120k
400k
2.14
200%
0
0
2.02
210
36k
39k
-0.58
PROT
24k
263
-0.00
70
12k
13k
00.62
0
0
00.44
150
90k
300k
1.09
KRGG
200k
50
30k
100k
1.07
25
50
100
50
100
0
50
75
100
CG
AP SGD
Marginal Likelihood Steps
Standard
Pathwise
Standard
warm start
Pathwise
warm start
Figure 21: Required number of solver iterations until reaching the tolerance 0.01 at each step of marginal likelihood optimisation and final predictive test log-likelihoods on different datasets. On the KEGGDIRECTED dataset, alternating projections with the standard estimator and without warm starting did not complete the experiment within 24 hours.
30

---

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 main claims made in the abstract and introduction accurately reflect the
contributions of Sections 3 to 5 and empirical results illustrated in Table T and Figure T0
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: Yes, the paper discusses the limitations of our work. In particular, the depen- dence on random features to sample prior functions discussed Section 3 and Figure and the bias due to warm starting is discussed in Section 4 and Theorem
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 "L_imitations" 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 Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
31

---

Answer: [Yes]
Justification: All theoretical claims and results made in the paper are formally derived in Appendix A including all assumptions and proofs.
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: Methods and experiments are described in Sections 3 to with implementation details, pscudocode, and link to source code being provided in Appendix B
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
32

---

reproducible, they 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 provide open access to our source code required to reproduce our experi- ments by uploading
an anonymous repository during the review phase and moving to public repository upon acceptance. The link to the repository is provided in Appendix B In terms of data, we only use publicly available regression datasets, which are also discussed in Appendix B
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please the NeurIPS code and data submission guidelines (HTTPS://NIPS.CC/PPP)
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" 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) for more details.
//NIPS. BOW.U//LL///EEDESSSSEEEESSEESSSS
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments 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 to understand the results?
Answer: [Yes]
Justification: All datasets and train/test splits used in this paper publicly available (7). Experimental settings are explained in Sections to with details being provided in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, 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: All experiments are repeated 10x using the train/test data splits provided by (7). except for the largest datasets, where only 5 splits were used due to computational costs. Tables in the main paper only report the mean over data splits due to space constraints. Tables the appendix report the mean + standard errors over data splits. In figures, mean + standard errors over data splits reported via error bars when appropriate.
33

---

Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi-
dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
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.
OK to report 1-sigma error bars, but one should state 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 out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: Details about compute resources are provided in Appendix B
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of individual experimental runs as 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 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics HHTTPS://NEURIPS /TT.://////.//T//CCCCCCII.C.C..T
Answer: [Yes]
Justification: The research conducted in this paper does not involve any human subjects participants and only uses curated, publicly available datasets (7). Additionally, the research conducted this paper considers generic regression tasks and does not pose a high risk in terms of safety, security, discrimination, surveillance, deception and harassment, environment, human rights, or bias and fairness. While specific applications of regression methods could potentially facilitate some of these harmful consequences, such applications are by no means topic of the research conducted in this paper. Furthermore, implementation details and other information which is required reproduce the research conducted in this paper is thoroughly provided in Appendix B
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 cireumstances that require deviation from the Code of Ethics.
34

---

The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations 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: In this paper, discuss and analyse techniques hyperparameter optimisation in Gaussian processes, which is foundational research and not tied to particular applications deployment. There is no direct path to generating content such fake speeches, images, or videos. Potential harm which could arise from intended use and correct functioning are thus minimal. Potential harm due to intended use but incorrect functioning very unlikely because Gaussian processes are specifically designed to produce calibrated uncertainty estimates. While any general purpose regression algorithm, including Gaussian process regression, technically has the potential to be used in harmful applications or technology, it is by no means topic of this paper. Furthermore, our paper conducts empirical studies on publicly available datasets which do not have a harmful downstream impact and which are
standard benchmark in the literature 9 30 15 16 33
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA 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), fairness 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 not 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 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 train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as 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 or unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a 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 a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper uses curated, publicly available regression datasets 7 and the proposed methods do not pose high risk for misuse or dual-use because they are for generic regression and, in particular, do not process or generate natural language, images, nor videos.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a 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 to usage guidelines or restrictions to access the model or implementing safety filters.
35

---

Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
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 existing assets
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?
Question:
Answer: [Yes]
Justification: The paper uses datasets from the popular UCI regression benchmark made available through GitHub repository provided in Appendix B
7 under the CC BY 4.0 license. This information is
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a partieular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, PAPERSWITHCODE. CWM/CATTO/ATA/ DAATETTS has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For cxisting datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is 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 the documentation provided alongside the assets?
Answer: [NA]
Justification: This paper does not release new assets except the code required to reproduce our experiments, which has already been discussed earlier in this checklist.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing 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: The paper does not involve crowdsourcing research with human subjects.
Guidelines:
36

---

The answer NA means that the paper does not involve crowdsoureing 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 much detail as possible should be included 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/teview based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on 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.
We recognize that the procedures for this may vary significantly between institutions and locations, and 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.
37


Reviews:
[{"review_id": "0c2aFNbLJw", "reviewer": "Reviewer_NcMH", "summary": "This paper investigates several iterative techniques for solving linear systems when applied to the problem of finding GP hyperparameters. Specifically, the following modifications to the method are suggested:\n- A \"pathwise\" sampling estimator for the Hutchinson trace estimator\n- Warm starting linear systems solvers\n\nAfter this, a thorough investigation is made into the effect of these choices on training GP hyperparameters.", "questions": "- Is there any information in the paper that gives a comparison to another method, or ideally a ground truth evaluation of an exact GP?", "limitations": "- No evaluations to any other methods, other than methods with iterative solvers.", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The method follows an existing schema for training GP hyperparameters, and performs an in-depth investigation into how different choices affect the algorithm.", "weaknesses": "The main weakness of this paper is that it does not compare the results of the approximate methods to a reliable exact method. This makes it impossible to determine from the paper whether the predictive performances are actually good.\n\nFor datasets of these sizes (<50k), it is feasible on modern hardware to run a Cholesky on the full dataset. This should be run with BFGS optimisation. Reporting these results will give a ground-truth comparison for how good the resulted test metrics are.\n\nIf these are close, then readers will know how much to value these speedups. If they are distant, readers will be able to determine how much they value a faster approximation that does not get close to optimal predictor.\n\nWithout a comparison to any other method, the value of the overall method cannot be established. I'm happy to discuss this and significantly increase my score if this is addressed sufficiently.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer NcMH for their time to read and review our work. In the following, we want to address their specific concerns and questions:\n\n---\n*The main weakness of this paper is that it does not compare the results of the approximate methods to a reliable exact method. This makes it impossible to determine from the paper whether the predictive performances are actually good.*\n\nWe thank Reviewer NcMH for mentioning BFGS. Indeed, this is the preferred algorithm for GP hyperparameter optimisation when its use is computationally tractable. We initially did not consider it due to our focus on the large-scale setting. **We conducted additional experiments using the Cholesky factorisation + BFGS optimiser on the POL, ELEVATORS, and BIKE datasets**. We were not able to run BFGS on larger datasets because our A100 80GB GPUs did not support the necessary memory requirements. In particular, we used jaxopt [1], and performed optimisation for a maximum of 100 iterations, a maximum of 30 linesearch steps per iteration, and used a stopping tolerance of 1e-5.\n\nAdditionally, **we provide results using the Adam optimiser + Cholesky factorisation** (instead of stochastic gradient estimation) with all other settings being identical to the iterative methods in our paper. This can be considered a pure comparison of the quality of estimated gradients. For both BFGS and Adam, all hyperparameters were initialised at the same values as for the iterative methods in our paper (namely 1.0) and also used the softplus transformation to enforce positive value constraints.\n\nFurthermore, **we also include results for SVGP** [2], a popular variational approximation, taken from Lin et al. [3], who used the same datasets and train / test splits, and performed optimisation until convergence using the Adam optimiser and 3000 inducing points.\n\nThe table below lists the mean over 10 splits of test root-mean-square-errors (RMSE) and test log-likelihoods (LLH), with \"Exact\" referring to Cholesky factorisation + backpropagation. In summary, BFGS achieves marginally better results on the POL and ELEVATORS datasets and significantly better results on the BIKE dataset. However, we conclude that the discrepancy on the BIKE dataset is due to Adam not converging in 100 iterations, because Adam with Cholesky factorisation + backpropagation achieves nearly identical performance compared to the iterative methods (which also use Adam). Therefore, the performance gap cannot be linked to the gradient estimation via iterative methods, but rather to the different optimisers. However, in the large-scale setting, memory requirements and stochastic gradient estimation make BFGS infeasible. Therefore, we conducted all our experiments with Adam. In comparison, SVGP generally yields worse results, particularly in terms of log-likelihood.\n\n|  |  |  | RMSE |  |  |  |  |  | LLH |  |  |\n|---|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|\n|  | POL | ELEV | BIKE | PROT | KEGG |  | POL | ELEV | BIKE | PROT | KEGG |\n| CG | 0.0750 | 0.3550 | 0.0326 | 0.5024 | 0.0837 |  | 1.2682 | -0.3856 | 2.1500 | -0.5871 | 1.0818 |\n| CG + ws | 0.0754 | 0.3562 | 0.0326 | 0.4909 | 0.0837 |  | 1.2716 | -0.3868 | 2.0674 | -0.6210 | 1.0818 |\n| CG + pw | 0.0750 | 0.3550 | 0.0327 | 0.5026 | 0.0837 |  | 1.2681 | -0.3856 | 2.1508 | -0.5871 | 1.0816 |\n| CG + ws + pw | 0.0758 | 0.3558 | 0.0329 | 0.4912 | 0.0836 |  | 1.2666 | -0.3856 | 2.0615 | -0.6214 | 1.0819 |\n|  |  |  |  |  |  |  |  |  |  |  |  |\n| AP | 0.0750 | 0.3550 | 0.0326 | 0.5024 | - |  | 1.2682 | -0.3856 | 2.1504 | -0.5871 | - |\n| AP + ws | 0.0754 | 0.3562 | 0.0325 | 0.4907 | 0.0837 |  | 1.2715 | -0.3868 | 2.0668 | -0.6214 | 1.0820 |\n| AP + pw | 0.0750 | 0.3550 | 0.0326 | 0.5027 | 0.0837 |  | 1.2681 | -0.3856 | 2.1503 | -0.5871 | 1.0818 |\n| AP + ws + pw | 0.0758 | 0.3558 | 0.0330 | 0.4912 | 0.0836 |  | 1.2666 | -0.3856 | 2.0616 | -0.6213 | 1.0817 |\n|  |  |  |  |  |  |  |  |  |  |  |  |\n| SGD | 0.0750 | 0.3550 | 0.0326 | 0.5026 | 0.0837 |  | 1.2681 | -0.3855 | 2.1535 | -0.5871 | 1.0816 |\n| SGD + ws | 0.0754 | 0.3562 | 0.0324 | 0.4894 | 0.0837 |  | 1.2708 | -0.3868 | 2.0692 | -0.6268 | 1.0822 |\n| SGD + pw | 0.0750 | 0.3550 | 0.0327 | 0.5027 | 0.0837 |  | 1.2682 | -0.3855 | 2.1524 | -0.5878 | 1.0821 |\n| SGD + ws + pw | 0.0757 | 0.3558 | 0.0332 | 0.4911 | 0.0839 |  | 1.2678 | -0.3854 | 2.0562 | -0.6217 | 1.0725 |\n|  |  |  |  |  |  |  |  |  |  |  |  |\n| Exact (BFGS) | 0.0714 | 0.3495 | 0.0317 | - | - |  | 1.2765 | -0.3693 | 3.2685 | - | - |\n| Exact (Adam) | 0.0750 | 0.3550 | 0.0326 | 0.5025 | 0.0837 |  | 1.2683 | -0.3856 | 2.1501 | -0.5868 | 1.0813 |\n| SVGP | 0.10 | 0.37 | 0.08 | 0.57 | 0.10 |  | 0.67 | -0.43 | 1.21 | -0.85 | 0.54 |\n\n---\n*Is there any information in the paper that gives a comparison to another method, or ideally a ground truth evaluation of an exact GP?*\n\nCurrently, **Figures 5 and 8** visualise the influence of the pathwise estimator and warm starts by comparing their optimisation trajectories to exact optimisation using Cholesky factorisation + backpropagation (and otherwise identical configurations).\n\nAdditionally, **Figures 11-13 in Appendix C** also compare iterative methods to exact optimisation with Cholesky factorisation + backpropagation by visualising the optimisation trajectories of selected hyperparameters. Furthermore, we are happy to include the quantitative results from the table above in the camera-ready version.\n\n---\nFinally, we hope that we successfully addressed all concerns and questions, and kindly encourage Reviewer NcMH to increase their score or to reach out again with follow-up questions if there are any unresolved concerns.\n\n---\n[1] Blondel et al. (2021), \"Efficient and Modular Implicit Differentiation\", *arXiv*.\n\n[2] Hensman et al. (2013), \"Gaussian Processes for Big Data\", *Uncertainty in Artificial Intelligence*.\n\n[3] Lin et al. (2024), \"Stochastic Gradient Descent for Gaussian Processes Done Right\", *International Conference on Learning Representations*.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ajz2GJdDF5", "reviewer": "Reviewer_Xgts", "summary": "The paper presents approaches that can speed up solving linear systems arising in the GP regression problem. The two basic ideas are a warm start and limiting the computation budget. Another novelty is the pathwise gradient estimator, which leads to fewer iterations needed for convergence. The detailed numerical experiments confirm the performance of the suggested approaches and raise questions on the connection between the accurate solving of linear systems and the performance of the GP regression method.", "questions": "1. The authors mentioned in the Appendix that they use precondition based on the Cholesky factorization. Please provide more details since, typically, such preconditioners, their structure, and properties significantly affect the performance of CG. \n2. Please add bold labeling the top timing and test log-likelihood values in Table 1 to simplify parsing such a large number of values\n3. Why are large datasets (lines 219-235) not included in Table 1 or a similar table where the runtime gain can be observed? Please add these results, too. The large datasets can illustrate the gain from the introduced approach even better.", "limitations": "The authors do not provide explicit limitations of the suggested approach.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The novel approach for gradient estimation is proposed and demonstrates significant speed-up.\n2. Different techniques to accelerate solvers for linear systems are proposed, and their efficiency is confirmed via extensive experimental evaluation. This evaluation includes target solvers CG, AP, and SGD.\n3. The manuscript is well-prepared, the motivation and suggestions are clear, and the results are well-explained.\n4. The presented techniques can be easily incorporated into the existing packages for GP regression and used to solve practically important problems.", "weaknesses": "In this paper, the authors use some heuristics to speed up linear solvers and suggest a new gradient estimator based on the transformation of the normal distribution. Although the suggested tricks demonstrate improvement, their theoretical interpretation and the possibility of obtaining even better results remain unclear. So, a more detailed discussion of further research directions related to the presented approaches would be interesting. The theoretical bounds for the potential gain obtained from linear solvers while preserving test error can make the work more solid and fundamental.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer Xgts for their time to read and review our work, and are grateful to hear that our \"motivation and suggestions are clear\", our manuscript is \"well-prepared\", and our experimental evaluation is \"extensive\". In the following, we want to address their specific concerns and questions:\n\n---\n*So, a more detailed discussion of further research directions related to the presented approaches would be interesting.*\n\nAt the end of Section 6, we mention that an interesting direction for future research would be to investigate why low relative residual norms do not always translate to good predictive performance. To elaborate on that, we believe that interesting research questions are:\n- Which metric, norm, or quantity actually corresponds to predictive performance (other than predictive performance itself)?\n- Can we develop a stopping criterion based on this particular quantity (instead of relative residual norm) to optimise for predictive performance?\n\n**We are happy to include this discussion at the end of Section 6 in the camera-ready version of our manuscript.**\n\n---\n*The authors mentioned in the Appendix that they use precondition based on the Cholesky factorization. Please provide more details since, typically, such preconditioners, their structure, and properties significantly affect the performance of CG.*\n\n**We will update our manuscript to provide the following details (and references) about the pivoted Cholesky preconditioner in Appendix B:**\n\nThe pivoted Cholesky factorisation [1] constructs a low-rank approximation to the full (but intractable) Cholesky factorisation by permuting rows and columns to prioritise the entries with the highest variance in a greedy way. Given the extensive and successful use of this preconditioner in the context of iterative Gaussian processes [2, 3, 4], we follow and refer to existing work.\n\n---\n*Please add bold labeling the top timing and test log-likelihood values in Table 1 to simplify parsing such a large number of values.*\n\n**We will update our manuscript to include a summary of the following explanation in the caption of Table 1:**\n\nIn Table 1, we compare different iterative methods when solving until the tolerance is reached. This results in all methods performing (nearly) identical in terms of log-likelihood (except for a few cases where the bias of the pathwise estimator due to random features leads to a slightly different outcome). Therefore, we did not bold these log-likelihood values because almost every value in the table would be bold. While the runtimes indeed differ, the comparison mainly evaluates the effects of warm starts and the pathwise estimator *per solver*. Thus, we bold the relative speed-ups per solver instead of the absolute runtimes.\n\n---\n*Why are large datasets (lines 219-235) not included in Table 1 or a similar table where the runtime gain can be observed? Please add these results, too. The large datasets can illustrate the gain from the introduced approach even better.*\n\n**We do provide runtimes for large datasets in Tables 7-10 in Appendix C, and we are happy to include these in main paper for the camera-ready version.**\n\nTable 1 considers small datasets and presents the predictive performances of runtimes of solving until reaching the tolerance. In this setting, the final performance is basically fixed and gains in runtime are the main focus. In contrast, for the large datasets, solvers run for a fixed maximum number of epochs, implying a (roughly) fixed maximum amount of time. While in some cases solvers eventually reach the tolerance in fewer epochs and less time, this is not the main emphasis of these experiments. Instead, the main goal is to show that warm starts and the pathwise estimator lead to better performance under a fixed compute budget due to accumulation of solver progress, which is visualised in Figure 10.\n\n---\nFinally, we hope that we successfully addressed all concerns and questions, and kindly encourage Reviewer Xgts to consider increasing their score or to reach out again with follow-up questions if there are any unresolved concerns.\n\n---\n[1] Harbrecht et al. (2012), \"On the low-rank approximation by the pivoted Cholesky decomposition\", *Applied Numerical Mathematics*.\n\n[2] Gardner et al. (2018), \"GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration\", *Neural Information Processing Systems*.\n\n[3] Wang et al. (2019), \"Exact Gaussian Processes on a Million Data Points\", *Neural Information Processing Systems*.\n\n[4] Wu et al. (2024), \"Large-Scale Gaussian Processes via Alternating Projection\", *Artificial Intelligence and Statistics*.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8rhXjKh7Ik", "reviewer": "Reviewer_89S9", "summary": "The paper considers the problem of simultaneously fitting a Gaussian process (GP) to data along with determining the hyperparameters (kernel width, noise variance) for the GP.  The overall algorithm is standard, consisting of an outer loop, which is a simple gradient update to the hyperparameters, and an inner loop that fits the GP to the data.  Fitting the data of course requires a system solve with the kernel matrix.  Computing the gradient (or an approximation thereof) for the hyperparamters requires multiple system solves.\n\nThe paper gives some numerical experiements and some rough reasoning to argue that tweaking the trace estimator used for the gradient estimate and warm starting the linear solver saves a tangible amount of computation in this process.", "questions": "How are the $\\xi$ in (11) generated and what is the computational cost relative to the system solves?", "limitations": "none", "rating": 5, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "Hyperparameter optimization for fitting GPs is an important and prevalent problem in machine learning.  \n\nThe numerical experiments are thorough and appear to be carefully done, and lend good supporting evidence to the efficacy of the method.", "weaknesses": "It is hard for me to see what the really new ideas are in this paper.  The the authors say, the most critical part of this whole process is the linear system solver.  They do not actually propose any new type of solver, but rather show that by modifying the way the gradient of the hyperparameters is estimated and using a warm start, the linear solver they happen to be using will take fewer iterations.  The idea of warm starting a linear solver using previous solutions is pervasive throughout numerical analysis, so it is no surprise that it gives you some gains.  The main novelty in the paper seems to be changing the trace estimator; that is, using equation (11) in place of (7).  \n\nSomething which is not really discussed, but seems important, is that the z in (7) are easier to generate than the $\\xi$ in (11).  The z can be random vectors with independent entries, the $\\xi$ need a prescribed covaraince structure that changes at every iteration.  Since this covariance matrix (the kernel matrix) is large, taking a random draw of the $\\xi$ seems nontrivial.  How this works, and why it does not have a significant effect on the number of computations needed, needs to be spelled out more clearly.\n\nThe numerical experiements are on interesting problems, but the paper would be stronger if the mathematical arguments were more precise.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer 89S9 for their time to read and review our work, and are excited to hear that our experiments are \"thorough\" and \"carefully done\". In the following, we want to address their specific concerns and questions:\n\n---\n*It is hard for me to see what the really new ideas are in this paper. [...] They do not actually propose any new type of solver [...]. The idea of warm starting a linear solver using previous solutions is pervasive throughout numerical analysis [...]. The main novelty in the paper seems to be changing the trace estimator [...].*\n\nIndeed, **we do not propose a new solver, nor do we claim to**; the title of our paper is “Improving Linear System Solvers […]”.\n\nWe want to emphasise that a major (if not our main) contribution consists of an **extensive empirical evaluation** which\n- quantifies the effectiveness of using **warm starts** and the **pathwise estimator**;\n- across **three different kinds of linear system solvers**;\n- for five small ( < 50k) datasets in the regime of **running solvers until convergence**;\n- and four large ( up to 1.8M) datasets in the regime of **limited compute budgets** of 10, 20, 30, 40, and 50 epochs.\n\nAs a result, we demonstrate **speed-ups of up to 72x** when running solvers until reaching the tolerance, and achieve **up to 7x lower average relative residual norms** under a limited compute budget. We believe that this evaluation will be of **great value to future researchers and practitioners**, and **should not be underrated** when judging the novelty and contributions of our work.\n\nWhile we agree that warm starts are commonly used, we also believe that **warm starts have not been considered for stochastic trace estimation in the context of marginal likelihood optimisation for GPs, nor compared across various solvers before** (if yes, we would appreciate a reference to the literature). Thus, we do not claim warm starts per se as a novel contribution, but rather our application and theoretical / empirical analysis in this particular setting.\n\n---\n*Something which is not really discussed, but seems important, is that the $z$ in (7) are easier to generate than the $\\xi$ in (11). The $z$ can be random vectors with independent entries, the $\\xi$ need a prescribed covaraince structure that changes at every iteration. Since this covariance matrix (the kernel matrix) is large, taking a random draw of the $\\xi$ seems nontrivial. [...] How are the $\\xi$ in (11) generated and what is the computational cost relative to the system solves?*\n\nWe agree that the $z$ are easier to generate than the $\\xi$, because the latter effectively requires a matrix square root of the kernel matrix $K$ of size $n \\times n$. To draw $\\xi$ efficiently, we use random features, which produce a low-rank square root with an unbiased approximation to the kernel, reducing the asymptotic time complexity from $\\mathcal{O}(n^3)$ to $\\mathcal{O}(nm)$, where $m$ is the number of random features. Since $m$ is constant and much smaller than $n$ (e.g. we used $m$ = 2k while $n$ > 1M), this effectively becomes an $\\mathcal{O}(n)$ operation, which is faster than CG with $\\mathcal{O}(n^2)$ per iteration, and similar to AP / SGD with $\\mathcal{O}(n)$ per iteration. Existing work used random features for this purpose and demonstrated strong empirical performance despite this approximation [1, 2]. A discussion of alternatives can be found in Section 4 of [1].\n\n**We mention random features in Section 2, lines 70-71, and discuss them further in Section 3, lines 143-149, and Figure 5. Additionally, a detailed description of how to sample random probe vectors is provided in Appendix B, lines 486-501.**\n\n---\n*The numerical experiements are on interesting problems, but the paper would be stronger if the mathematical arguments were more precise.*\n\nWe thank Reviewer 89S9 for describing the problems in our experiments as \"interesting\". To clarify any potential misconceptions about our mathematical arguments, we summarise them here:\n\nIn Section 3, we calculate the expected squared distance between zero initialisation and the solution of the linear solver. In particular, we show that, for the standard estimator, this distance depends on the spectrum of the regularised kernel matrix, while it is constant for the pathwise estimator. Thus, for the standard estimator, this distance tends to increase as the spectrum changes during hyperparameter optimisation, leading to substantially increased number of required solver iterations. In contrast, for the pathwise estimator, the constant distance translates to a roughly constant (and significantly lower) number of required solver iterations.\n\nIn Section 4, we investigate the concern that warm starts require deterministic probe vectors which introduce correlations into subsequent gradient estimates. Despite individual estimates being unbiased, these correlations introduce bias into the optimisation. Theorem 1 shows that, under certain assumptions, the marginal likelihood after optimisation with deterministic probe vectors converges in probability to the marginal likelihood of a true maximum as the number of probe vectors goes to infinity. The proof is provided in Appendix A, and we are working on an updated, simpler proof with better constants.\n\n**After the clarifications above, we are happy to provide further details upon request and would appreciate it if Reviewer 89S9 has specific suggestions for further mathematical arguments which they think might strengthen the paper.**\n\n---\nFinally, we hope that we successfully addressed all concerns and questions, and kindly encourage Reviewer 89S9 to increase their score or to reach out again with follow-up questions if there are any unresolved concerns.\n\n---\n[1] Wilson et al. (2021), \"Pathwise Conditioning of Gaussian Processes\", *Journal of Machine Learning Research*.\n\n[2] Lin et al. (2023), \"Sampling from Gaussian Process Posteriors using Stochastic Gradient Descent\", *Neural Information Processing Systems*.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "sUB0pHEbA9", "reviewer": "Reviewer_upE1", "summary": "This paper focuses on solving linear systems arising from large-scale Gaussian process hyperparameter optimisation.\n\nThe first technique proposed in this paper is a reformulation of the linear system for log determinant gradient estimation, which is called the \"pathwise estimator\".\nThe authors argue that this reformulation makes the all-zero initialization closer to the solution.\nThis reformulation also allows them to reuse the linear system solutions in the prediction time by pathwise conditioning (Wilson et al., 2021).\n\nThe second technique proposed in this paper is warm starting the linear solvers.\nThough, warm starts require using a set of **fixed** probe vectors, which necessarily introduces biases to the hyperparameter optimization.\nThe authors demonstrate empirically that the bias does not hurt in practice.", "questions": "- The conjugate gradient method estimates both the log marginal likelihood (by tridiagonalization) and its gradient.\nI am wondering if using the pathwise estimator for conjugate gradient produces an estimate for the log marginal likelihood.\nThis is useful monitoring the progress of hyperparameter optimization.", "limitations": "NA", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The writing is very clear.\n- By combining several techniques, this paper achieves large speed-up factors over a wide range of datasets.\nBesides the speed-up, it appears that the pathwise estimator with warm start also stabilizes the linear systems solvers in some cases.\n- The empirical evaluations are extensive with datasets ranging from 10 thousand data points to 1.8 million data points.\nThe extensive experiments with recent scalable GP methods could provide tips and tricks to practitioners.", "weaknesses": "The main weakness is that the techniques are not super novel per se.\n\nThere exists something similar to the pathwise estimator proposed in this paper.\nFor example, see Equation (7) in Gardner et al. (2018), where the probe vectors are sampled from $\\mathcal{N}(\\mathbf{0}, \\mathbf{P})$ with $\\mathbf{P} \\approx \\mathbf{H}_{\\mathbf{\\theta}}$.\nThough, the original motivation in Gardner et al. (2018) is different from this paper.\nNevertheless, they do appear to be similar at least on a conceptual level.\n\nIn addition, warming starting the linear solvers has been used in practice by Artemev et al. (2021) and Antoran et al. (2023).\n\nHowever, the main contribution of this paper is combining all techniques together and show their effectiveness on a wide range of datasets.\nFuture work could benefit from the the insights and empirical evaluations in this paper.\n\nMinor Weaknesses:\n- Theorem 1 is not very interesting.\nThis theorem is based on the fact that the gradient bias goes to zero as the number of probe vectors goes to infinity.\nMy impression is that the authors wrote this theorem to make the paper looks more \"theoretical\".\nFrankly, I think removing this theorem makes the paper look cleaner.\nA middle ground would be just presenting how fast the bias goes to zero as $s \\to \\infty$.\n\n- The name \"pathwise estimator\" might be somewhat confusing to some readers.\nThis name is usually reserved for the reparameterization trick gradient, which is evidently not the case for this paper.\nI suggest the authors think about alternative names for this estimator.\nAn possible option is to simply rephrase Section 3 as \"a reformulation of the linear systems\".\n\n- Line 172: footnote should be after the period.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer upE1 for their time to read and review our work, and are delighted to hear that our \"writing is very clear\", and \"empirical evaluations are extensive\". In the following, we want to address their specific concerns and questions:\n\n---\n*There exists something similar to the pathwise estimator proposed in this paper. For example, see Equation (7) in Gardner et al. (2018), where the probe vectors are sampled from $\\mathcal{N}(0, P)$ with $P \\approx H$. Though, the original motivation in Gardner et al. (2018) is different from this paper. [...]*\n\nWe agree that the probe vectors used by Gardner et al. (2018) have **conceptual similarities but the motivation is different**. In particular, Gardner et al. (2018) use $z \\sim \\mathcal{N}(0, P)$ probe vectors, where $P$ is constructed using a low-rank pivoted Cholesky decomposition, **to implement the preconditioner**, allowing their Lanczos quadrature and trace estimation to return the correct values. In contrast, our pathwise probe vectors are sampled using random features, and we used the pivoted Cholesky preconditioner *in addition* to the pathwise probe vectors.\n\nFurthermore, **we leverage the solution of the pathwise probe vector systems to construct posterior samples via pathwise conditioning** [1], which has not been done by Gardner et al. (2018). We do this by adapting the standard trace estimator to share linear solves with pathwise conditioning, effectively amortising expensive computations. In contrast, Gardner at al. (2018) solve additional linear systems to obtain the posterior predictive distribution.\n\n---\n*In addition, warming starting the linear solvers has been used in practice by Artemev et al. (2021) and Antoran et al. (2023).*\n\nIndeed, and **we state exactly this in our footnote on page 6**. However, we also point out in our footnote that **they considered slightly different settings**.\n\nIn particular, **Artemev et al. (2021) only apply warm starts to the mean system and not to any probe vectors**, because they do not perform stochastic trace estimation. Arguably, applying warm starts to the probe vectors is non-trivial due to the introduced correlations and bias, as discussed in Section 4.\n\n**Antoran et al. (2023) used warm starts for linear solvers in the context of finite-dimensional linear models with Laplace approximation rather than Gaussian processes.**\n\n---\n*Theorem 1 is not very interesting. This theorem is based on the fact that the gradient bias goes to zero as the number of probe vectors goes to infinity. [...] Frankly, I think removing this theorem makes the paper look cleaner. A middle ground would be just presenting how fast the bias goes to zero as $s \\to \\infty$.*\n\nTo avoid misconceptions about Theorem 1, we would like to point out that **it is not as simple as it may seem**. In particular, to obtain the result, **it is *not* sufficient to show that the error of the gradient estimator goes to zero as the number of probe vectors goes to infinity** (which would indeed not be very interesting), since the rate of convergence of the error could technically be different in different parts of the optimisation landscape, leading to asymptotically biased optima. Therefore, the actual result was quite non-trivial to obtain. The complete proof in Appendix A gives a non-asymptotic result with a convergence rate. To avoid excessive notation, we only provided a simplified asymptotic statement in Section 4.\n\nWe are currently working on a simplified proof with an improved convergence rate. **In general, we believe that this is a non-trivial result and would be inclined to keep Theorem 1 in the main paper. However, we are open to keeping it in Appendix A if Reviewer upE1 feels strongly about this.**\n\n---\n*The name \"pathwise estimator\" might be somewhat confusing to some readers. This name is usually reserved for the reparameterization trick gradient, which is evidently not the case for this paper. [...]*\n\nThis is an interesting observation! In fact, **our pathwise estimator inherits its name from pathwise conditioning** [1], **which can be interpreted as reparameterising a sample from the GP posterior as a deterministic transformation of a sample from the GP prior** (and thus its name)! To make the connection to the reparameterisation trick more evident, note that the sample from the GP prior itself can again be reparameterised as an affine transformation of a standard normal random variable.\n\n**We will update our manuscript to include a brief version of the explanation above.**\n\n---\n*The conjugate gradient method estimates both the log marginal likelihood (by tridiagonalization) and its gradient. I am wondering if using the pathwise estimator for conjugate gradient produces an estimate for the log marginal likelihood. This is useful monitoring the progress of hyperparameter optimization.*\n\n**Currently, the pathwise estimator does not estimate the marginal likelihood** (although this could be a future research endeavor!). However, it produces **posterior predictive samples via pathwise conditioning** [1], which can be used to **monitor the progress of hyperparameter optimisation by evaluating the predictive performance**. Conventionally, evaluating predictions is expensive because it requires additional linear solves. This becomes intractable for existing iterative methods without relying on additional approximations.\n\n---\n*Line 172: footnote should be after the period.*\n\n**We thank Reviewer upE1 for pointing this out and will correct this for the camera-ready version.**\n\n---\nFinally, we hope that we successfully addressed all concerns and questions, and kindly encourage Reviewer upE1 to consider increasing their score or to reach out again with follow-up questions if there are any unresolved concerns.\n\n---\n[1] Wilson et al. (2021), \"Pathwise Conditioning of Gaussian Processes\", *Journal of Machine Learning Research*.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers, ACs, SACs, PCs, organisers, and other volunteers for their time. NeurIPS 2024 would not be possible without their generous time commitments! In the following, we first give an overview of reviewer comments and then discuss the main concerns raised by reviewers.

---
**Overview**

- We are generally happy about the **large amount of positive feedback** received from the reviewers.
- In particular, we are glad to hear that our **"motivation and suggestions are clear"**, that our manuscript is **"well-prepared"**, and that our **"writing is very clear"**.
- In terms of our empirical contributions, we are excited to hear that reviewers described our **experimental evaluation** as **"extensive"**, **"thorough"**, and **"carefully done"**.
- In the following, we summarise the main concerns raised by reviewers, and juxtapose them with our perspective.

---
**No Comparison to Exact Optimisation**
- According to Reviewer NcMH: **"The main weakness of this paper is that it does not compare the results of the approximate methods to a reliable exact method."**
    - Therefore, we provide **new results for BFGS and Adam optimisers using Cholesky factorisation + backpropagation**, as requested by Reviewer NcMH.
    - **Tabulated results**, **experimental details**, and a **comprehensive discussion** are provided in the individual response to Reviewer NcMH.

---
**Novelty and Contributions**
- We want to emphasise that a major (if not our main) contribution consists of an **extensive empirical evaluation**.
    - Our evaluation quantifies the effectiveness of using **warm starts** and the **pathwise estimator**, separately and combined;
    - across **three different kinds of linear system solvers**;
    - for five small ($n$ < 50k) datasets in the regime of running solvers **until convergence**;
    - and four large ($n$ up to 1.8M) datasets in the regime of **limited compute budgets** of 10, 20, 30, 40, and 50 epochs.
    - As a result, we demonstrate **speed-ups of up to 72x** when running solvers until reaching the tolerance, and achieve **up to 7x lower average relative residual norms** under a limited compute budget.
    - We believe that this evaluation will be of **great value to future researchers and practitioners**, and **should not be underrated** when judging the novelty and contributions of our work.

- Some reviewers were **concerned about the novelty of warm starts**, because they are commonly used in the context of optimisation.
    - While we agree about their common use, we would like to highlight that, to the best of our knowledge, **warm starts have not been explored in the context of stochastic trace estimation for marginal likelihood optimisation in Gaussian processes**.
    - In this setting, warm starts **introduce correlations into otherwise unbiased gradient estimates**, leading to bias in the optimisation trajectory.
    - We provide a **theorem which supports the use of warm starts** in this setting, and demonstrate empirically that **warm starts are extremely effective in practice**.

- Reviewer upE1 pointed out that **concepts similar to our proposed pathwise estimator can be found in the existing literature**, such as Gardner et al. (2018) who sample probe vectors $z \sim \mathcal{N}(0, P)$ to implement a preconditioner $P^{-1}$.
    - While we also sample probe vectors from a similar distribution, we agree with Reviewer upE1 that our **motivation is different**.
    - In particular, we introduce a **novel connection based on pathwise conditioning**, which **amortises computations** between optimising hyperparameters and drawing posterior samples, and enables **efficient evaluation of predictive performance** during training.
    - Furthermore, we provide a **compelling theoretical argument** and **matching empirical evidence** which explain why the **standard estimator requires more solver iterations than the pathwise estimator** to reach the tolerance.

---
**Limitations**

- Reviewers mentioned that our **manuscript currently does not contain an explicit section which discusses the limitations** of our methods.
    - While this is true, **we do acknowledge limitations and discuss them** throughout the paper (as stated in the NeurIPS checklist).
    - In particular, we consider the **dependence on random features** to efficiently sample from a Gaussian process prior as the main limitation of the pathwise estimator
    - We briefly mention this dependence in **Section 2**, discuss its consequences in **Section 3** and **Figure 5**, and provide further details in **Appendix C**.
    - Additionally, we consider the **bias introduced into the optimisation trajectory** due to deterministic probe vectors required for warm starts as another limitation.
    - Theorem 1 shows that **optimal performance can theoretically still be achieved** despite this bias, and Figure 8 illustrates that this **bias is negligible in practice**.
    - We are happy to provide an **explicit discussion of these limitations in its own section** if the reviewers believe that this would be a valuable addition to the manuscript.

---

Last but not least, we again thank everyone involved for their time and contributions to making NeurIPS 2024 possible!

Meta-review: I think this is a solid paper with some weaknesses that are generally well addressed overall in the rebuttal. To summarize my thoughts:

- The inclusion of exact results where possible clearly demonstrates that the authors' speed ups are "valid" in the sense that their approximate methods are achieving nearly exact results in terms of RMSE and LLH. Additionally, in my experience in general is that the strength of BFGS is often oversold compared to Adam provided care is taken with GPs, and the results presented here are roughly in line with what I've seen.

- I don't generally see comparing a method whose sole purpose is to speed up iterative methods with the iterative methods it speeds up as a weakness. Given that the authors' results are ultimately the same as the underlying iterative method, whether or not these methods perform better than e.g. SVGPs is a total red herring: that is probably the fault of the underlying iterative method, not the authors' techniques in this paper. I ultimately don't think that this paper is the correct venue to litigate the merit of iterative methods as a whole.

- On the weaknesses side (and I mention this in the possibly vain hope that the authors consider this for the camera ready deadline), I'd like to note that a trivial way around the issue of comparing to exact GPs would be to simply run on larger datasets. In particular, I'm surprised that -- given fairly large speedups to already relatively scalable methods -- we aren't running on much larger datasets here? Many of the original papers (e.g., [32]) run on much larger datasets. Why aren't we showing impressive 72x speedups on massive datasets that should now take very little time to train on?

Overall I think the paper is good and I'm strongly in favor of acceptance, but I *do* think the paper could be *great* if the authors would just run on larger datasets and solve point 3.
Decision: Accept (poster)