Stochastic Gradient Descent for Gaussian Processes Done Right

Thank you very much for your review! We are delighted to hear that our results are significant—thank you very much for this comment. Let us address your questions:

1. Early behaviour in optimisation iterations

In Fig. 1 you note that the primal gradient makes more progress in K^2-norm while the dual one in K-norm... in the left panel of Fig. 1 it seems that for a few iterations in the beginning of optimisation, the primal gradient was also making more progress than dual in the K-norm.

Regarding Figure 1, specifically the ordering of performance right at the beginning of optimisation: we’ve thought about this, and we genuinely do not know! We think that this behaviour is relatively rare, as it occurs way out of convergence, with relatively few iterations spent in this regime, which should make it unlikely to affect behaviour in practice very much.

2. Kernel hyperparameter optimisation / Bayesian model selection

The GP hyper-parameters ... are typically estimated by maximising the marginal log-likelihood using gradient-based methods. Do you think it could be possible run the gradient-based hyper-parameters inference jointly with the posterior inference that you discussed in this paper?

Our method can indeed be easily combined with standard approaches to kernel hyperparameter optimisation. We have added a demonstration of this in Appendix C.7 (we use an approach similar to Gardner et al. 2018), and include some empirical results to that end where we match performance of marginal likelihood optimisation with exact linear system solves (see the figure here). Thank you for pointing out that this may be interesting to our readers to highlight!

Thank you for your thorough review. We are glad you think the paper is easy to follow and well-written.

First of all we would like to address the overall comment regarding our technical contributions.

...technical contributions seem to be below the bar of ICLR

All applied kernel/GP work we are familiar with, which uses SGD methods, makes suboptimal algorithmic design choices (e.g. Lin et al. 2023 (NeurIPS2023) and Dai et. al. 2014 (NeurIPS2014)). On this basis, we believe that a simple algorithm that combines recent theoretical advancements in a principled way, and is the first to reliably outperform conjugate gradient-based algorithms used in state-of-the-art GP software packages, can be of value to the ICLR community. We consider the simplicity of our approach a strength!

We go on to address your remaining concerns and questions.

1. Instability and divergence of baselines.

More figures/tables to explicitly show the weakness/instability of the baseline methods are expected.

We have added a number of further comparisons where we show failure modes of the baselines we consider, namely SGD from Lin et al. 2023, CG, and SVGP.

1. In Figures 12 and 13 in Appendix C.6, also available here and here, you can see how the primal SGD baseline diverges very quickly at even moderately high learning rates.
2. In Figure 6 in Appendix C.1 (here), you can see failure modes for CG (on the infill asymptotics task) and SVGP (large-domain asymptotics task). Our dual SGD method performs robustly in both tasks.
3. In Figure 11 in Appendix C.5 (here), you can see that CG struggles in ill-conditioned settings, where the dual SGD method is quite robust.

2. Algorithmic insights

Sharing more insights into the algorithm settings, such as the choice of geometric averaging, the effect/influence on the sparsity of the unbiased estimator \hat(g)(\alpha), etc, are expected.

Sparsity of the random coordinate estimate plays no role in our recommendation of this gradient estimate. Rather we point it out to help the reader understand what this estimate looks like.

Regarding averaging, our main point is that iterate averaging may not be necessary for convergence under multiplicative noise (Section 3.3), and can be completely done away with. Indeed, arithmetic iterate averaging slows down convergence empirically (see Figure 3).

However, we find that geometric averaging, which is much softer a form of averaging (places much more weight on recent iterates), is useful in practice. We thus recommend its use.

3. Theoretical analysis

A theoretical convergence analysis is expected (not only some figures).

Our paper is targeted at GP researchers and practitioners. While using SGD for GP inference has recently been found to be very promising (Lin et al. 2023 and Dai et. al. 2014), previous work makes a number of suboptimal choices when it comes to the details of the algorithm's formulation. Our contribution is handling these details carefully, highlighting to researchers and practitioners how to make SGD inference even more effective, particularly on harder problems like the graph kernel comparison.

In the full-batch (non stochastic) case, the analysis of the algorithm is standard (Nesterov 1983), and it is known that the dual objective will result in faster convergence than its primal counterpart due to conditioning. See, for instance Section 3.7, page 289, of Bubeck 2015, a textbook on convex optimisation.

Our full algorithm makes use of an estimator with multiplicative noise. The convergence of Stochastic GD under multiplicative noise is currently at the forefront of theoretical research and beyond the scope of our more practical work. We have extended our discussion of pertinent work in this space in Section 3.4. We hope that our work's demonstration of the empirical performance achievable with SGD under multiplicative noise will motivate research in this area by the optimisation community.

(continued below)

8. Effect of batch size

How do different batch sizes affect the overall convergence?

We conducted additional experiments with varying batch sizes. Please see our figures here and here and Appendix C.4 for further details.

Smaller batch sizes degrade the final performance of the additive noise estimator because the amount of additive noise is inversely proportional to the batch size. In contrast, the multiplicative noise estimator produces nearly identical optimisation traces with different batch sizes, as long as the batch size is not small enough to cause divergence. Therefore, increasing the batch size consistently improves the final performance of the additive noise estimator while it does not impact the performance of the multiplicative noise estimator, given that the algorithm converges.

In terms of performance vs time, a larger batch size requires additional compute per step. Thus, the additive noise estimator performance is proportional to the batch size and the multiplicative noise estimator performance per time is inversely proportional to the batch size, as long as the algorithm converges. Therefore, the multiplicative noise estimator should use the smallest batch size which does not cause divergence. Selecting such a batch size a priori is, in general, non-trivial. We have included this discussion in Appendix C4.

9. Non-stationary kernels

It is better to add a test where samples are generated by a non-stationary kernel, to show the ability of the random coordinate estimate. (to distinguish with the random Fourier features)

Thank you for the question. Our optimisation method does not require the existence of closed-form random features: for computing the posterior mean, we can work with any kernel. We only introduce feature-based gradient estimators due to their use in previous literature (Lin et al. 2023 and Dai et. al. 2014). We ultimately recommend against their use for gradient estimation due to producing additive rather than multiplicative noise.

We do use random features for efficiently sampling from the GP prior, which is needed, for instance, in the Bayesian optimisation experiments. We similarly use random hashing for the Tanimoto kernel in the graph experiment. In non-stationary kernels, efficient sampling can often be done on a case-by-case basis, but we prefer not to focus on it as it is orthogonal to our contributions.

10. $\beta$ vs. $\beta n$

what is the difference between the $\beta$ in the main texts and the $\beta n$ in Algorithm 1?

The $\beta n = \beta \times n$ is simply the learning rate $\beta$ times $n$, the size of the training data set. This formulation allows us to keep the step-size $\beta$ times the data set size $n$ constant across data sets. We mostly use $\beta n$ throughout our writing, since the algorithm ultimately depends on the respective product.

11. Figure 3 green dashed line

The green dashed line is missing in Figure 3.

Thank you for pointing this out! We have updated the caption to reflect that this refers to the olive green line.

Summarising: thank you again for your feedback, which has made our experimental evaluation much more comprehensive. Please let us know what you think of the updated results!

4. SGD baseline divergence at large step-sizes

In Figure 1, we can not see the primal gradient descent becomes unstable and diverges for >0.1. Please show the unstable or compare the evaluated conditional numbers. Under higher step sizes, why does the gradient descent of the primal return nan (any possible reasons)?

We include the requested figures showing the divergence of primal SGD with high step-sizes, here and here, and Appendix C.6. This behaviour quickly exhausts the numerical floating point precision, leading to overflows and thus NaNs.

We did not include these in the original text because the loss curves diverge so quickly that they cannot be viewed simultaneusly with the optimisation trajectories of the non-diverging methods. We run our algorithm for 100k iterations. After only 10 iterations, amplifying oscillations are visible and after 20 iterations an exponential trend is already present.

5. Varying step-sizes

Figure 2 shows the random coordinate estimate with a step size equal to 50. what is the performance on varied step sizes?

We have conducted further experiments under the same settings with varying step-sizes for both random coordinate estimators (the one with additive and the one with multiplicative noise). Please see the corresponding figure here and more details in Appendix C.3.

The performance of the additive noise estimator plateaus earlier as the step-size increases. Smaller step-sizes converge slower but eventually reach better solutions, such that a trade-off is necessary: the step-size is proportional to the speed of convergence but inversely proportional to the final performance. In contrast, the multiplicative noise estimator consistently improves both convergence speed and final performance as the step-size is increased, as long as the step-size is not large enough to cause divergence.

6. Step size in Figure 3

What is the step size used to generate the Figure 3?

The step-size used to generate Figure 3 was $\beta n = 50$. We added this to the figure caption now.

7. Random feature plot line not matching in Figures 2,3 + rising blue line

It seems less than 50 and has a competitive rate compared to the random feature estimate shown in Figure 2. Extra clarification and comparison would be better..

Can any explanation of the rising part (the blue dashed line) in the middle figure in Figure 2 be given?

Thank you for bringing up this source of confusion. Only the dashed light blue line should be nearly identical in both cases. Both plots use a step-size of 50, but different batch sizes. However, we have discovered the lines are different due to the floating point precision used when running the experiment. Figure 2 was generated with 32-bit floating point, which becomes unstable when the error is very small (~ 1e-7). We have now replaced this with 64-bit floating point, so this is now fixed. Please see the fixed plot here. To improve clarity, we also added "random features" and "random coordinates" to the labels in Figure 2.

(continued below)

Thank you very much for your time reviewing our paper! We are very happy to hear our work was well-written, and focused on an interesting problem. Let us address the key questions:

1. Kernel and hyperparameter estimation.

It would be useful to emphasise that this work is useful when the Kernel is already known.

Thank you for pointing this out: We have further clarified this in our introduction with the language "As a result of the careful choices behind our method, we significantly improve upon the results of Lin et al. 2023 for stochastic gradient descent in the fixed kernel setting.".

Our work solely focuses on improving SGD inference and thus employs the setting of Lin et al., where the kernel is fixed.

Comments on whether these methods would be useful in hyperparameter estimation would be useful.

Our method can indeed be easily combined with standard approaches to kernel hyperparameter optimisation using the usual trick of applying the Hutchinson trace estimator on the marginal likelihood's gradient (see for instance Gardner et al. 2018). As a proof of concept, we have added a demonstration of this to Appendix C.7, and include some empirical results to that end where we match performance to standard optimisation of the marginal likelihood (see here).

2. Ease of implementation

The claim that the method can be implemented in a few lines of code should be demonstrated.

Thank you for your suggestion! We support our claim that SDD is easy to implement, by including a very simple Jupyter Notebook that implements our full method in NumPy from scratch and runs it on a toy problem. Without counting data set generation and kernel definition, a basic version of our algorithm can be implemented with only 20 lines of code.

3. Visualisation of samples

The paper would benefit from a visualisation comparing samples from a GP using SDD to an exact GP fit to show that the samples lie within the confidence interval.

We agree this would strengthen the paper! We have added a visual comparison to exact GP, and to our CG and SVGP baselines, on two toy examples considered in Lin et al. 2023. See the figure here and Appendix C.1.

4. Kernel of the form $\sum z_j z_j^T$

What are the implications of limiting the kernel to the form $\sum_{j=1}^m z_j z_j^T$?

Thank you for drawing us to this potential misunderstanding. We wish to clarify that our method does not require the kernel to be finite-dimensional. We only mention this family of approaches in the context of random feature gradient estimators. We study these due to their use in previous literature (Lin et al. 2023 and Dai et. al. 2014), but ultimately recommend against their use due to producing additive rather than multiplicative noise. To avoid confusion, we have added a line in Section 3.2: "For the sake of exposition (and exposition only, this is not an assumption of our algorithm), we assume that we are in a $m$-dimensional (finite) linear model setting,...".

5. Effect of ill-conditioning

How does ill conditioning affect the performance of the method?

Thank you for your question! In general, methods to solve linear systems are very dependent on the respective condition numbers. However, for Gaussian processes, Lin et al. 2023 show that SGD is significantly less sensitive to conditioning than CG, which can take an arbitrary amount of time to converge (and obtain good results) for sufficiently ill-conditioned systems.

For SGD, conditioning determines the maximum step-size which can be used without diverging. Since the dual gradient presents better conditioning, as it depends on the kernel matrix eigenvalues as opposed to their squares, we can use much larger step-sizes. This can be seen in Figure 1 in the main paper, and in a new figure which we added to the appendix here.

We further demonstrate this in Appendix C.5 where we compare both methods on the pol data set while setting the noise standard deviation to 0.001. The figure is also shown here.

Thank you for taking the time to read our work and providing feedback. We are delighted to hear that our paper was structured well and presented clearly and logically! Below we address the key concerns:

1. Geometric averaging

Notably, the authors introduce a new averaging scheme called geometric averaging.

We want to clarify that we did not introduce geometric averaging, which is widely used in deep learning. Our contribution is the observation that, for the estimators considered, geometric averaging outperforms arithmetic averaging due to the presence of multiplicative noise. This choice is unorthodox because almost all previous analysis of SGD uses gradient estimators with additive noise. In that setting, arithmetic averaging (as opposed to geometric) is optimal (Dieuleveut et. al. 2019).

2. Theoretical demonstration of convergence

The authors do not provide a theoretical justification to verify the convergence of the proposed method. Nevertheless, it is likely that convergence can be ensured under mild conditions, as the optimization techniques employed are standard and well-established in the community and literature.

Our paper is targeted at GP researchers and practitioners. While using SGD for GP inference has recently been found to be very promising (Lin et al. 2023 and Dai et. al. 2014), previous work makes a number of suboptimal choices when it comes to the details of the algorithm's formulation. Our contribution is handling these details carefully, highlighting to researchers and practitioners how to make SGD inference even more effective, particularly on harder problems like the graph kernel comparison.

In the full-batch (non stochastic) case, the analysis of the algorithm is standard (Nesterov 1983), and it is known that the dual objective will result in faster convergence than its primal counterpart due to conditioning. See, for instance, Section 3.7, page 289, of Bubeck 2015, a textbook on convex optimisation.

Our full algorithm makes use of an estimator with multiplicative noise. The convergence of Stochastic GD under multiplicative noise is currently at the forefront of theoretical research and beyond the scope of our more practical work. We have extended our discussion of pertinent work in this space in Section 3.4.

3. Significance

While the remaining studies and techniques are also important, they are somewhat incremental and standard. Consequently, I am uncertain about whether the paper's contribution alone justifies its publication in ICLR.

We agree that our paper uses insights from the theory literature. Indeed, this is highlighted throughout our text. To our knowledge, these insights have not been studied in a unified way by the kernel/GP community before. Even worse, all applied kernel/GP work we are familiar with, which use SGD methods, make suboptimal algorithmic design choices. On this basis we believe that a simple algorithm that combines recent theoretical advancements in a principled way, and is the first alternative which reliably outperforms conjugate gradient-based algorithms used in state-of-the-art Gaussian process software packages, can be of value to the ICLR community.

We would also like to note that in machine learning, there is a long history of combinations of small improvements eventually leading to breakthrough performance: the most prominent example, arguably, is the development of deep neural networks, which required small tweaks to initialisation, activation functions, architecture, optimisation, and other seemingly incremental factors in order to result in today’s effectiveness.

Thank you for your comment! We would like to clarify that "Nesterov" in Figure 7 refers to our proposed method. To prevent further confusion, we will change the corresponding label to "SDD". Further, we would like to reiterate that Figure 7 illustrates that our method reaches solutions with several magnitudes smaller squared $K$- and $K^2$-norms compared to AdaGrad, RMSprop and Adam.

Thank you for your review! We are delighted to see the significance of our graph neural network demonstration pointed out in your review! We would also like to highlight that our method reliably outperforms conjugate gradients, which is the standard optimisation method for this problem. We address your key concerns below.

1. Theoretical guarantees

"This paper only provides numerical experiments... would be good if rigorous theoretical guarantees could be proved, at least for some special cases."

Our paper is targeted at GP researchers and practitioners. While using SGD for GP inference has recently been found to be very promising (Lin et al. 2023 and Dai et. al. 2014), previous work makes a number of suboptimal choices when it comes to the details of the algorithm's formulation. Our contribution is handling these details carefully, highlighting to researchers and practitioners how to make SGD inference even more effective, particularly on harder problems like the graph kernel comparison.

In the full-batch (non-stochastic) case, the analysis of the algorithm is standard (Nesterov 1983), and it is known that the dual objective will result in faster convergence than its primal counterpart due to conditioning. See, for instance section 3.7, page 289, of Bubeck 2015, a textbook on convex optimization.

Our full algorithm makes use of an estimator with multiplicative noise. The convergence of Stochastic GD under multiplicative noise is currently at the forefront of theoretical research and beyond the scope of our more practical work. We have extended our discussion of pertinent work in this space in Section 3.4. We hope that our work's demonstration of the empirical performance achievable with SGD under multiplicative noise will motivate research in this area by the optimisation community.

2. Algorithmic details.

"The authors stress too much on the algorithm details. which can be deferred to the appendix for a major part of them while trying to leave some room for theoretical analysis."

Thank you for this comment. Our aim was to show how to make the already promising SGD-based learning algorithms for Gaussian processes more effective by giving a careful treatment to algorithmic details. To do so, it is very important to precisely state the algorithm we recommend, particularly for readers coming from the Gaussian process rather than optimisation community.

At the same time, we recognize that papers which focus more on the side of theoretical analysis are also valuable. We believe it is best to have both kinds of papers available in the literature, so that readers can more precisely find what they are looking for.

3. Comparison with other optimisers.

different optimizers such as AdaGrad, Adam, etc. The authors should also try these methods.

We now include an empirical comparison on the POL dataset of AdaGrad vs RMSprop vs Adam vs our method. We show that our method compares favourably to these methods by a significant margin: see linked figure and Appendix C.2.

Thank you for pointing this comparison out. We originally opted to omit these algorithms, as they are usually designed to tackle problems with a non-constant curvature. GP regression has a constant curvature: it yields a quadratic objective where Nesterov-type momentum is theoretically rate-optimal (Dieuleveut et. al. 2017 and Jain et. al. 2018). However, given the prevalence of other optimisation algorithms, we agree that a comparison to them is valuable to the larger community!

(continued in next comment)