# Explainable Gaussian processes: a loss landscape perspective

Maximilian P Niroomand, Luke Dicks, Edward O Pyzer-Knapp, David J Wales

<figure><img src="image_1.png" /></figure>

## To cite this version:

Maximilian P Niroomand, Luke Dicks, Edward O Pyzer-Knapp, David J Wales. Explainable Gaussian processes: a loss landscape perspective. Machine Learning: Science and Technology, 2023, 5, 10.1088/2632-2153/ad62ad. hal-04838780

HAL Id: hal-04838780

https://hal.science/hal-04838780v1

Submitted on 15 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

PAPER • OPEN ACCESS

# Explainable Gaussian processes: a loss landscape perspective

To cite this article: Maximilian P Niroomand et al 2024 Mach. Learn.: Sci. Technol. 5 035025

View the article online for updates and enhancements.

You may also like

- Faster variational quantum algorithms with quantum kernel-based surrogate models
Alistair W R Smith, A J Paige and M S Kim

- Data efficiency and extrapolation trends in neural network interatomic potentials
Joshua A Vita and Daniel Schwalbe-Koda

- Neural network Gaussian processes as efficient models of potential energy surfaces for polyatomic molecules
J Dai and R V Krems

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

21 November 2023

REVISED

14 June 2024

ACCEPTED FOR PUBLICATION

12 July 2024

PUBLISHED

23 July 2024

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_6.png" /></figure>

PAPER

# Explainable Gaussian processes: a loss landscape perspective

Maximilian P Niroomand¹,*, Luke Dicks² ID, Edward O Pyzer-Knapp² and David J Wales¹ ID

¹ University of Cambridge, Cambridge, United Kingdom

² IBM Research Europe, Daresbury, United Kingdom

* Author to whom any correspondence should be addressed.

E-mail: mpn26@cam.ac.uk

Keywords: loss landscapes, Gaussian processes, machine learning, physics-inspired

## Abstract

Prior beliefs about the latent function to shape inductive biases can be incorporated into a Gaussian process (GP) via the kernel. However, beyond kernel choices, the decision-making process of GP models remains poorly understood. In this work, we contribute an analysis of the loss landscape for GP models using methods from chemical physics. We demonstrate $v$-continuity for Matérn kernels and outline aspects of catastrophe theory at critical points in the loss landscape. By directly including $v$ in the hyperparameter optimisation for Matérn kernels, we find that typical values of $v$ can be far from optimal in terms of performance. We also provide an *a priori* method for evaluating the effect of GP ensembles and discuss various voting approaches based on physical properties of the loss landscape. The utility of these approaches is demonstrated for various synthetic and real datasets. Our findings provide insight into hyperparameter optimisation for GPs and offer practical guidance for improving their performance and interpretability in a range of applications.

## 1. Introduction

The interpretation and understanding of model decision-making is one of the most important unsolved problems in machine learning. Practitioners in fields from healthcare to cybersecurity remain reluctant to deploy machine learning solutions because the decision process remains opaque.

Gaussian Processes (GPs) are a well-known class of nonparametric, supervised learning models that introduce a Bayesian framework to machine learning. The Bayesian framework allows the construction of a confidence measure around predictions, and therefore GPs overcome one of the most prominent limitations of other, non-Bayesian methods, namely uncertainty quantification. Uncertainty in this context is a general term and below we specifically refer to the variance of a prediction. While uncertainty quantification provides an important step towards model reliability, the more fundamental issue of interpretability remains unresolved, and a significant barrier to the widespread adoption of GPs in critical and sensitive applications.

Mathematically, a GP is a collection of random variables, any finite subset of which is jointly Gaussian distributed [1]. Unlike a neural network, where the learned parameters predict some output in combination with the functional form of the model, a GP learns a distribution over functions to fit the data. In parametric machine learning, a distribution $p(y|X, \theta)$ is employed to explain some input data $X$, with the optimal set of coefficients $\theta$ inferred by maximum likelihood or maximum *a posteriori* estimation. In nonparametric methods, such as GPs, a weighted distribution over functions $p(f|X) = N(f|\mu, \sigma)$ is learned, where the parameters are the function itself. Unlike parametric methods, updating the GP with new information directly changes the shape of the function, not parameters that specify the function. Given some data, a posterior distribution $p(f_*|X_*, X, f)$ for new predictions $f_*$ at position $X_*$ can be computed.

GPs allow the incorporation of inductive biases via the covariance kernel. The choice of kernel, together with a mean function $\mu(x)$ (which is usually just set to 0 such that $\mu(x) = 0$), is the critical design choice when constructing a GP. However, the number of kernels being used in practice is surprisingly small [2], even though methodology for more elaborate kernel selection exists [3]. Driven partly by availability in standard GP Python packages, but also by convention, only a handful of different kernels are commonly used. In this contribution we explore alternative kernels and provide insights into the benefits and
disadvantages of using non-standard choices. In particular, we study the Matérn kernel, a general formulation that encompasses most other popular kernels by variation of its smoothness parameter, $\nu$.

Obtaining a variance estimate in GPs does not come for free and the most prominent disadvantage of using GPs over parametric machine learning methods is computational cost. The training of GPs scales naively as $O(n^3)$ for $n$ data points due to inversion of the covariance matrix. Various approaches exist to reduce the complexity [4], which usually revolve around the concept of sparse GPs [5] and inducing point methods [6] or, more recently, nearest-neighbour approaches [7]. However, with computational resources continuing to grow, and an increased demand for explainable machine learning, the growing interest in GPs will likely continue.

The study of loss landscapes in machine learning aims to increase interpretability of black box models and identify new ways to think about the solutions. Due to the computational cost associated with evaluating large areas of the loss landscape, relatively little work exists on characterisation [8] and visualisation [9, 10]. Here, we study loss landscapes using tools from the physical sciences [11], where molecular energy landscapes are employed to predict thermodynamic and kinetic properties of a system. We will treat ML loss landscapes analogously, as in previous work [12–14], which identified various similarities between molecular energy landscapes and machine learning loss landscapes; a tutorial review on these applications is available in [15].

In GPs, the loss function is usually the negative log marginal likelihood (NLML). Since GPs are a non-parametric method, the training phase does not include learning a set of model parameters, in contrast to neural networks. Rather, the solution space of the loss function and hence its arguments, is the set of hyperparameters $\theta \in \Theta$ that characterise the kernel. Current research in loss landscapes focuses almost exclusively on supervised parametric methods. Such landscapes usually exhibit rich and high-dimensional structures due to the vast number of parameters. Loss landscapes of non-parametric methods, such as GPs, are inherently simpler, due to a substantially smaller number of parameters. In a GP loss landscape, the global minimum represents the set of kernel hyperparameters, $\theta$, that provide the most accurate solution given training data, and locating this solution is the principal challenge for GP fitting. Understanding the solution space for a GP therefore provides an understanding of robustness and inference accuracy of the underlying GP model.

## 1.1. Related work

**Loss landscapes:** The study of loss landscapes in machine learning has received increasing attention in recent years. Insights into the structure and geometric properties of the loss landscape have been employed to improve model performance [13, 16], adversarial robustness [17], interpretability [18, 19], and generalisability [20–22] of neural networks. An overview focusing on the role of loss landscapes in machine learning, and the link to the physical sciences, is given in [15]. Empirically, loss landscape characteristics have also been considered to explain fundamental aspects of machine learning, such as the quality of minima [23] and the structure of the solution space as a function of hyperparameters [12]. Unlike parametric methods, the parameters in GPs are the function itself [1] and only a selected few hyperparameters specifying the kernel must be learned. Here, this problem amounts to minimising the NLML surface with respect to the hyperparameters. Previous work has considered loss landscapes in the context of Bayesian optimisation. However, we believe that the present contribution is the first account of a detailed study of GP loss landscapes.

**GP interpretability:** The interpretation of machine learning models is one of the most active fields of research today [24]. In GPs specifically, various approaches have been employed to interpret inference beyond kernel choices. Interpretability has several facets that stem from hyperparameter optimisation, such as feature importance, robustness of the hyperparameter choices, and how difficult optimal values are to locate. Previous approaches include kernel decomposition methods and local explanations to mirror linear models [25, 26], while others employ sensitivity analysis, a more classical approach used across various machine learning methods [27]. Some kernels, referred to as automatic relevance detection (ARD) kernels, may also be more readily interpretable than others. ARD kernels are anisotropic, meaning that they have a dedicated length scale parameter for each dimension of the input data. Studying the individual length scales may then be used for feature selection, as for example in [28]. Lastly, additive GPs [29] can be used for similar purposes, decomposing the latent function into a sum of simpler functions to quantify contributions from individual input variables. Yet, to the best of our knowledge, no methods exist that consider the loss landscape towards GP interpretability.

**GP ensemble methods:** Ensemble methods that combine multiple predictors by methods such as bagging [30], are commonly used to improve model performance. The intuition that a combination of single predictors should outperform a single model has been confirmed theoretically and empirically [13, 31]. Both papers ascribe the performance increase in bagging ensembles to a set of qualitatively different modes of the
loss landscape, akin to different minima. GP ensembles have been used in online learning applications [32], usually with Bayesian ensemble weights [33], and as GP experts with different kernel functions [34]. However, widespread adoption of GP ensembles is hindered by high computational costs and often only minor improvements in terms of inference accuracy. In contrast to neural network ensembles, none of the above methods consider features of the loss landscapes, as, for example, suggested in Fort, Hu, and Lakshminarayanan [16] or [13]. In these reports, loss landscape features are employed to determine whether the substantial cost of using ensembles is justified. Hence, a method to determine ensemble performance a priori is highly desirable.

**Matérn kernel:** The mathematical form of the Matérn kernel is introduced in the following section. Usually, a standard implementation of a half-integer version of the kernel is employed by machine learning practitioners [35]. A helpful account of covariance kernels in general can be found in [36]. For their flexibility, Matérn kernels have also been studied in the context of non-Euclidean, Riemannian geometries [37] and are commonly employed for physical applications due to their reduced smoothness compared to other kernels [38, 39].

## 1.2. Contribution and organisation

In this work, we apply state-of-the-art methods from chemical physics to improve the understanding and interpretability of GPs. To the best of our knowledge, we are the first to study and visualise loss landscapes of the NLML function for GPs. Specifically, we outline a novel approach to making decisions in GP model design, such as the kernel choice, using methods from the energy landscapes field. Our interest lies in the analysis of loss surface topology for various standard GP problems and its direct applications towards improving the usability of GPs. In summary, we

* Outline how methods and knowledge from chemical physics can be employed in interpretability for Bayesian machine learning.

* Provide novel insights into non-standard kernels by directly including $\nu$ in the optimisation problem and explore the evolving hyperparameter optimisation problem.

* Describe the use of geometric and physical features of the loss landscapes to increase GP ensemble performance and reduce computation time.

# 2. Methods

## 2.1. GP

Let $f: \mathcal{X} \to \mathbb{R}$ be a random function on some data $\mathcal{X}$. A GP is a collection of random variables $f(x)$, $x \in \mathcal{X}^d$ such that for any finite collection of points, and any $d$, the joint distribution of variables in the vector $f(\mathbf{x})$ is a multivariate Gaussian distributed with prior mean vector $\mu$ and covariance matrix $\Sigma$. A latent function, $f$, can hence be assigned a GP prior $f \sim GP(\mu(x), k(x, x'))$, with $\mu(x) = 0$ being a common choice for the mean function and $k(x, x')$ being some covariance kernel. Thus, the latent function in a GP is represented by samples from an infinite multivariate Gaussian distribution, parameterised by a mean and covariance function. Given the GP prior $f$, together with some training data $\mathcal{D} = \{\mathcal{X}, y\}$, and allowing for noise by setting $y_i = f(x_i) + \epsilon_i$, $\epsilon_i \sim N(0, \sigma^2)$, the posterior $p(f|\mathbf{y})$ is conveniently also Gaussian distributed with

$$ \mu(\mathbf{x}_*) = \mathbf{\Sigma}_{\mathbf{x}_*, \mathbf{x}} (\mathbf{\Sigma}_{\mathbf{x}, \mathbf{x}} + \sigma^2 \mathbf{I})^{-1} \mathbf{y} \\ \sigma^2(\mathbf{x}_*) = \mathbf{\Sigma}_{\mathbf{x}_*, \mathbf{x}_*} - \mathbf{\Sigma}_{\mathbf{x}_*, \mathbf{x}} (\mathbf{\Sigma}_{\mathbf{x}, \mathbf{x}} + \sigma^2 \mathbf{I})^{-1} \mathbf{\Sigma}_{\mathbf{x}, \mathbf{x}_*} \quad (1) $$

for some new testing point $\mathbf{x}_*$, where $I$ is the identity matrix. By virtue of the Bayesian approach, the posterior can be updated for incoming data, with the values of $\mu(\mathbf{x}_*)$ specifying the prediction and $\sigma^2(\mathbf{x}_*)$ the variance at point $\mathbf{x}_*$. In training a GP, we commonly minimise the NLML function

$$ \log p(y|\mathcal{X}, \theta) = -\frac{1}{2} y^T \Sigma^{-1} y - \frac{1}{2} \log |\Sigma| - \frac{n}{2} \log 2\pi, \quad (2) $$

with respect to a set of hyperparameters $\theta$, and the hyperparameters that globally minimise this function provide the optimal GP fit. The NLML surface can contain many local minima and the key challenge to GP fitting is location of its global minimum. The NLML is a suitable objective function, since it contains both a model complexity term and a data fit term which can counterbalance for a well-performing model.

## 2.2. Matérn kernel

One of the most important decisions when designing a GP is the choice of kernel or covariance function, which computes the correlation (covariance) between any pair of points in the training dataset $\mathcal{D}$. Element
$k_{ij}$ of a covariance matrix $\Sigma$ contains the covariance between points i and j. The general Matérn kernel is defined as

$$k_{ij} = k(x_i, x_j) = \sigma^2 \frac{1}{\Gamma(\nu) 2^{\nu-1}} \left( \frac{\sqrt{2\nu}}{\ell} d(x_i, x_j) \right)^{\nu} K_{\nu} \left( \frac{\sqrt{2\nu}}{\ell} d(x_i, x_j) \right) \quad (3)$$

where $\Gamma(\cdot)$ is the Gamma function, $K_\nu$ is the modified Bessel function of the second kind of order $\nu$ and $d(\cdot, \cdot)$ is the Euclidean distance between two points. The kernel is parameterised by an amplitude $\sigma$ and a lengthscale $\ell^d$, which for the anisotropic kernel we use here is a vector of length $d$ for input data $x \in \mathbb{R}^d$.

While the Matérn function can theoretically be evaluated for all values of $\nu$, for which it is $[\nu] - 1$ differentiable, there exist special solutions at half integers. The commonly used 1/2, 3/2 and 5/2 Matérn kernels are realised as

$$M_{1/2}(d) = \sigma^2 \exp\left(-\frac{d(x_i, x_j)}{\ell}\right), \\ M_{3/2}(d) = \sigma^2 \left(1 + \frac{\sqrt{3}d(x_i, x_j)}{\ell}\right) \exp\left(-\frac{\sqrt{3}d}{\ell}\right), \\ M_{5/2}(d) = \sigma^2 \left(1 + \frac{\sqrt{5}d(x_i, x_j)}{\ell} + \frac{5d(x_i, x_j)^2}{3\ell^2}\right) \exp\left(-\frac{\sqrt{5}d(x_i, x_j)}{\ell}\right). \quad (4)$$

The Matérn $M_{1/2}$ kernel is identical to the absolute exponential kernel. The smoothness of the Matérn kernel increases with $\nu$, and for $\nu \to \infty$, it becomes equivalent to the RBF kernel. Hence many of the most commonly used kernels are directly incorporated into the Matérn kernel when we consider wide ranges of $\nu$. In practice, low values of $\nu$ may be chosen because the reduced smoothness can fit complex data more accurately. Derivatives of the Matérn kernel with respect to its smoothness parameter $\nu$ require derivatives of the modified Bessel function of the second kind $K_\nu$ w.r.t. $\nu$, for which numerical derivatives are usually employed. To increase efficiency and accuracy of these derivative computations, Geoga et al [40] have recently presented an implementation using automatic differentiation. However, computing non-half integer kernels remains expensive and a bottleneck in motivating the use of such values of $\nu$.

## 2.3. Loss landscape exploration

To characterise a GP loss landscape, an efficient way to explore this complex, high-dimensional functional surface is required. This challenge is analogous to the characterisation of solution spaces for molecules in chemical physics. The energy landscape framework [11] can be used to efficiently locate stable molecular configurations, explore the transitions between them, and visualise potential energy surfaces. Characterising minima and transition pathways is critical in exploring and characterising loss landscapes, and the algorithms adapted from the energy landscape framework to achieve this goal are introduced below.

We represent a loss landscape by its stationary points, where the gradient vanishes, specifically local minima and the transition states that connect them. Transition states are maximal in one direction and minimal in all orthogonal directions, i.e. they exhibit a single negative eigenvalue of the second-derivative (Hessian) matrix. In physical systems transition states characterise energy barriers between the two minima that are connected by steepest-descent paths parallel and antiparallel to the Hessian eigenvector corresponding to the unique negative eigenvalue. For loss function surfaces they report on the intermediate topology between local minima in the same way.

Local minima are explored using the basin-hopping algorithm [41–43]. Basin-hopping global optimisation involves taking steps to a different region of the landscape by perturbing the current local minimum in a Markov chain. After local minimisation, this step is accepted or rejected using a Metropolis criterion of $P \propto \exp(-\Delta\mathcal{L}/c)$ for some constant $c$ and the difference in loss value between the old minimum and the new minimum $\Delta\mathcal{L}$. If accepted, the new local minimum becomes the current minimum and the process continues. Due to the theoretical unboundedness of the GP loss landscape space, the algorithm stops after a convergence criterion is reached. Such a criterion may be a fixed number of basin-hopping steps, or more commonly, be defined by some interval $n$ for which no new lower minimum was identified in the last $n$ basin-hopping steps. By iteratively performing this search, multiple low-lying minima of the loss landscape are obtained.

Various methods for finding transition states are discussed in [11] including the doubly-nudged elastic band and hybrid eigenvector-following approaches used in the present work. We find transition states between all minima located via basin-hopping and the minima they connect. From these stationary points we can visualise the surface topography using disconnectivity graphs [44, 45]. Disconnectivity graphs
provide a coarse-grained, low-dimensional representation of the landscape, characterised by the minima and transition states. Each terminal leaf node in a disconnectivity graph is a minimum, and the splitting of parent nodes is defined by the transition states between minima. The vertical axis represents the loss value, here the NLML. The order of nodes on the horizontal axis is determined by the splitting of daughter nodes from the parent node in the level above, putting nodes that contain the lowest minima in the middle. The nodes are spaced on the horizontal axis to allow for all the subsequent branches as the loss value decreases.

## 3. Results and discussion

To show the value of characterising the loss landscape using methods from chemical physics, we present selected examples of how these methods can improve our understanding of GPs. The examples are certainly not exhaustive, and we note that an extension to other kernels, to which these methods generalise without any restrictions, is a natural avenue for future work. Here, we present results that focus on hyperparameter choices for the Matérn kernel and loss landscape guided GP ensembles.

### 3.1. Matérn kernel for changing $\nu$

In this section we analyse the changing topology of the loss landscape when modifying the Matérn $\nu$ parameter. We constructed the loss landscape for each $\nu$ in the interval [0.5, 4.0; 0.1] using the methods described in section 2.3. The loss landscape at each $\nu$ is specified over all remaining kernel hyperparameters. We selected the range [0.5, 4.0] since it encompasses all the commonly used Matérn kernels $\nu \in \{0.5, 1.5, 2.5\}$ and also explores the regime corresponding to larger $\nu$ values. The transition between different $\nu$ values permits the exploration of the marginal likelihood across the space of common kernels. In practice, practitioners routinely select a single popular kernel irrespective of the given application, but here we examine the loss landscape changes that occur due to kernel specification. We performed this analysis for the $d$-dimensional Schwefel function defined by

$$f(x) = 418.9829d - \sum_{i=1}^{d} x_i \sin \sqrt{|x_i|}, \quad (5)$$

for which we constructed a dataset composed of 400 data points in the [−100, 100] hypercube of a 3d Schwefel function. The Schwefel function is commonly used to study GP optimisation problems, and provides a way to test GPs on a relatively complex, but well-defined problem.

In figure 1, we compare the loss values of all minima, located via basin-hopping, to the corresponding mean squared error (MSE). The emerging pattern is identical for the two plots, and the sequential changes between different values of $\nu$ are clear. Minor differences are observed, which is expected, since loss does not necessarily correlate perfectly with test performance, due to differences in the distribution of points in training and testing datasets. We observe a continuous, smooth change in the loss value of minima as a function of increasing $\nu$. The clustering highlights which minima are qualitatively the same, and are preserved across all landscapes. Most minima are present for almost the entire $\nu$-range. The continuous, smooth changes in the loss landscape with changing $\nu$ are characterised by smooth changes of minima.

This view is strongly reinforced by the disconnectivity graphs in figure 2, which provide an alternative visualisation of the landscapes. Disconnectivity graphs allow us to faithfully represent the barriers in addition to the minima, and reveal the global organisation of the surface. The complete set of landscapes at each $\nu$ are given in appendix A, and here we highlight some key examples. The same minima are displayed as in figure 1, but now with additional connectivity information.

The loss landscapes for sequential values of $\nu$, visualised as disconnectivity graphs, underline the $\nu$-dependent continuity. Changes in the landscapes are largely continuous and smooth, with the exception of minima that disappear at $\nu \le 2.1$ and $\nu \le 3.5$. In both cases, an additional minimum exists for higher $\nu$ but disappears when $\nu$ decreases below the threshold.

As $\nu$ decreases towards the critical points at which we observe a disappearance of minima, the barriers separating the disappearing minimum with an adjacent one decrease smoothly. This behaviour occurs through a fold catastrophe, where a transition state and minimum coalesce. Catastrophe theory provides the mathematical machinery to describe the families of landscapes with different topologies that appear for continuous parameter changes, as for $\nu$ [46, 47]. Previously, Leonardos and Piliouras [48] have considered fold catastrophes for machine learning loss landscapes in Q-learning, which we extend here to nonparametric machine learning and GPs. These non-smooth, radical transitions between subsequent values of $\nu$ highlight the importance of hyperparameter tuning in GPs, and specifically the importance of non-standard (half-integer) values of $\nu$. Minima that disappear in a fold catastrophe are a concern to the practitioner who may observe large discontinuous changes in model performance for small perturbations in
<figure><img src="image_7.png" /><figcaption>**Figure 1.** (Left) Loss value (NLML) and (right) mean squared error for all identified minima at all values of *v*. *k*-means clustering is performed on the minima, where *k* is chosen to be the maximum number of minima at any *v*. The data shown here were obtained from the 3d Schwefel function.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 2.** Loss landscape fold catastrophes illustrated in disconnectivity graphs. Each leaf node of the graphs is a minimum on the NLML surface. The disappearance of minima from (a) $v = 2.1 \rightarrow v = 2.0$ and (b) $v = 3.5 \rightarrow v = 3.4$ corresponds to a reduction in the number of leaf nodes.</figcaption></figure>

*v*. Additional experiments to confirm that, outside the fold catastrophes highlighted above, individual values of *v* do have smoothly changing solution spaces are included in appendix A.

## 3.2. *v* in hyperparameter optimisation

We have analysed the effect of the Matérn *v* on loss landscapes, and the corresponding sharp changes in GP accuracy, through minima addition, removal, or reranking. Here, rather than computing the landscape for a fixed *v*, we include it in the hyperparameter optimisation. Instead of fixing *v* and optimising the NLML surface for the remaining hyperparameters, we optimise *v* as part of *θ*. The lack of previous work on this idea can perhaps largely be attributed to the increase in computational cost of computing non-half-integer solutions and the absence of this method in standard Python packages. One of the few works that study non-half-integer values of *v* is the study of the bandit problem in reproducing kernel Hilbert spaces by Janz et al [49].

For the 3-dimensional Schwefel function and a varying number of data points with repeated sampling, we obtain the results shown in table 1. In general, the best test MSE improves substantially with an increased number of training points, whilst the training MSE is highest for a small number of datapoints, clearly indicating overfitting in the low-data regimes. The main result, however, is that the best value for *v* stabilises for increasing *n* in the region of *v* = 4.6. Importantly, when we compare the *n* = 900 case with a standard *v* = 2.5 kernel, we find that including *v* in the optimisation performs around 18% better, since the MSE of a standard Matérn 5/2 kernel is 2.8 for the same dataset. Hence, we provide empirical evidence that the commonly used parameterisations of *v* are far from optimal in terms of inference accuracy in various scenarios. These results transform the discussion around kernel choice into a standard speed-accuracy trade-off problem. Using half-integer *v* values is substantially faster, yet may also be far less accurate.
<table><caption>**Table 1.** GP performance on the 3d Schwefel dataset with $v$ included in the hyperparameter optimisation. $n$ is the number of datapoints, $t$ the CPU (12 core Intel X5650) time, in hours, to complete landscape exploration. The MSE is the best MSE of all identified minima, the corresponding $v$ is the $v$-value of the minimum with the best MSE.</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">#minima</th><th rowspan="2">t</th><th colspan="2">test</th><th colspan="2">train</th></tr><tr><th>MSE</th><th>v</th><th>MSE</th><th>v</th></tr></thead><tbody><tr><td>100</td><td>42</td><td>0.8</td><td>280.7</td><td>2.29</td><td>1.9 × 10<sup>-10</sup></td><td>1.01</td></tr><tr><td>200</td><td>14</td><td>9.3</td><td>125.1</td><td>5.15</td><td>1.1 × 10<sup>-08</sup></td><td>5.15</td></tr><tr><td>300</td><td>10</td><td>10.1</td><td>43.5</td><td>1.55</td><td>1.6 × 10<sup>-09</sup></td><td>3.53</td></tr><tr><td>400</td><td>10</td><td>18.5</td><td>28.3</td><td>4.83</td><td>1.9 × 10<sup>-08</sup></td><td>4.83</td></tr><tr><td>500</td><td>15</td><td>35.2</td><td>19.1</td><td>3.77</td><td>1.8 × 10<sup>-08</sup></td><td>3.77</td></tr><tr><td>600</td><td>10</td><td>67</td><td>7.7</td><td>4.74</td><td>4.1 × 10<sup>-08</sup></td><td>4.74</td></tr><tr><td>700</td><td>7</td><td>79</td><td>6.9</td><td>4.34</td><td>7.0 × 10<sup>-08</sup></td><td>4.34</td></tr><tr><td>800</td><td>8</td><td>94.2</td><td>2.8</td><td>4.67</td><td>9.1 × 10<sup>-08</sup></td><td>4.67</td></tr><tr><td>900</td><td>13</td><td>176.1</td><td>2.3</td><td>4.45</td><td>1.2 × 10<sup>-07</sup></td><td>4.45</td></tr><tr><td>1000</td><td>7</td><td>181.2</td><td>2.6</td><td>4.87</td><td>2.7 × 10<sup>-07</sup></td><td>4.87</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 3.** GP ensemble methodology for a model prediction problem. Various minima of the loss landscape are obtained via landscape exploration. At inference time, each minimum is individually employed to provide an answer. Each answer is given a weight by one of the methods described below. The final answer from the ensemble is a weighted combination of all the predictions.</figcaption></figure>

### 3.3. Physics-inspired ensemble learning

Combining the outputs of multiple predictors has become common practice in machine learning. Instead of relying on a single model, multiple ML models perform inference on the same task and a weighted aggregation of these individual inferences forms the final prediction. There are relatively few applications of GP ensembles, due both to computational cost and insufficient model performance. We discuss a novel ensembling approach that aims to tackle both issues.

The ensemble scheme considered here is based on the loss landscapes described in previous sections, and a detailed account of the ideas behind landscape-inspired ensembles is given in [15]. Usually, the multiple ML models that form an ensemble result from repeated stochastic gradient descent beginning from new random initialisations of the model. Each minimum of the loss landscape is a locally-optimal set of hyperparameters that model the training data. In contrast to this conventional ensemble generation process, the landscape-inspired methodology guarantees unique fully converged minima, mitigating the risk of including repetitions of the same minimum. A visual depiction of this scheme is given in figure 3. Furthermore, by employing methods from the energy landscapes framework, we exploit knowledge and insights from a different scientific field to construct, and weight, GP ensembles in a novel way. In various experiments, we show that these ensembles substantially outperform any single minimum of the same loss landscape by a significant margin for various synthetic and real-life datasets. Furthermore, we show that the choice of voting or weighting scheme between individual classifiers is extremely important and greatly impacts model accuracy.
<figure><img src="image_10.png" /><figcaption>**Figure 4.** Percentage improvement over single best minimum of GP ensembles for fitting the (left) 3d and (right) 4d Schwefel function for various parameterisations of *v*. Different weighting strategies are highlighted in respective colours.</figcaption></figure>

### 3.3.1. Ensemble methodology and performance

Combinations of the individual predictions from all the considered minima can be performed by simple majority vote or uniform-weight based prediction across all the considered minima. However, we have also tested various other weighting schemes, inspired by common methods in the physical sciences, and the results are shown in table 1. The three schemes considered are weighting by loss value, occupation probability, and spectral norm of the Hessian.

**Loss value:** Weighting ensemble minima by loss value makes the prediction a function of performance on the training set. We also considered weighting by inverse training MSE, but we did not observe a visible difference between these schemes. Weighting by loss value performs very well for all datasets, as shown in table 2. This performance highlights both the value of the best minimum in the loss landscape and the fact that evidently, additional information contained in other minima contributes to an improvement during inference.

**Occupation probability:** The occupation probability is a function of the log product of positive Hessian eigenvalues at a given local minimum; a detailed account of the computation is included in appendix B. This quantity is the analogue of the equilibrium distribution for a set of states in statistical mechanics, which depends on a temperature parameter, as explained in the appendix. The occupation probability includes a measure of the width or flatness of a minimum. Minimum flatness is often seen as a proxy for robustness to perturbations in the input data [50–52]. A higher occupation probability may therefore indicate a wider and more robust minimum. Thus, weighting by occupation probability is expected to give more robust minima a higher weight, leading to better test MSE. We can see significant improvements in test MSE for this weighting method in table 2.

**Spectral norm of Hessian:** The Hessian norm at a minimum (detailed in appendix B) provides a geometric interpretation of flatness, which, like the occupation probability, is expected to correlate with robustness. Hence, it is no surprise that occupation probability and Hessian norm perform similarly, with the Hessian performing slightly better on some datasets.

The performance of GP ensembles for different values of *v* for 3- and 4-dimensional Schwefel data are shown in figure 4. As above, a largely smooth progression across values of *v* is observed, more so in the less complex 3d Schwefel problem. It is noteworthy that ensembles outperform the single best minimum in almost all cases, and that landscape-guided weighting schemes almost always outperform simple majority vote combinations. This result strengthens the argument for using carefully crafted GP ensembles, particularly when designed under the consideration of loss landscape features.

Table 2 describes the best performance of each ensemble method for various datasets and shows how computing geometric features of the loss landscape can help to improve the model at inference time. The QSAR biodegradation dataset [53] is a collection of 1055 chemical molecules characterised by 41 molecular
<table><caption>**Table 2.** MSE improvement of GP ensemble as percentage compared to single best minimum from global optimisation. Improvements for each of the weighting schemes considered are shown, with best performing scheme highlighted.</caption><thead><tr><th></th><th>Schwefel 3d</th><th>Schwefel 4d</th><th>QSAR</th></tr></thead><tbody><tr><td>Unweighted</td><td>26.8</td><td>12.3</td><td>-12.6</td></tr><tr><td>NLML</td><td><strong>31.1</strong></td><td>15.5</td><td>19.1</td></tr><tr><td>Occupation probability</td><td>30.9</td><td>15.5</td><td>11.0</td></tr><tr><td>Spectral norm of Hessian</td><td>30.6</td><td><strong>16.1</strong></td><td><strong>22.9</strong></td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 5.** Correlation between MSE improvement achieved by ensemble methods and the number of minima in a loss landscape. The underlying data is aggregated over multiple regression problems, both real-world and synthetic.</figcaption></figure>

descriptor attributes with the aim of predicting the bioconcentration factors, which here serves as a proxy for toxicity. For each of the datasets, both synthetic and from the real world, weighting the minima by any landscape feature outperforms the unweighted combination by around 20%–25%.

### 3.3.2. To ensemble or not to ensemble

The decision to use an ensemble of models in a given ML problem involves a trade-off between accuracy and computational cost. Ensembles often improve inference, but require substantial additional computational resources, especially if a large number of individual predictors are queried. Estimating the value of an ensemble, relative to a single model, prior to construction, will allow better resource allocation. Hence, a means of a priori estimating whether ensembles are beneficial would be useful. Here, we present and test an empirical method to estimate the benefit of using ensembles by considering features of the loss landscape.

Figure 5 reveals a positive correlation between MSE improvement by ensembles and the number of minima in the loss landscape. When the number of minima is small, in our cases < 10, ensembles perform worse than the single best minimum. Some minima provide poor solutions that corrupt the weighted average in the ensemble generation process. The small numbers of minima in GPs compared to neural networks [12, 18] can, as discussed above, be largely attributed to the lower dimensionality of the optimisation space. In GPs, only a handful of hyperparameters are optimised when minimising the NLML, as opposed to billions in large neural networks.

The plots in figure 5 showing weighting-strategy correlation provide further evidence that ensembles require a large number of minima to average out bad predictors. All the weighting schemes have low correlation, indicating that they successfully manage to select against bad minima, while the unweighted ensemble exhibits extremely strong correlation. Therefore, GP ensembles using simple majority voting are usually insufficient for models that feature small numbers of minima in the loss landscape. This result also
strengthens the argument for landscape-based weighting schemes in GP ensembles, such as those shown in figure 4. These more elaborate weighting schemes discriminate against bad minima and, hence, are vital for accurate inference.

## 4. Conclusions

Understanding the loss landscapes of GPs has received surprisingly little attention, which is likely due to the lack of readily available methodology and high computational cost. We have demonstrated that knowledge about the loss landscape structure can improve GP model performance, provide increased confidence in the trained model, and decrease computational cost. Naturally, computing large areas of the loss landscape to sample minima comes at a cost, but this cost is relatively low compared to extensively sampling and computing solutions for various minima in common ensemble methods. However, a speed-accuracy/ explainability trade-off is inherent to the loss landscape approach for studying GPs. The same trade-off occurs when considering non-standard (half integer) parameterisations of the Matérn kernel, which may lead to substantially lower inference errors at higher computational cost. Lastly, we have provided empirical evidence that transitions between adjacent values of $v$ in Matérn kernels occur smoothly, which provides a novel insight into the structure of the solution space for GPs.

An exciting area of future work in GP loss landscapes will be improved hyperparameter sampling methods based on geometric properties of the loss landscape using methods from chemical physics. A fully Bayesian treatment of a GP would replace elimination of single point estimates for hyperparameters by sampling from a distribution of the hyperparameters. Currently, this sampling is often performed using Hamiltonian Monte Carlo (HMC) [54] methods, which are computationally expensive. Such HMC methods could potentially be replaced by a landscape approximation constructed from the Hessian around minima, and transition state information. This approximate landscape could be used as the inverse sampling distribution for a Bayesian treatment of hyperparameter sampling at a substantially lower cost compared to prevalent methods.

## Data availability statement

No new data were created or analysed in this study.

## Acknowledgments

MPN acknowledges funding from Downing College, Cambridge. LD and EOP-K would like to acknowledge this work was supported by the Hartree National Centre for Digital Innovation—a collaboration between Science and Technology Facilities Council and IBM. LD also acknowledges the financial support of the EPSRC via a knowledge transfer fellowship with Grant Number RG90413 MBAG/660 EP/R511675/1. DJW gratefully acknowledges an International Chair at the Interdisciplinary Institute for Artificial Intelligence at 3iA Cote d'Azur, supported by the French government, with reference number ANR-19-P3IA-0002, which has provided interactions that furthered the present research project.

## Appendix A. Minima for varying $v$

In figure A.1 we present the loss landscapes at all values of $v$, visualised using disconnectivity graphs. These graphs contain the same minima as in figure 1, but display the additional information about the barriers between minima that reveal the intermediate structure of the surface. We observe a smooth change in these graphs for changing $v$, with the number of minima retained and the loss values of minima smoothly changing too. There are two changes to the number of minima at $v: 2.0 \rightarrow 2.1$ and $v: 3.4 \rightarrow 3.5$, which result from fold catastrophes. Around the catastrophe the barrier between minima decreases until it vanishes, where the transition state and the higher-valued minimum merge at a non-Morese point.

We can show that differences between minima do not exist at the level of individual values of $v$. Different values of $v$ cannot be distinguished when minima from various $v$ runs are shown in a single plot, figure A.2.

This result should not be a surprise, since the Matérn function is smooth and, hence, one expects smooth transitions in NLML values for adjacent, subsequent parameterisations of the kernel. Furthermore, the right plot of figure A.2 underlines the fact that a 2-dimensional representation is insufficient to characterise the hyperparameter space which is also expected, since every hyperparameter in a GP should be relevant, in contrast to an overparameterised neural network.
<figure><img src="image_12.png" /><figcaption>**Figure A.1.** Disconnectivity graphs visualising loss function landscapes for sequential values of $v$ when fitting the 3d Schwefel function with $n$ points. Each leaf node of the graphs is a minimum on the NLML surface. The disappearance of minima from (a) $v = 2.1 \rightarrow v = 2.0$ and (b) $v = 3.5 \rightarrow v = 3.4$ corresponds to a reduction in the number of leaf nodes.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure A.2.** PCA of individual minima for all considered values of $v$ coloured by 1) $v$ value and 2) $k$-means clustering of the minima in high-dimensional space. No clusters are recognisable to distinguish between different values of $v$.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

## Appendix B. Computing occupation probability and Hessian spectral norm

The occupation probability, $p_\alpha(T)$, for minimum $\alpha$ is computed here in the harmonic approximation to the vibrational density of states, related to the vibrational entropy for a molecular system. $p_\alpha(T)$ is a function of the log product of positive Hessian Eigenvalues, given by:

$$p_{\alpha}(T) = \frac{n_{\alpha} \exp(-L_{\alpha}/k_{B}T)/\bar{\nu}_{\alpha}^{\kappa}}{\sum_{\gamma} n_{\gamma} \exp(-L_{\gamma}/k_{B}T)/\bar{\nu}_{\gamma}^{\kappa}} \quad (6)$$

where $n_\alpha$ is the number of permutational isomers, $\nu_\alpha$ the geometric mean normal mode vibrational frequency, $L_\alpha$ the loss value (potential energy in molecular sciences) for minimum $\alpha$, $T$ the temperature (fictitious in ML systems), $k_B$ is the Boltzmann constant, and $\kappa = 3N - 6$ is the number of vibrational degrees of freedom. Simplified, in ML terms, this becomes:

$$p_{\alpha}(T) = \frac{\exp(-L_{\alpha}/k_{B}T)/\bar{\nu}_{\alpha}^{\kappa}}{\sum_{\gamma} \exp(-L_{\gamma}/k_{B}T)/\bar{\nu}_{\gamma}^{\kappa}} \quad (7)$$

which depends on the loss value of the minima and $\nu_\alpha$, a metric related to the product of positive Hessian eigenvalues, and thus basin geometry, which for minimum $\alpha$ is $\ln \prod_i^\kappa \nu_{\alpha,i}^2$.

The spectral norm of a matrix $A$, shown as $||A||$, is the maximum singular value of $A$. Usually, this is the square root of the maximum eigenvalue. However, in the case of a symmetric matrix like the Hessian $H$, the spectral norm $||H||$ is simply the maximum eigenvalue.

## Appendix C. Normalising weights for GP ensembles

Normalisation of weights for GP ensembles is required to improve the quality of the ensembles. In figure A.3 we compare normalised and non-normalised weighting schemes. Normalisation substantially improves the ensembles.

In figure A.4 we make this result even more explicit by showing the absolute difference in percentage improvement between normalised and unnormalised weighting. For all values of $\nu$, the improvement by normalisation is striking.

Simply weighting the minima by (inverse) training MSE does not yield good results, as we see in figure A.5. In fact, weighting by inverse training MSE substantially decreases performance compared to the single best minimum. The explanation is that for any minima considered below, the training MSE lies
<figure><img src="image_16.png" /><figcaption>**Figure A.4.** Difference between normalised and non-normalised weighting of GP ensembles. Let the improvement of normalised ensembles for some weighting scheme be $\delta\mathcal{L}_n^w$ and the improvement of unweighted ensembles be $\delta\mathcal{L}_u^w$, then this figure is computed as $\delta\mathcal{L}_n^w - \delta\mathcal{L}_u^w$ for all $w$.</figcaption></figure>

between $10^{-8}$ – $10^{-13}$. Hence, scaling predictions by these weights means essentially setting the weight for one minimum to be unity, and all others to be 0.

Thus, weighting by training MSE achieves the opposite of the desired effect. Only a single minimum is considered, but the best minimum by training MSE is not necessarily the one that generalises the best. This point is underlined by figure A.4, which shows that normalising the weights has no influence on the, already poor, performance of weighting by testing MSE. The values are all small, yet differ by orders of magnitude, which makes weighting with a systematic scheme nearly impossible.
<figure><img src="image_17.png" /><figcaption>**Figure A.5.** Percentage improvement over single best minimum for GP ensembles for 3d and 4d Schwefel data. Different weighting strategies are highlighted in the respective colours.</figcaption></figure>

## ORCID iDs

Maximilian P Niroomand https://orcid.org/0000-0002-7189-0456

Luke Dicks https://orcid.org/0000-0002-5278-4412

David J Wales https://orcid.org/0000-0002-3555-6645

## References

[1] Williams C K and Rasmussen C E 2006 Gaussian Processes for Machine Learning vol 2 (MIT Press)

[2] Stephenson W T, Ghosh S, Nguyen T D, Yurochkin M, Deshpande S K and Broderick T 2021 Measuring the robustness of Gaussian processes to kernel choice (arXiv:2106.06510)

[3] Duvenaud D, Lloyd J, Grosse R, Tenenbaum J and Zoubin G 2013 Structure discovery in nonparametric regression through compositional kernel search *Int. Conf. on Machine Learning* (PMLR) pp 1166–74

[4] Liu H, Ong Y-S, Shen X and Cai J 2020 When Gaussian process meets big data: a review of scalable GPs *IEEE Trans. Neural Netw. Learn. Syst.* 31 4405–23

[5] Lawrence N, Seeger M and Herbrich R 2002 Fast sparse Gaussian process methods: the informative vector machine *Advances in Neural Information Processing Systems* vol 15

[6] Titisias M 2009 Variational learning of inducing variables in sparse Gaussian processes *Proc. 12th Int. Conf. on Artificial Intelligence and Statistics* (PMLR) pp 567–74

[7] Allison R, Stephenson A and Pyzer-Knapp E O 2023 Leveraging locality and robustness to achieve massively scalable Gaussian process regression *Advances in Neural Information Processing Systems* vol 36, ed A Oh, T Naumann, A Globerson, K Saenko, M Hardt and S Levine (Curran Associates, Inc.) pp 18906–31

[8] Sun R, Li D, Liang S, Ding T and Srikant R 2020 The global landscape of neural networks: an overview *IEEE Signal Process. Mag.* 37 95–108

[9] Bosman A S, Engelbrecht A and Helbig M 2020 Visualising basins of attraction for the cross-entropy and the squared error neural network loss functions *Neurocomputing* 400 113–36

[10] Li H, Xu Z, Taylor G, Studer C and Goldstein T 2018 Visualizing the loss landscape of neural nets *Advances in Neural Information Processing Systems* vol 31

[11] Wales D J 2003 *Energy Landscapes: Applications to Clusters, Biomolecules and Glasses* (Cambridge University Press)

[12] Verpoort P C, Lee A A and Wales D J 2020 Archetypal landscapes for deep neural networks *Proc. Natl Acad. Sci.* 117 21857–64

[13] Niroomand M P, Morgan J W, Cafolla C T and Wales D J 2022 On the capacity and superposition of minima in neural network loss function landscapes *Mach. Learn.: Sci. Technol.* 3 025004

[14] Dicks L and Wales D J 2022 Elucidating the solution structure of the k-means cost function using energy landscape theory *J. Chem. Phys.* 156 054109

[15] Niroomand M P, Dicks L, Pyzer-Knapp E O and Wales D J 2024 Insights into machine learning models from chemical physics: an energy landscapes approach (EL for ML) *Digit. Discovery* 3 637–48

[16] Fort S, Hu H and Lakshminarayanan B 2019 Deep ensembles: a loss landscape perspective (arXiv:1912.02757)

[17] Eustratiadis P, Gouk H, Li D and Hospedales T 2022 Attacking adversarial defences by smoothing the loss landscape (arXiv:2208.00862)

[18] Niroomand M P, Cafolla C T, Morgan J W and Wales D J 2022 Characterising the area under the curve loss function landscape *Mach. Learn.: Sci. Technol.* 3 015019

[19] Niroomand M P and Wales D J 2023 Physics-inspired interpretability of machine learning models *ICLR 2023 Workshop on Physics for Machine Learning*
[20] Baldassi C, Lauditi C, Malatesta E M, Perugini G and Zecchina R 2021 Unveiling the structure of wide flat minima in neural networks Phys. Rev. Lett. 127 278301

[21] Ruiz-Garcia M, Zhang G, Schoenholz S S and Liu A J 2021 Tilting the playing field: dynamical loss functions for machine learning Int. Conf. on Machine Learning (PMLR) pp 9157–67

[22] Choromanska A, Henaff M, Mathieu M, Arous G B and LeCun Y 2015 The loss surfaces of multilayer networks Artificial Intelligence and Statistics (PMLR) pp 192–204

[23] Wu L et al 2017 Towards understanding generalization of deep learning: perspective of loss landscapes (arXiv:1706.10239)

[24] Rudin C, Chen C, Chen Z, Huang H, Semenova L and Zhong C 2022 Interpretable machine learning: fundamental principles and 10 grand challenges Stat. Surv. 16 1–85

[25] Cheng L, Ramchandran S, Vatanen T, Lietzén N, Lahesmaa R, Vehtari A and Lähdesmäki H 2019 An additive Gaussian process regression model for interpretable non-parametric analysis of longitudinal data Nat. Commun. 10 1798

[26] Yoshikawa Y and Iwata T 2021 Gaussian process regression with interpretable sample-wise feature weights IEEE Trans. Neural Netw. Learn. Syst. 34 5789–803

[27] Martinez-Ferrer L, Piles M and Camps-Valls G 2020 Crop yield estimation and interpretability with Gaussian processes IEEE Geosci. Remote Sens. Lett. 18 2043–7

[28] Liu K, Li Y, Hu X, Lucu M and Widanage W D 2019 Gaussian process regression with automatic relevance determination kernel for calendar aging prediction of lithium-ion batteries IEEE Trans. Ind. Inform. 16 3767–77

[29] Duvenaud D K, Nickisch H and Rasmussen C 2011 Additive Gaussian processes Advances in Neural Information Processing Systems vol 24

[30] Breiman L 1996 Bagging predictors Mach. Learn. 24 123–40

[31] Lakshminarayanan B, Pritzel A and Blundell C 2017 Simple and scalable predictive uncertainty estimation using deep ensembles Advances in Neural Information Processing Systems vol 30

[32] Lu Q, Karanikolas G, Shen Y and Giannakis G B 2020 Ensemble Gaussian processes with spectral features for online interactive learning with scalability Int. Conf. on Artificial Intelligence and Statistics (PMLR) pp 1910–20

[33] Deng Z, Zhou F, Chen J, Wu G and Zhu J 2022 Deep ensemble as a Gaussian process approximate posterior (arXiv:2205.00163)

[34] Lu Q, Karanikolas G V and Giannakis G B 2022 Incremental ensemble Gaussian processes IEEE Trans. Pattern Anal. Mach. Intell. 45 1876–93

[35] scikit-learn developers 2023 sklearn.Gaussian_process.kernels.matern, scikit-learn 1.2.1 documentation (available at: http://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.Matern.html)

[36] Genton M G 2001 Classes of kernels for machine learning: a statistics perspective J. Mach. Learn. Res. 2 299–312

[37] Borovitskiy V et al 2020 Matérn Gaussian processes on riemannian manifolds Advances in Neural Information Processing Systems vol 33 pp 12426–37

[38] Prakash A K, Xu S, Rajagopal R and Noh H Y 2018 Robust building energy load forecasting using physically-based kernel models Energies 11 862

[39] Doctor Z, Farr B, Holz D E and Pürrer M 2017 Statistical gravitational waveform models: what to simulate next? Phys. Rev. D 96 123011

[40] Geoga C J, Marin O, Schanen M and Stein M L 2022 Fitting matern smoothness parameters using automatic differentiation (arXiv:2201.00090)

[41] Li Z and Scheraga H A 1987 Monte carlo-minimization approach to the multiple-minima problem in protein folding Proc. Natl Acad. Sci. 84 6611–5

[42] Wales D J and Doye J P 1997 Global optimization by basin-hopping and the lowest energy structures of lennard-jones clusters containing up to 110 atoms J. Phys. Chem. A 101 5111–6

[43] Wales D J and Scheraga H A 1999 Global optimization of clusters, crystals and biomolecules Science 285 1368–72

[44] Becker O M and Karplus M 1997 The topology of multidimensional potential energy surfaces: theory and application to peptide structure and kinetics J. Chem. Phys. 106 1495–517

[45] Wales D J, Miller M A and Walsh T R 1998 Archetypal energy landscapes Nature 394 758–60

[46] Zeeman E C 1976 Catastrophe theory Sci. Am. 234 65–83

[47] Wales D J 2001 A microscopic basis for the global appearance of energy landscapes Science 293 2067–9

[48] Leonardos S and Piliouras G 2022 Exploration-exploitation in multi-agent learning: catastrophe theory meets game theory Artif. Intell. 304 103653

[49] Janz D, Burt D and González J 2020 Bandit optimisation of functions in the matérn kernel RKHS Int. Conf. on Artificial Intelligence and Statistics (PMLR) pp 2486–95

[50] Hochreiter S and Schmidhuber J 1997 Flat minima Neural Comput. 9 1–42

[51] Hinton G E and Van Camp D 1993 Keeping the neural networks simple by minimizing the description length of the weights Proc. Sixth Annual Conf. on Computational Learning Theory pp 5–13

[52] Zhang Y, Saxe A M, Advani M S and Lee A A 2018 Energy–entropy competition and the effectiveness of stochastic gradient descent in machine learning Mol. Phys. 116 3214–23

[53] Mansouri K, Ringsted T, Ballabio D, Todeschini R and Consonni V 2013 Quantitative structure–activity relationship models for ready biodegradability of chemicals J. Chem. Inf. Model. 53 867–78

[54] Saatçi Y, Turner R D and Rasmussen C E 2010 Gaussian process change point models Int. Conf. on Machine Learning