{'fraction_non_alphanumeric': 0.04189815575791446, 'fraction_numerical': 0.013314047141589139, 'mean_word_length': 4.484944103858637, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GaussED: A Probabilistic Programming Language for Sequential Experimental Design

Sequential algorithms are popular for experimental design, enabling emulation, optimisation and inference to be efficiently performed. For most of these applications bespoke software has been developed, but the approach is general and many of the actual computations performed in such software are identical. Motivated by the diverse problems that can in principle be solved with common code, this paper presents GaussED, a simple probabilistic programming language coupled to a powerful experimental design engine, which together automate sequential experimental design for approximating a (possibly nonlinear) quantity of interest in Gaussian processes models. Using a handful of commands, GaussED can be used to: solve linear partial differential equations, perform tomographic reconstruction from integral data and implement Bayesian optimisation with gradient data.


Introduction
This paper concerns the development of a probabilistic programming language (PPL) for sequential experimental design (SED). A PPL is an attempt to streamline the process of performing computation with a statistical model (Goodman, 2013). SED is often associated with computational workflows that are complicated and cumbersome, as one is required to iterate between designing an experiment (to augment a dataset with a new datum) and performing inference for a specified quantity of interest (based on the augmented dataset). Thus SED is well-placed to benefit from the development of a high-level PPL. The research challenge here is to identify a class of statistical models that are sufficiently general to include important applications of SED, while being sufficiently narrow to permit both inference and SED to be efficiently and automatically performed. This paper aims to address two important open problems in PPL for SED: P1 automate SED for Gaussian process (GP) models with general nonlinear quantities of interest, in the setting of continuous linear functional data (e.g. function values, gradients, integrals); P2 circumvent the requirement for the user to specify an acquisition function for SED, in the spirit of AutoML (Hutter et al., 2014).
In limiting attention to the relatively narrow class of GP models in P1, we aim to develop more powerful algorithms than would have been possible in a more general-purpose PPL. The setting of P1 includes SED for the important tasks of emulating computer models (Kennedy and O'Hagan, 2001), performing Bayesian optimisation (Shahriari et al., 2015), and running probabilistic numerical methods (Hennig et al., 2015). Bespoke PPLs have been developed for these individual tasks, but many of the actual computations performed in such software are identical. Indeed, in Section 3 we demonstrate how a single PPL can: solve partial differential equations using a probabilistic numerical method, perform tomographic reconstruction from integral data, implement Bayesian optimisation with gradient data, and emulate a complex computer model. Such a PPL enables advances in computational methodology to be immediately brought to bear on diverse application areas where SED is performed. Existing PPLs for SED require the user to specify an acquisition function, which is used to select the next experiment and serves to control the exploration-exploitation trade-off. Unfortunately, the process of determining an effective acquisition function requires domain expertise and, while several choices have been documented in the literature (see e.g. Wilson et al., 2018, for acquisition functions in Bayesian optimisation), many problems that fall into the setting of P1 have not received such detailed treatment. In removing the technical burden of prescribing the acquisition function in P2, we may sacrifice a degree of performance relative to dedicated software for tasks such as Bayesian optimisation, for which bespoke acquisition functions have been developed. However, empirical results in this paper suggest that the loss of performance may be modest, and in turn we are able to considerably expand the applicability of the PPL.

Our Contribution
In this paper we present GaussED, a simple PPL coupled to a powerful experimental design engine for performing SED in the nonparametric GP context. GaussED achieves the aims P1 and P2, just outlined. To achieve P1, and to ensure that GaussED can handle data arising from general continuous linear functionals, we present a rigorous probabilistic treatment of conditioning for GPs. This enables us to, for example, prevent attempts to condition on a derivative that does not exist under the GP model. To achieve P2 and circumvent the user-specification of an acquisition function, we adopt a classical but surprisingly overlooked decisiontheoretic approach to SED, which requires only the quantity of interest and a loss function to be specified. The loss function quantifies the loss incurred when the true quantity of interest is approximated, a notion that is meaningful in the applied context and comparatively straightforward to elicit. The computational backend for GaussED comprises a spectral GP, a reparametrisation trick, and stochastic optimisation over the experimental design set.

Related Work
Several general-purpose PPLs have been developed for Bayesian parameter inference in parametric models (e.g. Wood et al., 2014;Carpenter et al., 2017;Bingham et al., 2019), often based on Markov chain Monte Carlo or variational approximations in the backend. Specialised PPLs have been developed for inferring parameters that minimise a predictive loss (e.g. using neural networks; Paszke et al., 2019), often based on automatic differentiation and stochastic gradient descent. For inference in nonparametric models, specialised PPLs have been developed for GP models (e.g. Rasmussen and Nickisch, 2010;Matthews et al., 2017), including for numerical applications (ProbNum, 2021).
The combination of PPL and SED for general parametric models has received attention in Rainforth (2017, Chapter 11) and Ouyang et al. (2016); Kandasamy et al. (2018), who provided a high-level syntax for Bayesian SED. Several application-specific PPL have been also been developed for SED in parametric models (e.g. Liepe et al., 2013). The focus of much of the research involving parametric models centres around the computational challenge of conditioning random variables on observed data, a problem that is often difficult (Olmedo et al., 2018).
SED for nonparametric models has received considerable attention in the context of Bayesian optimisation; see the review of Shahriari et al. (2015). However, existing PPLs are specialised to this single task. More closely related to the present paper, Paleyes et al. (2019) developed a PPL called Emukit, in which computer model emulation, Bayesian optimisation, and a number of probabilistic numerical methods are automated. However, Emukit focuses on function-value data as opposed to general continuous linear functionals (c.f. P1) and requires the user to specify a suitable acquisition function (c.f. P2).
Outline: The remainder of the paper is structured as follows: Section 2 presents a detailed technical description of GaussED. Section 3 described the syntax of GaussED and presents diverse applications of SED, for which bespoke code had previously been developed but whose automation is essentially trivial using GaussED. The potential and limitations of GaussED are summarised in Section 4.

Methodology
This section presents the statistical and computational methodology used in GaussED. First, in Section 2.1, the notation and mathematical set-up are introduced. The elements of SED are outlined in Section 2.2 and a classical, but surprisingly overlooked, approach to SED is presented in Section 2.3. This decision-theoretic approach circumvents the requirement to specify an acquisition function and, moreover, enables state-of-the-art stochastic optimisation to be employed in SED, as explained in Sections 2.4 and 2.5. The hyperparameters of the GP model are estimated online during SED, as explained in Section 2.6.

Notation and Set-Up
Let F be a normed vector space of real-valued functions on some domain X ⊆ R d . The problems that we consider involve a latent function f ∈ F, associated with a high computational cost, and the task is to approximate a (possibly nonlinear) quantity of interest q(f) using SED. The experiments are represented 1 as continuous linear functionals δ : F → R and may, for example, include pointwise evaluation δ(f) = f(x) of the latent function f at a specified location x ∈ X , pointwise evaluation of a gradient, or evaluation of an integral, such as a Fourier transform. A limited computational budget motivates the careful selection of informative experiments δ 1 , . . . , δ n . SED is often preferred 2 over a priori experimental design, since it allows data δ 1 (f), . . . , δ n−1 (f), which have already been observed, to inform the design of the next functional δ n .
Bayesian statistics provides a general framework in which SED can be performed. To this end, let (Ω, S, P) be a probability space and consider a random variable f : Ω → F. This serves as a statistical model for the latent f, and encodes a priori knowledge, such as the smoothness of f. To notate the distribution of f , we first define the pre-image of a set B ⊆ F as f −1 (B) := {ω ∈ Ω : f (ω) ∈ B} and we let f # P denote the pushforward of P through f ; i.e. the probability distribution on F that assigns, to each Borel set B ⊆ F, the mass f # P(B) := P(f −1 (B)). The distribution of f will be denoted P f := f # P in the sequel. Our presentation allows for general priors for f until Section 2.5, at which point we will assume f is a GP. Throughout we adopt the convention that f refers to the latent function of interest, f is a random variable model for f, and f is a generic element of the set F.

Sequential Experimental Design
SED iterates between designing an experiment δ n , to augment a dataset with a new datum δ n (f), and performing inference for a specified quantity of interest, based on the augmented dataset δ n (f) := (δ 1 (f), . . . , δ n (f)) . Let D ⊆ F indicate the design set, where F is the topological dual space of F, containing the continuous linear functionals on F. The design set D will depend on the problem at hand, and contains only the experiments that can actually be performed. At iteration n, SED selects an experiment δ n from the design set in order that an acquisition function is maximised 3 : The role of the acquisition function A is to control the exploration-exploitation trade-off, but the computational convenience of computing (1) is also important. Much research has been dedicated to exploring choices for A, and the statistical and computational properties of the associated sequence (δ n ) ∞ n=1 . Specific applications, where interest is not necessarily in f but rather a derived quantity of interest q(f), have developed bespoke acquisition functions that balance computational cost with accurate approximation of the quantity of interest, in particular in Bayesian optimisation (see Table 1 in Wilson et al., 2018). This presents a major problem (P2) for the development of a general purpose PPL for SED, since in general we cannot expect a user to specify a suitable acquisition function for the problem at hand.
As a first step toward solving P2, we consider a Bayesian approach to the design of an acquistion function. To this end, let P f (·|δ n (f)) denote the conditional distribution (or posterior ) of f obtained by setting the values δ n (f ) equal to the observed data δ n (f). From a mathematical perspective, the proper construction of a conditional distribution for an infinite-dimensional random variable f is non-trivial; we suppress further discussion in the main text but refer the reader to Appendix A for full mathematical detail. A Bayesian approach to the design of an acquisition function is then to let U : R n−1 × R → R be a utility function, to be specified, and to seek an experiment for which the current expected utility is maximised. The utility U (δ n−1 (f), δ(f)) represents the value to the user of observing the datum δ(f). Thus the design of an acquisition function can be reduced to the design of a utility function. A popular default choice for U is the information gain (Lindley, 1956) which quantifies the extent to which observation of the datum δ(f) changes a posteriori belief; here KL denotes the Kullback-Leibler divergence. For related approaches and discussion see the recent survey in Kleinegesse and Gutmann (2021). However, in the setting where data are exactly observed, the two distributions in (3) will be mutually singular and the Kullback-Leibler divergence will not exist. This renders information-based acquisition functions such as (3) unsuitable for our PPL. Instead, we propose to revisit a classical but often overlooked idea from experimental design, next.

A Decision-Theoretic Approach
A general approach to construction of a utility U is provided by Bayesian decision theory in the parameter inference context 4 . Let L : F × F → R denote the loss L(f, g) when estimating the function (or parameter ) f by g. Then we can take U to be the negative Bayes' risk − min g∈F L(g, g ) dP f (g |δ n−1 (f), δ(f)), which corresponds to the negative expected loss when the Bayes decision rule g is used. Compared to an acquisition function or a utility function, it can be more straightforward to specify a suitable loss function L, since no consideration of the design set is required. Although appealing in terms of its generality, the presence of the optimisation over g has historically rendered this utility unappealing from a computational viewpoint, and motivated more convenient choices, such as (3), that have since become canonical (see the survey in Chaloner and Verdinelli, 1995). However, we argue that the presumed intractability of loss-based utilities might need to be revisited in light of modern and powerful stochastic optimisation techniques. Indeed, for loss functions of the form L(f, g) = q(f) − q(g) 2 , indicating that one has a quantity of interest q(f) taking values in a normed space 5 , under mild conditions (4) is equal to The required regularity conditions and a formal proof are contained in Appendix B. At first glance it is unclear why this observation is helpful, since we have replaced an optimisation problem with an integration problem, and integration is typically more difficult than optimisation. However, this formulation turns the experimental design problem to find δ n into a double expectation and, if the design set D has enough structure for calculus, then gradient-based stochastic optimisation can be applied. The restriction to squared error loss is not as limited as it may first appear, since one has the freedom to specify the quantity of interest q(f) in such a way that application of squared error loss to q(f) captures salient aspects of the task at hand. Concrete examples of this are provided in Section 3.2.
( 6) This acquisition function does not permit a closed form in general. Several numerical methods have been proposed for maximisation of acquisition functions in the literature, including Bayesian optimisation (Overstall and Woods, 2017;Kleinegesse and Gutmann, 2019), non-gradient based Monte-Carlo methods, and approximation strategies. Similar to the approach 6 of Wilson et al. (2018), here we consider the use of stochastic optimisation techniques (Robbins and Monro, 1951) for selecting an experiment δ for which (6) is approximately maximised. For an overview of stochastic optimisation, see Kushner and Yin (2003); Ruder (2016). First we perform a reparametrisation trick (Kingma and Welling, 2014), expressing using a deterministic transformation η of a random variable ω that is δ-independent. Section 2.5, below, details how we applied the reparametrisation trick to a GP model. Now, suppose further that the elements of the design set can be parametrised as D = {δ z } z∈R m ⊆ F. Assuming sufficiently regularity for the following calculus to be well-defined, an unbiased estimator of the gradient of the acquisition function is , where the g i are independent random variables with distribution P f (·|δ n−1 (f)) and the ω ij are independent random variables with distribution P. This is an instance of nested Monte Carlo. The optimal balance between N and M for a fixed computational budget is discussed in Rainforth et al. (2018); for a continuously differentiable gradient, an optimal choice 7 is N ∝ M 2 .
GaussED exploits state-of-the-art spectral GPs to perform the reparametrisation trick, as presented next.

Spectral Approximation of GPs
Up to this point our discussion applied to general statistical models P f for the latent function f. In the remainder GPs will be used, since they facilitate closed form conditional distributions, as appearing in (6). The purpose of this section is twofold; to briefly introduce GPs and to describe how the reparametrisation trick can be performed.
A random variable f taking values in a normed vector space F is Gaussian if, for every continuous linear functional δ : F → R, the random variable δ(f ) is a Gaussian on R; see Definition 2.41 in Sullivan (2015). It follows that the statistical properties of a GP are characterised by its mean function µ(x) := E[f (x)], x ∈ X , and covariance function k(x, y) := C[f (x), f (y)], x, y ∈ X , and we write f ∼ GP(µ, k). GPs admit conjugate inference, meaning that for a continuous linear functional δ ∈ D, the conditional distributions P f (·|δ(f)) are also Gaussian, with mean and covariance functions that can be computed in closed form; see Appendix C.1.
For the reparametrisation trick, we aim to write a GP as a deterministic transformation f = η(ω) of a random variable ω, such that the distribution of ω does not depend on µ or k. However, being a nonparametric statistical model, an infinite-dimensional ω will in general be required. This motivates the use of an accurate finite-dimensional approximation of a GP at the outset, i.e. for the prior P f . A truncated Karhunen-Loeve expansion (see e.g. Theorem 11.4 in Sullivan, 2015) in principle provides such a transformation, however this requires computation of the eigenfunctions of k, and linear functionals thereof, which will in general be difficult. The solution adopted in GaussED is to use the finite-rank approximation to isotropic GPs introduced in Solin and Särkkä (2019) are independent, s is the spectral density of k, and (φ i , λ i ) are the pairs of eigenfunctions and eigenvalues of the Laplacian ∆ over the domain X ; see Appendix D for detail. The approximation converges as m → ∞, with small values of m often sufficient for accurate approximation; see Riutort-Mayol et al. (2020). GaussED puts the user in control of m, since m is the principal determinant of computational complexity in the experimental design engine, aside from the computations involving the latent function f itself.

Hyperparameter Estimation
To this point we assumed that a GP model can be specified at the outset. In reality one is usually prepared only to posit a parametric class of GPs whose parameters (called hyperparameters) are jointly estimated. In GaussED the hyperparmaters of the GP are estimated at each iteration n ≥ n 0 of SED, using the available dataset δ n (f), after an initial number n 0 ∈ N of data have been observed. Maximum likelihood estimation is employed, facilitated using automatic differentiation and Adam (Kingma and Ba, 2015). The role of n 0 is to guard against over-confident inferences, since maximum likelihood tends to overfit when the dataset is small; see e.g. Chapter 5 of Rasmussen and Williams (2006). In GaussED, the default value is taken as n 0 = 10.
This completes our description of GaussED. Our attention turns, next, to demonstrating and assessing its capabilities.

Demonstration
The aims of this section are to validate GaussED and to highlight the diverse and non-trivial applications that can be tackled. GaussED is based on Python and utilises the automatic differentiation capabilities of Pytorch (Paszke et al., 2019). Source code and documentation for GaussED can be downloaded from https://github.com/MatthewAlexanderFisher/GaussED.
Full details for each of the following examples are provided in Appendix F. An investigation into the sensitivity of the computational methodology to initial conditions, the choice of stochastic optimisation method, and the number of basis functions m, can be found in Appendix G.

Probabilistic Solution of PDEs
Our first example concerns the probabilistic numerical solution of Poisson's equation with Dirichlet boundary conditions; the intention is to validate our methodology on a problem that is well-understood. SED for such problems was investigated with bespoke code in Cockayne et al. (2016). The PDE we consider is defined on X = [−1, 1] 2 and takes the form Our quantity of interest is the solution f and the black-box source g is assumed to be associated with a computational cost, so that numerical uncertainty quantification is required. For this demonstration we simply took as a test bed. The latent f was modelled as a GP f with mean zero and Matérn covariance with smoothness parameter ν = 3 + 1 2 , ensuring the corresponding GP samples are almost surely contained in C 3 (X ), implying the evaluations of the Laplacian of f are continuous linear functionals (see Appendix C.2). The design set D, parameterised by x ∈ X , consists of functionals of the form δ(f) = ∆f(x). It is known that an optimal experimental design in this case is space filling (Wendland, 2004;Novak and Woźniakowski, 2010), as quantified by the fill distance and this fact will be used to validate GaussED. The syntax of GaussED is demonstrated in Figure 1, and consists of specifying a covariance function (k), a quantity of interest (qoi), an observation model (obs), here the Laplacian (Laplace), a loss function (loss), a design (d) initialised with an initial_design, and an acquisition function (acq). BayesRisk is the default acquisition function from (6), but GaussED retains the capability for alternative acquisition functions in the event that they can be user-specified. The experiment object (experiment) then collates these objects together to perform n = 150 iterations of SED, optimising hyperparameters as specified in Section 2.6. Results are shown in Figure 2 and required only the 8 lines of code shown in Figure 1. The number of basis functions used was m = 30 2 , we computed n 0 = 10 iterations of SED before beginning hyperparameter optimisation and a total of 9 CPU hours were invested to ensure that all n = 150 instances of stochastic optimisation converged. The fill distance is lower-bounded by Θ(n −1/2 ), and Figure 2c demonstrates that this optimal rate is empirically achieved by GaussED. This validates our approach to SED.

Tomographic Reconstruction
Our next example is tomographic reconstruction from x-ray data (Mersereau and Oppenheim, 1974). The aim is to reconstruct a latent function f : , is a parameterisation of a line with endpoints r(a), r(b) ∈ ∂X . SED for this problem was recently addressed, using bespoke code, in Burger et al. (2021) and Helin et al. (2021). Following Burger et al. (2021), an experiment consists of a set of 9 parallel line integrals across X , with lines a perpendicular distance of 0.03 apart. As a toy example, we consider tomographic reconstruction of an indicator function For our statistical model f we used a stationary GP with Matérn covariance and smoothness parameter ν = 2 + 1 2 , and the non-linear quantity of interest was q(f) = exp(3f) which, when combined with squared error loss, serves to prioritise the reconstruction of the ball in SED. See Appendix F.2 for full detail.  Results are shown in Figure 3 and only 32 lines of code were required. In this experiment, we used m = 28 2 basis functions and began optimising hyperparameters at SED iteration n = 1. In total, 2.5 CPU hours were required. SED using GaussED provides improved reconstruction compared to a random design (right panel). As an additional comparison, we also performed SED with the linear quantity of interest q(f) = f and a space-filling design was obtained. Exploratory investigation of this kind is straight-forward in GaussED.

Gradient-Based Bayesian Optimisation
Our next example uses Bayesian optimisation to perform parameter inference via maximum likelihood, and for this we consider the Lotka-Volterra model dp where p(t), q(t) > 0 are the predator and prey populations, respectively, at time t and α, β, γ and δ are free parameters to be inferred. To facilitate visualisation of experimental designs we consider inferring only α and β, which we collect in a single parameter vector x = (α, β). For this demonstration we restrict attention to X = [0.45, 0.9] × [0.09, 0.5], to avoid failure of the numerical integrator applied to (8). The remaining parameters, γ and δ, are then taken as fixed. Our latent function f is the log-likelihood, denoted f = log L, arising from a particular dataset of noise-corrupted observations described in Appendix F.3. Our quantity of interest is the maximum likelihood estimator q(f) = max x∈X f(x). The design set D contains pointwise evaluation functionals δ 1 x (f) = log L(x) and gradient evaluation functionals δ 2,i x (f) = ∇ xi log L(x), and at each iteration of SED we evaluate (δ 1 for some x ∈ X , mimicking the information provided when (8) is solved using an adjoint method. Through a suitable sequence of evaluation functionals, SED aims to approximate the maximum likelihood estimator.
Results are shown in Figure 4 and only 17 lines of code were required. In this experiment we used m = 35 2 basis functions, we computed n 0 = 10 iterations of SED before beginning hyperparameter optimisation and 1.5 CPU hours were required. For reference, results based on gradient ascent and L-BFGS (Nocedal, 1980) are also displayed. All algorithms were initialised at the midpoint of the domain X and run for n = 30 iterations. Bayesian optimisation with gradient data outperformed the first order optimisation methods in this example, where attention is focused on performance after a small number of likelihood evaluations, to mimic more challenging applications in which the likelihood is associated with a more substantial computational cost.

Discussion
This paper introduced GaussED, a simple PPL coupled to a powerful engine for SED. Through four experiments we illustrated the diverse applications that can be automatically solved using GaussED. However, automation of SED comes at a cost: Firstly, GaussED is restricted both to continuous linear functional data and to GPs, limiting the potential for more flexible statistical models to be employed. Alternative PPLs, such as Emukit, offer more modelling flexibility but require acquistion functions to be manually specified. Secondly, in automating the specification of an acquisition function in GaussED, there may be a loss in performance terms compared to bespoke solutions for specific tasks. Our experiments involving Bayesian optimisation in Section 3.3 were encouraging, however, and suggested that such performance gaps, if they do exist, may be acceptably small. One role for GaussED in these settings is to provide an off-the-shelf benchmark for SED, against which more sophisticated methods can be compared.

Supplement
The supplement is structured as follows: • Appendix A contains the mathematical preliminaries for the subsequent sections Appendix B and Appendix C.
• Appendix B presents the conditions for the equivalence of (4) and (5) as advertised in Section 2.3.
• Appendix C presents the formal background of conditioning on continuous linear data for Gaussian processes and presents properties of the Matérn covariance function.
• Appendix D presents a derivation of the Gaussian process model that forms the foundation of GaussED.
• Appendix E discusses computational aspects of GaussED. In particular, we discuss linear algebra solvers, different approaches to sampling from the posterior, and we present a complete description of how GaussED attempts to optimise the acquisition function in SED.
• Appendix F contains full details of the experiments presented in Section 3. Appendix F.1 details the partial differential equation (PDE) experiment presented in Section 3.1. Appendix F.2 details the tomographic reconstruction experiment presented in Section 3.2. Appendix F.3 details the Lotka-Volterra experiment presented in Section 3.3.
• Appendix G presents further empirical evaluation of GaussED. Appendix G.1 presents an empirical investigation of the stochastic optimisation methods used in SED. Appendix G.2 presents an empirical investigation on how the number of basis functions affects the quality of inference.

A Mathematical Preliminaries
In this section we present the mathematics required to ensure that the conditioning of stochastic processes in the main text is well-defined (Appendix A.1), as well as recalling the concept of a Fréchet derivative (Appendix A.2).

A.1 Conditioning as Disintegration
In finite dimensions, conditioning of random variables can be performed using the density formulation of Bayes' theorem. However, typical stochastic processes will be infinite-dimensional, meaning that (Lebesgue) densities do not exist in general. This necessitates a level of mathematical abstraction to ensure that conditional probabilities are well-defined. The appropriate notion, for this work, is that of disintegration, defined next. Let (F, S F ) and (Y, S Y ) be measurable spaces and let δ be a measurable function from F to Y. Recall that δ −1 (S) = {f ∈ F : δ(f ) ∈ S} denotes the pre-image of S ∈ S Y . Let P be a probability measure on (F, S F ) and recall that δ # P denotes the pushforward measure (δ # P)(S) := P(δ −1 (S)) on Y.
Definition 1. The collection {P(·|y)} y∈Y is called a δ-disintegration of P if 1. P(δ −1 (y)|y) = 1 for δ # P almost all y ∈ Y and, for each measurable function g : A disintegration is a particular instance of a regular conditional distribution which also satisfies property (1) in Definition 1; see Chang and Pollard (1997). A basic theorem on the existence and δ # P almost everywhere uniqueness of disintegrations is given in Parthasarathy (2005, p147). Two disintegrations will be identified if they coincide δ # P almost everywhere, and we will therefore refer to the δ-disintegration of P. The concept of disintegration makes precise what it means to "condition GPs on data", as discussed in Appendix C.1.

A.2 Fréchet Derivatives
Recall that F was defined as a normed vector space, meaning that the notion of a Fréchet derivative can be exploited. A function q : F → R d is called Fréchet differentiable at f ∈ F if there exists a bounded linear operator A : F → R d such that If such an operator exists it can be shown to be unique, called the Fréchet derivative of q at f, and denoted Dq(f) = A. To emphasise that the Fréchet derivative is an operator, we occasionally write Dq(f)(·) in the sequel. A Fréchet derivative Dq(f ) is said to have full rank if Dq(f)(g) = 0 implies g = 0.
The chain rule for Fréchet derivatives takes the form As a concrete example, that we use later, consider a(f) = q(f) to be the quantity of interest and b(x) = x − q(g) 2 for all x ∈ R d and some fixed g ∈ F. Then Db(x)(·) = 2 x − q(g), · is a linear operator from R d to R and we have which is a linear operator from F to R. Further background on Fréchet derivatives can be found in Berger (1977, Section 2.1C). An important technical result on Fréchet derivatives, that we will use in the sequel, is when the interchange of a Fréchet derivative and an integral can be permitted: Proposition 1. Let F be complete (i.e. a Banach space) and (Ω, S, P) be a probability space. Let : F ×Ω → R satisfy the following: is Fréchet differentiable, with derivative Proof. A special case of Kammar (2016).

B Regularity Conditions for the Decision Theoretic Formulation
The aim in this section is to establish sufficient conditions for the equivalence of (4) and (5) as advertised in Section 2.3. To achieve this, we will use the notion of a Fréchet derivative from in Appendix A.2. Our sufficient conditions are presented in Appendix B.1. A short discussion of the strength of these conditions is contained in Appendix B.2.

B.1 From Optimisation to Expectation
Firstly, we rigorously establish an infinite-dimensional analogue of the classical result that the posterior mean is a Bayes act for squared error loss: Assume that F is complete (i.e. a Banach space) and that: (A2) the Fréchet derivative Dq(f) has full rank at all f ∈ F; Then any solution to satisfies q(f) = q(g) dP f (g|δ n (f)).
Proof. From an application of Proposition 1 with (f, ω) = L(f, g(ω)), where g : Ω → F is a random variable with distribution P f (·|δ n (f)), we deduce that our assumptions on L and q (A1) are sufficient for the Fréchet derivative of r to exist. Thus, a minimiser f of (10) satisfies Dr(f) = 0. To evaluate Dr we exploit the integrability assumption (A3) on q to differentiate under the integral, which is also justified from Proposition 1: Next we apply the chain rule for Fréchet derivatives in the form of (9), yielding , Dq(f)(·) .
Proof. Let f ∈ F solve (10). Then consider the algebraic identity Using this identity, the loss function can be expressed as From linearity of the inner product we have that where we have used the integrability assumption (A3) on q to bring the integrals into the inner product, and we have used Proposition 2 to conclude that each argument is equal to 0. Finally, from the fact that g and g are identically distributed, we have which completes the argument.

B.2 Verifying the Assumptions
The main assumption in Proposition 2 is (A2); the requirement that Dq(f) has full rank for all f ∈ F. As we explain below through a worked example, (A2) is non-trivial but may often be satisfied with only minor modification to the SED task in hand. As a worked example, suppose F is a Hilbert space containing smooth, real-valued functions defined on a compact set X ⊂ R d . Suppose that we are interested in the quantity of interest Then (A2) is not satisfied in general, because q(f + g) = q(f) for all g in the linear subspace G = {g ∈ F : X g(x)dx = 0}. It follows that Dq(f)(g) = 0 for all g ∈ G, so that Dq(f) does not have full rank whenever G is non-trivial. However, (A2) is satisfied if we restrict attention to the normed vector space F c spanned by the elements of F \ G, since then Dq(f)(g) = X g(x)dx and thus Dq(f)(g) = 0 with g ∈ F c implies g = 0. This illustrates that, with a small amount of technical care, the assumptions of Proposition 2 can often be satisfied.

C Properties of Gaussian Processes
In this section, we present the formal background of conditioning on continuous linear data for Gaussian processes and detail properties of the Matérn covariance function.

C.1 Disintegration of Gaussian Measures
Let X be a compact subset of R d for some d ∈ N and let C r (X ) denote the vector space of r-times continuously differentiable real-valued functions on X equipped with the norm where the maximum ranges over multi-indices α ∈ N d 0 with |α| = α 1 +· · ·+α d ≤ r and f (α) (x) : . In what follows we consider disintegration in the case where F = C r (X ), equipped with the Borel σ-algebra, and Y = R. For an operator δ and a bivariate function k(·, ·), denote δk(·, ·) to be the action of δ on the first argument of k, and denoteδk(·, ·) to be the action of δ on the second argument of k.
Lemma 1. Let P be a Gaussian measure on C r (X ) with mean function m : X → R and covariance function k : X × X → R. Let δ : C r (X ) → R be a continuous linear functional. For each y ∈ R, define P(·|y) to be a Gaussian measure with mean and covariance function Then {P(·|y)} y∈R is a δ-disintegration of P.
The fact that the elements P(·|y) of the disintegration are again Gaussian enables the repeated application of Lemma 1, for example to condition on n ≥ 1 continuous linear functionals δ n = (δ 1 , . . . , δ n ) , as exploited in the main text. Constructed in this way, it can be verified that the elements P(·|y n ) of the resulting disintegration, with y n ∈ Y n , are invariant to the order in which the disintegrations are performed.

C.2 Sufficient Conditions for Disintegration of Matérn Processes
To exploit Lemma 1 in practice it is sufficient to verify that samples from the Gaussian measure are almost surely contained in C r (X ). Such analysis is technical but specific results are available for derivative data in the context of the tensor product Matérn covariance model that we primarily use in this work. Indeed, let P be a Gaussian measure with mean function m and covariance function k, such that m ∈ C r (X ) and where K ν denotes the modified Bessel function of the second kind and ν i := r + 1 2 . Then Theorem 2 of Wang et al. (2021a) establishes that samples are almost surely contained in C r (X ). Moreover, maps of the form δ(f ) = f (α) (x), |α| ≤ r are continuous linear functionals from C r (X ) to R, since |δ(f )| ≤ f C r (X ) for all f ∈ C r (X ). Thus, in this case Lemma 1 can be used to condition the tensor product Matérn process in (12) on the derivative data f (α) (x), safe in the knowledge that the conditional process will be well-defined.

D Spectral Approximation
This section presents an informal derivation of the spectral GP approximation of Solin and Särkkä (2019). The following utilises properties of the Fourier transform, which are first briefly recalled.

D.1 Properties of the Fourier Transform
In the following we use F to denote the Fourier transform operator and use the notationf := F (f ) to denote the Fourier transform of f . In the following we use the convention of using the angular frequency domain. Therefore, for square-integrable f : R d → R, we have Recall that, when an operator T satisfies F (T f )(ω) = m(ω)f (ω), the operator T is called a multiplier operator and the corresponding m is called the multiplier of T . As a trivial example, the identity operator T f = f is a multiplication operator, with associated multiplier 1. A more elaborate example, that is used in the subsequent section, is the Laplace operator ∆ := ∂ 2 , acting on sufficiently smooth functions f , is also a multiplier operator with multiplier (− ω 2 ) n . This can be seen by induction on the previous formula, By the convolution theorem, every multiplier operator T with multiplier m T , has an associated convolution kernel k T := F −1 (m T ) that satisfies the following where denotes convolution. Thus a multiplier operator is, in this sense, equivalent to a convolution operation. We now state two important results that define the intimate connection between covariance functions and the Fourier transform. The first result is known as Bochner's theorem (Rudin, 1990).
Theorem 1 (Bochner's theorem). A stationary covariance function, i.e. a covariance function of the form k(x, y) = k(x − y), k : R d → R, can be written as the inverse Fourier transform of a finite positive measure µ such that k(0) = µ(R d ). That is The measure µ is called the spectral measure of k and the density of µ, if it exists, is called the spectral density s(ω) of k. In the case where the spectral density s(ω) of a stationary covariance function k exists, k and s exist as Fourier duals. This result is known as the Wiener-Khintchine theorem (Khintchine, 1934).
Theorem 2 (Wiener- Khintchine theorem). Suppose that the spectral density s : R d → R of a stationary covariance function k : R d → R exists, then In the proceeding section the Wiener-Khintchine theorem and the equivalence between a multiplier operator and an associated convolution operation are both used to establish a correspondence between the covariance operator of a stationary kernel k and its spectral density s. This is the foundation upon which the spectral GP approximation of Solin and Särkkä (2019) is established.

D.2 Spectral Gaussian Processes
For every covariance function k, there exists an associated Hilbert-Schmidt integral operator, termed the covariance operator, When k is stationary, the resulting covariance operator takes the form of a convolution By the convolution theorem, we can then write the operator in the form F (Kf ) =kf and so K is a multiplier operator with multiplierk. By Theorem 2, the multiplier of K is the spectral density s =k of k.
Assuming now that the covariance function is isotropic and so satisfies k(x, y) = k( x − y ), the corresponding spectral density s of k can be written as a function of ω only and so s(ω) = S( ω ), for an appropriate function S. As a further manipulation, we can write s as a function of ω 2 only, s(ω) = ψ( ω 2 ). Assuming that ψ possesses a Taylor expansion, we can write with each µ i ∈ R. Inspired by the multiplier − ω 2 of the Laplacian in (13) and by utilising the above Taylor expansion, we can write the Fourier transform of the covariance operator of an isotropic kernel in the form By continuity of F , taking the inverse Fourier transform of the above yields a polynomial expansion form of the covariance operator The remaining step is to approximate the negative Laplacian operator. To achieve this, we write the convolution kernel k −∆ of the negative Laplacian as a Mercer expansion. To this end, we consider the following eigenvalue problem of the Laplacian over a compact domain X ⊆ R d , with boundary ∂X , with Dirichlet boundary conditions Over a suitable domain contained within L 2 (X ), the negative Laplacian is a positive definite Hermitian operator and so we can provide a Mercer expansion of the convolution kernel k −∆ of the negative Laplacian, utilising the eigenfunctions φ i . Similarly, we can provide a Mercer expansion of the convolution kernel of (−∆) n , noting that each φ i is again an eigenfunction, but now with corresponding eigenvalue λ n i . This can be seen by iteratively applying −∆ to the eigenvalue problem (15). Therefore, we have Plugging the preceding formula into equation (14) yields the following: Comparing the above form of Kf (x) to its original definition Kf (x) = k(x − y)f (y) dy implies that we can approximate k as follows where, in the final step, we utilised our Taylor expansion of the spectral density s of k and set ω 2 = λ j for each j ∈ N. Refer to the original work Solin and Särkkä (2019) for convergence analyses of the given approximation. Therefore, the resulting Gaussian model assumes the following truncated basis expansion where c i ∼ N (0, s( √ λ i )) and the φ i and λ i are the corresponding eigenfunctions and eigenvalues of the Laplacian over a compact domain X with Dirichlet boundary conditions φ i (x) = 0 on ∂X .
When the domain is the unit hypercube, X = [0, 1] d , the resulting eigenfunctions and eigenvalues can be explicitly computed as where j = (j 1 , . . . , j d ) ∈ Z d m . Taking m sinusoidal functions in each dimension yields m d eigenfunctions in total. For computational purposes, in GaussED the domain X of the Gaussian model is taken as a d-dimensional The required eigenfunctions can be obtained by a simple rescaling of the previous formula.

E Computational Details of GaussED
In this section we provide details of certain aspects of the computational approaches of GaussED. In Appendix E.1, we derive the relevant conditional distributions of the spectral Gaussian process model detailed in Section 2.5, under general linear information. In

E.1 Conditioning
In this section, we both derive and discuss GaussED's approach to conditioning and sampling from the posterior. For completeness, we present the derivation of the conditional distributions of the Gaussian process model detailed in Appendix D. For the sake of generality we consider a general truncated basis model, which takes the form of where the c i are pairwise independent Gaussian variables and the φ i form our basis functions. Suppose that we have a vector of n continuous linear functionals δ n = (δ 1 , . . . , δ n ) ∈ D n , such that each δ i belong to the design set D (see Section 2.2). We form the conditional distribution f | δ n (f) as follows, letting c = (c 1 , . . . , c m ) , we have c δ n (f ) where K cc = C(c, c) ∈ R m×m , K cδ = C(c, δ n ) ∈ R m×n , K δc = K cδ and K δδ = C(δ n (f ), δ n (f )) ∈ R n×n . The conditional distribution can be computed using standard finite-dimensional formulae as Since the components of c are pairwise independent, we have K cc = Λ = diag (V(c 1 ), . . . , V(c m )). Furthermore, since δ n is a vector of linear functionals, we have, for each i ∈ {1, . . . , n}, that δ i f = m j=1 c j δ i φ j . Therefore, we have and so K δc = (δΦ)Λ. Thus all the required quantities can be explicitly evaluated.

E.2 Sampling
To sample from the posterior process f (·) | δ n (f), we can sample from the conditional distribution c | δ n (f) and then utilise the basis expansion of f in (D.2). To achieve this, we are required to perform a matrix square root of the posterior covariance matrix Σ δ , and we recall that, when conditioning on exact information, the resulting Σ δ is singular in general. The standard solution of performing a singular value decomposition (SVD) is unsuitable, since the Σ δ often have repeated singular values, which are incompatible with existing implementations of automatic differentiation that assume uniqueness of the singular values (Papadopoulo and Lourakis, 2000;Paszke et al., 2019). Although there have been recent efforts to address this (Wang et al., 2021b), the resulting algorithms are computationally prohibitive in our setting. An alternative method of sampling from f (·) | δ n (f) is called Matheron's update rule (Wilson et al., 2021, Corollary 4). Matheron's update rule takes the form The advantage of Matheron's update rule over the preceding approach is that we are not required to compute the square root of Σ δ ; this is the default approach used in GaussED.

E.3 Optimising the Acquisition Function
As discussed in Section 2.4, we utilise stochastic optimisation methodology to optimise the acquisition function. Unfortunately, the acquisition functions often exhibit multiple local optima, implying that it is unlikely that the optimiser will find a global optima. There are many approaches to reduce this probability, for instance by running the optimiser at different initialisations in parallel. In GaussED, the default approach is to sample uniformly from the design set, then evaluate the acquisition function at each of the sample points, before proceeding to initialise the optimiser at the best obtained point (i.e. Monte Carlo optimisation is used to initialise a stochastic optimisation method). This was the approach used in all the experiments of Section 3. Since our design sets are based on intervals 8 , we perform a standard reparameterisation to obtain a global optimisation problem in R d . This is achieved by using a scaled logistic function of the form where, for x, a, b ∈ R d , we consider logit to be applied component-wise.

F Experimental Details
In this section we present full details for the experiments presented in Section 3. All experiments can be reproduced using source code available at https://github.com/MatthewAlexanderFisher/GaussED.

F.1 Probabilistic Solution of PDEs
Approximating the Loss: Following from Section 3.1, recall that the quantity of interest was the function f, implying the loss takes the form Since there is not a closed-form solution to this integral when g is a Gaussian process, we proceed by approximating the integral through a cubature rule. For this experiment, we performed a Riemann sum over a uniform 15 × 15 grid over the domain X = [−1, 1] 2 .
Gaussian Model: For this experiment, we used a mean-zero Gaussian process f with Matérn covariance function with smoothness parameter ν = 3.5. The Dirichlet boundary conditions of the PDE were automatically enforced by the spectral GP approximation, applied to the domain X = [−1, 1] 2 (c.f. Equation (16)).
Optimisation: For both the optimisation of the acquisition function and performing maximum likelihood estimation, we used the Adam stochastic optimisation methodology (Kingma and Ba, 2015).
Using the methodology discussed in Appendix E.3, at each iteration of SED, we sampled 100 points uniformly from the design set and computed the corresponding values of acquisition function, using the default values of N = 81 and M = 9 in the stochastic gradient estimator of Section 2.4. We then proceeded by initialising the stochastic optimiser at the sample point which minimised the acquisition function. The learning rate used was the default value of 10 −1 and the optimiser was run for 1000 iterations, at each step of SED.
Using the methodology as discussed in Section 2.6, we began optimising the amplitude λ and the lengthscale after n 0 = 10 iterations of SED. This n 0 = 10 is the default value in GaussED. The initial parameter values were taken as the default values of λ = 1 and = 0.2. The learning rate used was the default value of 10 −3 and the optimiser was run 1000 iterations, at each step of SED.
Code: The code used to run the experiment can be seen in Figure 1 and discussed in Section 3.1.

F.2 Tomographic Reconstruction
Approximating the Loss: Following from Section 3.2, recall that the quantity of interest was the function exp(3f), implying the loss takes the form We follow the same approach of Appendix F.1 and approximate the integral through a Riemann sum, now over a uniform 25 × 25 grid over the domain X = [−1, 1] 2 .
Gaussian Model: For this experiment, we utilised a stationary Gaussian process model f with Matérn covariance with smoothness parameter ν = 2. The Gaussian model is defined on the domain [−1.05, 1.05] 2 , since the boundary conditions of the resulting GP do not necessarily agree with the boundary conditions of the quantity of interest.
Quantity of Interest: Recall from Section 3.2 that the quantity of interest was of the form Since this quantity of interest defines a circle within the domain X = [−1, 1] 2 , it is possible to find a closed form solution to the line integrals of f for given parameters values (θ, x, y). However, for ease of implementation and to allow our approach to be easily generalised to more complex examples, we computed the line integrals of f by performing a Riemann integral over a uniform mesh consisting of 200 evaluations from f.

Optimisation:
All settings used were the same as the previous experiment detailed in Appendix F.1, apart from the following settings: We began optimising the amplitude λ and lengthscale parameters of the Gaussian process at step n 0 = 0. The initial parameter values were taken as λ = 0.5 and = 0.4.

Code:
The GaussED code used to run this experiment is presented in Figure 5. The structure of the code is quite different to the code used in the other experiments (Figure 1 and Figure 7). This is due to the fact that the design object (d) is not instantiated by the EvaluationDesign class. Note that for both the PDE experiment (Section 3.1) and the Bayesian optimisation experiment (Section 3.3), the design sets D consisted of evaluations of the Gaussian process, δ(f ) = f (x), or its derivatives δ(f ) = ∂ i f (x). In situations such as these, the EvaluationDesign class may be used. For this example, however, the observed data consists of line integrals. Therefore, in this more general situation, we must specify two further functions: Given a parameterisation D θ of the design set, the first function must take in a sequence of parameters θ 1 , . . . , θ n and return the corresponding (δΦ) ij = δ θi φ j matrix, where the φ j are the eigenfunctions of (17). This is reflected in the code ( Figure 5) in the function d_func, which, for each design set parameter constructs the corresponding line integral for a given number of basis functions (m). The second function we must specify must be able to, given a parameter θ, sample from the process δ θ f | δ n , where δ n is data gathered from SED. This is directly reflected in the code ( Figure 5) in the function d_sample. Note that, in Figure 5 we omit the get_line_int_gps function. This is a function that, given a parameter value θ and Gaussian process f , returns the corresponding δ θ f object. We do this because get_line_int_gps is complexified due to the parameterisation of the line function r(x) and the calculation of the limits of integration a, b in the line integral b a f (r(x)) dx.
We, therefore, omit get_line_int_gps for clarity. A second major difference, is the use of an output warp (OutputWarp). Due to the non-linear nature of the output warp, the resulting object qoi is only able to sample from the prior and posterior. Note that the syntax for specifying a output deformation of a GP is the same as specifying other transformations (e.g. see Figure 1 and Figure 7).
Another difference is that the Gaussian model specified in the PDE experimental code (Figure 1) agreed with the boundary conditions of the PDE; here, however, we specify the domain (gp.set_domain) as [−1.05, 1.05] × [−1.05, 1.05]. Since we took the domain of the Gaussian process to be larger than the domain on which the task is defined, we must also specify the domain of the design object (d.set_domain), which otherwise, by default, would be taken as the same the Gaussian model (gp).
Finally, note that the acquisition function (acq), as discussed previously, is instantiated with a nugget value of 10 −2 and the experiment object (experiment) is instantiated with m = 28 2 basis functions. This is in contrast to the code for the PDE example (Figure 1), which used the default value of m = 30 2 basis functions.

F.3 Gradient-Based Bayesian Optimisation
Approximating the Loss: Recall from Section 3.3 that our quantity of interest is q(f) = max x∈X log L(x). Thus, our loss function takes the form In order to optimise the samples, we used a grid-based optimiser using a uniform 40 × 40 grid over the domain of interest [0.45, 0.9] × [0.09, 0.5].
Gaussian Model: For this experiment, we used a mean-zero stationary Gaussian model f with Matérn covariance, with smoothness parameter ν = 3. Since our GP satisfies the boundary conditions in (16), which are unrelated to the task at hand, we took the domain of the GP to be [0.4, 0.95] × [0.04, 0.55], which is wider than the domain on which the task is defined.

Optimisation:
All settings were as the previous experiment detailed in Appendix F.1, apart from the following settings: For both the optimisation of the acquisition function and performing maximum likelihood estimation, we used the Adam stochastic optimisation methodology (Kingma and Ba, 2015).
Using the methodology as discussed in Appendix E.3, at each iteration of SED, we sampled 100 points times uniformly from the design set and computed the corresponding values of acquisition function, using the default values of N = 81 and M = 9 in the stochastic gradient estimator of Section 2.4. We then proceeded by initialising the stochastic optimiser at the sample point which minimised the acquisition function. The learning rate used was the default value of 10 −1 and the optimiser was run 1000 iterations, at each step of SED. Furthermore, in order to increase the numerical stability of linear algebra operations, we used a nugget term of value 10 −5 .
For this experiment, we began optimising the amplitude λ and lengthscale at step n 0 = 10. The initial kernel parameter values were taken as the values of λ = 1 and = 0.1. The initial parameter values of the spatial deformation were taken as θ x = (1, 0, 0, 1) and θ y = (1, 0, 0, 1), thus specifying the initial spatial deformation as the identity function. The learning rate used was the default value of 10 −3 and the optimiser was run 1000 iterations, at each step of SED. Code: The GaussED code used to run this experiment is presented in Figure 7. The structure of the program is very similar in nature to the PDE experiment of Section 3.1. The first difference is that, at each step of SED, we evaluate multiple functionals δ from the design D. This is directly reflected in Figure 7, where the design object (d) is constructed by the statistical model f and its first derivatives ([gp, gp_d1, gp_d2]).
The second difference is that, at each step SED, we perform a maximisation, rather than an integral, of sample paths when estimating the acquisition function. In the code for the PDE experiment ( Figure 1) the integral of posterior samples is hidden within the loss object (L2(qoi)), which, by default, performs a Riemann sum over a uniform mesh if the quantity of interest qoi is function valued. Therefore, in Figure 7 we specify a numerical method that acts on samples from f . In this instance, we perform a grid search (maximise_method) over a uniform 40 × 40 mesh (mesh) over the domain of optimisation.
Another difference is that the Gaussian model specified in the PDE experimental code (Figure 1), agrees with the boundary conditions of the PDE and therefore the domain of the GP is taken as the default value [−1, 1] 2 . In Figure 7, we must specify the domain (gp.set_domain) as [0.4, 0.95] × [0.04, 0.55]. Since we took the domain of the Gaussian process to be larger than the domain over which we wish to maximise, we must also specify the domain of the design object (d.set_domain), which otherwise, by default, would be taken as the same the Gaussian model (gp).
The final difference is that, in order to increase the numeric stability of linear algebra operations in the SED, we specify a nugget term (nugget) of value 10 −5 in the acquisition function acq.

G Evaluating Computational Aspects of GaussED
In this section we empirically investigate computational aspects of GaussED. In Appendix G.1, we explore the role of the optimisation methodology and how this affects the experimental design as well as the quality of output. In Appendix G.2, we investigate how the number of basis functions used, for a given problem, affects the quality of posterior inference.

G.1 Investigating the Efficacy of Stochastic Optimisation
In this section, we investigate the effect of the random seed on the quality of the experimental design and, further, investigate the effect of changing the stochastic optimisation approach itself. To explore these aspects of GaussED, we repeat the Bayesian optimisation with gradient data experiment presented in Section 3.3. Recall that, in all the demonstrations in Section 3, we utilised the Adam stochastic optimisation method (Kingma and Ba, 2015).
Results on the effect of the random seed can be seen in Figure 8 and Figure 9. The obtained designs imply that our approach of SED is sensitive to the initial conditions. Although the specific design is sensitive, the overall performance and qualitative nature of the designs are approximately independent of random seed.
Results on the effect of stochastic optimisation methodology can be seen in Figure 10 and Figure 11. In each of these experiments, the random seed was fixed, and so we are only comparing the effect of different optimisation methodologies. In each experiment, the learning rate was set at 10 −1 and the other parameter values were taken as their default values, as specified in PyTorch (Paszke et al., 2019).

G.2 Investigating the Effect of the Number of Basis Functions
Picking an appropriate number of basis functions for a given problem is an important means to reduce computational cost in GaussED. In this section, we investigate how the number of basis functions may affect the quality of posterior inference. To this end, it is sufficient to consider the behaviour of posterior sampling in dimension d = 1, since the behaviour will naturally extend to higher-dimensions due to the exponential scaling of the number of basis function due to (17).
In the event where the number of basis functions is smaller than the number of linearly independent data, the resulting posterior will not be well-defined in general. The introduction of a nugget term on the diagonal of the covariance matrix, implicitly assuming noisy Gaussian observations, is a pragmatic solution that is widely-used. However, the success of this strategy depends crucially on an appropriate amount of regularisation being introduced.
Results on the effect on the number of basis functions and the nugget term are presented in Figure 12. Through visual inspection, by m = 20 basis functions, it appears that the posterior process has converged  In the event where m is smaller than the number of data points conditioned upon, the corresponding nugget term σ 2 is also displayed.