Residual Deep Gaussian Processes on Manifolds

Provision of background on GPs in Riemannian geometry. (Weakness 1)

Thank you for raising this issue.

Indeed, it would be very helpful to provide the reader with an easy way to access a more complete introduction to GPs in Riemannian geometry. Fortunately, works that have initially developed these notions, and those that have later refined it, provide excellent introductions to the subject. Thus, instead of augmenting the appendix, we suggest that it is appropriate to instead direct the reader to these resources, which we now do in footnote 3.

Realising the code for public use.

Thank you for this suggestion.

We plan to release the updated code base, with the additional experiments done for the discussion phase, shortly after the discussion closes.

Motivation for GVF choice. (Weakness 1)

Thank you for bringing this issue to light.

We used the projected GVF construction in Section 4.4, since we believe that it is generally more useful for practitioners, as it does not require choosing a coordinate frame. Despite that, wanting to showcase both constructions, we included results from coordinate-frame GVFs in Section 4.2, where we found little difference between models in our preliminary testing.
We agree that our results would be more complete with both constructions in both experiments; however, we limited their scope to avoid overloading the presentation, especially since we did not want to emphasise the comparison between coordinate-frame and projected GVFs. Indeed, one of our key conclusions was to focus research on intrinsic approaches, such as the Hodge kernel, as they are superior not only in shallow GPs but also in deep models.

Parameter count matching. (Weakness 2)

Thank you for raising this issue.

In experiment 4.1 we attempt to make the comparison fair by matching the number of trainable parameters between models. We give the precise numbers of variational parameters and kernel parameters per level in paragraph "Variational Parameters" in the Appendix section A.3, which we have now expanded to include precise formulae for improved clarity. Additionally, we ensured that if we cannot match the number of parameters exactly, the baseline model would be given the highest number, to prevent artificial handicapping for our model.

Baseline model with additional layers could outperform residual deep GPs. (Weakness 2)

Thank you for this very apt observation.

To investigate the possibility of the baseline outperforming all manifold models, we extended experiment 4.1 for the baseline and Hodge models with up to 10 layers and doubled the number of randomised runs to 10 to increase confidence. We report the results in Figure 9 and Figure 10. We find that the baseline performance saturates more slowly giving the illusion that it might exceed that of the residual deep GP with deeper layers. However, what is especially clear with more data, it does saturate and never exceeds the performance of our model - neither in terms of NLPD nor MSE. In fact, with only three layers, our model surpasses the baseline across all depths.

It is true that the MSE advantage of our model is smaller than its NLPD advantage. However, it is important to stress that NLPD assesses both the predictive uncertainty and the predictive mean, making it a more complete metric for GP evaluation than MSE. Indeed, if we were prioritising MSE, then we probably would not use GPs at all.

Euclidean baselines in Sections 4.2 and 4.3. (Weakness 2)

We thank the reviewer for bringing forth this issue.

Firstly, however, it is important to mention that the results of these experiments are not in a vacuum with respect to Euclidean baselines, since they are both based on papers which show superior performance of shallow manifold GPs over Euclidean GPs. Nevertheless, we agree that a baseline of Euclidean deep GPs would bring additional context to Sections 4.2 and 4.3. Due to time restrictions we provided such a baseline only in experiment 4.3; however, we hope that the results will be of good help in assessing any advantage of residual deep GPs over Euclidean deep GPs overall.

We report a comparison of our model with the baseline in terms of NLPD and MSE in Figures 14 and 15 respectively. We find that our model achieves a lower NLPD than the baseline for all altitudes and model depths. Furthermore, the lowest mean MSE is attained by the residual deep GP for every altitude level; however, for 4 and 5 layers at altitudes of 2.0 km and 0.1 km, and 2 layers at altitude of 5.0 km, our model falls very slightly short of the baseline. Thus, as in Section 4.1, while Euclidean models may attain similar MSE to residual deep GPs, our models provide better-calibrated uncertainty estimates, resulting in better NLPD performance.

Dear Reviewer ezvQ,

Could you please confirm whether our rebuttal has adequately addressed your concerns?

We are happy to provide any further details.

Additional details on limitations of residual deep GPs

Thank you for this excellent suggestion.

It is true that in most of our experiments residual deep GPs saturate quickly with increasing depth, especially when data is scarce. We should stress, however, that this accords with the findings of Salimbeni and Deisenroth (2017) for Euclidean deep GPs, who additionally show that deeper layers become useful with much more data (up to a billion data points). We find an analogous trend of reduced saturation with increased data density in Figure 3; however, it would be useful to investigate whether this pattern holds for much larger real world datasets. We now highlight this limitation as a potential research direction in the Conclusion section.

When to choose residual deep GPs over manifold neural networks?

Thank you for this question.

Unfortunately, the literature on manifold neural networks is scarce and practically non-existent for their probabilistic formulations. Thus, any judgement of preference for one over the other we would make would be chiefly speculative. A comparison between these two frameworks would be certainly valuable; however, due to the additional work necessary on the side of neural networks, we believe it is better suited for future research.

Thank you for your responses. Regarding the second point, one simple practical version of probabilistic manifold neural networks would be using deep ensembles, but I agree with the authors that additional work on the manifold neural network side would be necessary for a meaningful comparison.

I would encourage the authors to include the discussions in their rebuttal regarding limitations of residual deep GPs and perhaps the future research of manifold neural networks in the revision of the paper.

Benefits of residual deep GPs without a known manifold. (Question 1)

Thank you for this question.

Our results do show that residual deep GPs are effective models for manifold-input data, outperforming shallow manifold GPs and Euclidean deep GPs with a geometry-aware input layer. Moreover, we do suggest that, if an auxiliary manifold that does not degrade prediction quality is found, our model could serve as a faster alternative to Euclidean deep GPs. However, we do not intend to imply that our model ought to be used where the manifold of the data is unknown---as kindly remarked by reviewer a17f, we deliberately make this limitation known in Section 4.4. If there is any part of the text that is misleading, we are eager to correct it.

Overcoming training issues and questionable uncertainty estimates. (Question 1)

Thank you for this question.

Our results suggest that training and uncertainty estimation are not an issue---at least insofar as they are not an issue for shallow GPs. Indeed, all experiments show that our models outperform shallow GPs in MSE and, most importantly, NLPD. In particular, in Section 4.1 we find that initialising the variance of GVFs to a small value, e.g. 0.0001, along with the residual structure, helps the model recover a shallow solution where data is scarce. Additionally, in Section 4.3, we show that this improvement in metrics translates to an improvement in visual interpretability of predictions. Indeed, we see that predictive uncertainty of residual deep GPs captures information not only about how far a test point is from a training point (as in the case of shallow GPs) but also about the behaviour of the nearby training values---e.g., regions with large differences between nearby wind velocities have high predictive uncertainty.

Motivation for modelling on the hyper-sphere. (Question 2)

Thank you for this question.

The choice of hypersphere as the auxiliary manifold in Section 4.4 is motivated by the work of Dutordoir et al. (2020). To justify this choice, they note that GPs representing limits of neural networks have kernels which normalise their inputs, and thus can be seen as kernels on hyperspheres.

We focus on hyperspheres in our work as a whole because: 1. they are spaces of practical relevance to many real-world problems; 2. they make the application of manifold-specific techniques quite simple. These two advantages were crucial for us, since our goal was to evaluate the real-world relevance of our model and examine how it fits in with existing manifold GP machinery. Nevertheless, our choice to focus on hyperspheres does not limit the applicability of our model to other manifolds.

Combining known manifolds in composition. (Question 3)

Thank you for this interesting suggestion.

Indeed, if we understand correctly, it is suggested to use different manifolds at subsequent layers. This is a very interesting problem---one that we have thought about, but is beyond the scope of our current paper. One can immediately imagine constructing such sequences of manifolds through inclusion, projection, or other more complicated maps; however, we believe this would require extensive exploration and additional theoretical consideration better suited for a follow-up work.

Relation of residual deep GPs to wrapped GPs. (Question 4)

Thank you for suggesting that this topic should be expanded upon. We hope to have addressed your question in the general response.

Expand caption to Figure 2. (Question 5)

Thank you for this suggestion.

Adding an in-caption explanation of the GVF constructions would certainly help avoid jumping between pages. We will attempt to make space for and expand the caption once all other changes are finalised at the end of the discussion process.

How to choose when to switch the models and is the 180th iteration optimal? (Weakness 2)

Thank you for this question.

With the switching strategy, we aimed to show that by only using a deep GP when the shallow GP gets stuck, we can achieve large performance improvements without a runtime penalty. Indeed, the choice of the 180th optimisation step is not optimal in terms of attaining the lowest regret at step 200 - instead, it was chosen so that it is evident from the flat regret curve that the shallow GP is stuck by the time we switch to a deep model. Choosing when to switch is non-trivial and largely depends on the compute budget, since, as seen below, we can often achieve better results by employing only deep GPs. What our results do suggest is that an adaptive strategy, such as switching periodically to a deep GP when the shallow GP is stuck, could achieve a good balance of compute to performance.

Using composite GPs from iteration 1. (Weakness 2)

Thank you very much for this suggestion.

Employing residual deep GPs from the very beginning is certainly possible, since at low data densities the model recovers a shallow GP. The only difference is that it would recover a sparse shallow GP, whereas in Section 4.2 we use an exact shallow GP. To see the actual behaviour, we performed 15 deep-GP-only BO runs with the irregular target function. Through this process we found that the combination of hyper-priors and quasi-Newton optimisation employed in Jaquier et al. (2022) is useful for effective exploration with the shallow GP and, indeed, with the deep GP too (we remark about this behaviour in Appendix A.4).

We report the results of this experiment in Figure 12. In 12 runs, the deep GP found a point with significantly lower regret than shallow GPs - often in fewer than 100 iterations. In 3 runs, however, the deep GP got stuck in a local minimum. The reasons for this are not fully clear. It may be caused by the approximate nature of the shallow GP recovered by the deep GP or by using the L-BFGS optimiser instead of BFGS as for the exact shallow GP - regardless, additional experiments would be needed to explain this behaviour. Despite this, if sample-efficiency is key, Bayesian optimisation with deep GP from the start can be more effective than the synthetic strategy beginning with a shallow GP.

Switching to a different shallow GP. (Weakness 2)

Thank you for this question.

To find out, we performed 10 randomised BO runs of the irregular function, switching to a shallow GP with the manifold squared-exponential kernel at iteration 180. We found no improvement in the quality of minima found, despite the change of kernel.

Hyper-parameters and hyper-priors. (Weakness 2)

Thank you for this question.

The hyper-priors, which we set following Jaquier et al. (2022), are only used in the training of the models and do not play a role at inference time. In particular, we do not marginalise the hyper-parameters.

Increased variance of regret after switching to a deep model. (Weakness 2)

Thank you for this astute observation.

Fortunately, the reasons for this increase in the variance of regret after switching to residual deep GPs can be easily elucidated, which we now do in the paragraph "Additional results and analysis" of Appendix A.4. This increase in variance is due to the fact that there is a large variation in the initial 180 points acquired by the shallow GP between different runs. The locations of these points influence the quality of fit of the deep GP and, if the points were particularly close and dense around the optimum, the deep model often improves after 1 or 2 iterations, but with fewer points around the optimum we find the fit poorer and the deep model requires a higher number of iterations and, although rarely, sometimes 20 is not enough. This variance in the number of steps needed after the transition to find a better optimum leads to the increased variance in regret.

Additional benchmark functions. (Weakness 2)

Thank you for this suggestion.

We agree that additional benchmark functions would help examine the behaviour of Bayesian optimisation with our model more broadly. In fact, it would be very interesting to dive deeper into geometry-aware Bayesian optimisation with residual deep GPs, the different strategies for their use with shallow GPs, and their use in various optimisation landscapes; however, we believe that a more thorough exploration is better suited for a follow-up work. Thus, we decided to limit the number of examined functions, balancing the sizes of different experiments for a more well-rounded examination of our model.

Comparison to wrapped GP regression. (Weakness 1 and Question 1)

Thank you for suggesting that this matter be discussed in more detail. We hope to have answered your question in the general response.

Why do residual deep GPs not improve upon shallow GPs in BO on Ackley. (Weakness 2)

Thank you for bringing this point to our attention.

Indeed, because the Ackley function exhibits significant non-stationarity, residual deep GPs do model it better in regions where enough data is present to capture the function's complexity. In particular, with a grid of points on the sphere, the residual deep GP we used models the entire function better than the shallow GP, both in terms of MSE and NLPD as seen in the Table below for 100, 200, and 400 points. As in the original experiment, we set the inducing locations to the training data locations.

Yet, in spite of these results, the lack of improvement of the deep model over the shallow GP in Bayesian optimisation is not an anomaly. This is because typically at iteration 180 the shallow model has already found the true minimum mode; thus, to make any improvements, it is only necessary to model the function in a region around that mode. Since in that region the function is smooth, the shallow GP is as effective as the deep GP at modelling it. We designed the irregular function to showcase what can happen when a minimum is not in a smooth region but rather is approached by nearing a singularity of the objective function. In this case, when the models are switched the shallow model has also typically found the mode but is unable to accurately model the region near it due to the singularity---here the deep GP can yield an improvement.

We hope that the modified phrasing of the last sentence in the "Results" paragraph of Section 4.2 conveys the explanation above.

Optimiser choice. (Weakness 3)

Thank you for bringing this issue to light.

Firstly, it is important to stress that both our model and the Euclidean baseline are trained using the Adam optimiser, making the comparison between them valid in its own right.

Nevertheless, our choice of optimiser is not arbitrary. Chiefly, we choose Adam as it alleviates memory issues via mini-batching, while the L-BFGS implementation in large libraries such as Scipy or Optax are classical and do not support stochastic optimisation---unlike some of the recent L-BFGS variants (Byrd et al. 2014, Moritz et al. 2016). Thus, Adam is more convenient both for us and for anyone drawing on our results in the future. Secondly, Adam is used by Salimbeni and Deisenroth (2017), making our results more readily comparable with their experiments, which gives useful context to our results, since residual deep GPs generalise their linear mean model. It is true that matching our optimiser with Salimbeni and Deisenroth (2017) causes a mismatch with the work of Dutordoir et al. (2020). However, we believe that the preceding advantages to outweigh this disadvantage, especially since, even with L-BFGS, we would not be able to replicate the experimental setup of Dutordoir et al. (2020) because the collapsed ELBO objective they use is intractable for deep GPs.

Wall-clock time to a given accuracy. (Weakness 3)

Thank you for this suggestion.

We believe that, if our model could achieve the accuracy of the Euclidean model, a wall-clock time to accuracy metric would be most appropriate. However, we chose the number of training iterations so that both models approximately reach convergence; thus, in this setting residual deep GPs cannot achieve the performance of Euclidean deep GPs---no matter the wall-clock time. For this reason we choose to report per-iteration wall time and performance separately. Indeed, although we report a negative result, with this separation we hope to make clear that, if the performance improved for a different (compact) auxiliary manifold, we could have an alternative to Euclidean deep GPs that is about twice as fast.

We have clarified the reason behind the number of training iterations in Appendix A.

Implications for vector- and function-valued mappings. (Question 2)

Thank you for this interesting question.

For finite-dimensional vector spaces, we have an answer. Since such spaces are isomorphic to the Euclidean manifold of the same dimension, our model reduces to the linear mean deep GP of Salimbeni and Deisenroth (2017). For function spaces which are compact or non-compact symmetric Riemannian manifolds we can apply our model directly, though non-compact manifolds, e.g. that of positive semi-definite bilinear forms, require sampling from the spectral density of the kernel when approximating its series expression (Azangulov and Smolensky et al. 2023). For infinite-dimensional function spaces, this might be an interesting research direction.