{'fraction_non_alphanumeric': 0.0324254215304799, 'fraction_numerical': 0.005590589919048258, 'mean_word_length': 4.548387096774194, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse Spectrum Gaussian Process for Bayesian Optimisation

We propose a novel sparse spectrum approximation of Gaussian process (GP) tailored for Bayesian optimisation. Whilst the current sparse spectrum methods provide good approximations for regression problems, it is observed that this particular form of sparse approximations generates an overconfident GP, i.e. it predicts less variance than the original GP. Since the balance between predictive mean and the predictive variance is a key determinant in the success of Bayesian optimisation, the current sparse spectrum methods are less suitable. We derive a regularised marginal likelihood for finding the optimal frequencies in optimisation problems. The regulariser trades the accuracy in the model fitting with the targeted increase in the variance of the resultant GP. We first consider the entropy of the distribution over the maxima as the regulariser that needs to be maximised. Later we show that the Expected Improvement acquisition function can also be used as a proxy for that, thus making the optimisation less computationally expensive. Experiments show an increase in the Bayesian optimisation convergence rate over the vanilla sparse spectrum method.


Introduction
Bayesian optimisation is a leading method for global optimisation for expensive black-box functions [1,2,3]. It is widely used in hyperparameter tuning of massive neural networks [4], some of which can take days to train, It has also been used for optimisation of physical products and processes [5] where one experiment can take days, and experiments can also be expensive in terms of materials cost. However, in many scenarios it is possible that the objective function is not very costly. For example, consider the use cases involving simulation software. They are often used in the early of a product design process to reduce a massive search space to a manageable one before real products are made. Whilst evaluation, a few thousands time may be feasible, but evaluating a million time is not. We term this problem as semi-expensive optimisation problems. Such problems cannot be handled by the traditional global optimisers. Bayesian optimisation will also struggle, because its main ingredient, Gaussian process, does not scale well beyond few hundred observations, rendering it incapable to exploit the semi-expensive nature of the objective functions. In this paper we address the scalability issue of the Gaussian processes for Bayesian optimisation for semi-expensive objective functions.
Scalability issues for Bayesian optimisation has been previously addressed in two main ways: 1) by replacing Gaussian process with a more scalable Bayesian model, e.g. using Bayesian neural network [6] or Random Forests [7], or 2) by making sparse approximations for the full Gaussian process. The latter is often desirable as it provides the option of using many different types of covariance functions, offering a flexible approach to function modeling. There are many sparse models in literature such as fully independent training conditional approximation (FITC) [8,9], variational free energy (VFE) [10] and sparse spectrum Gaussian process (SSGP) [11]. FITC is a method of systematically approximating the full GP kernel function, and thus it corresponds to exact inference. VFE is derived from a variational approach where one can obtain a lower bound of the marginal likelihood by minimising the KL divergence between the sparse model and the true model. SSGP is a sparse Bayesian linear regression model of pairs of trigonometric basis functions and is a only sparse model so far that induces a stationary covariance function so that we can sample the SSGP posteriors over its global maxima. Hence, we use SSGP as our baseline in this work. These models perform well on regression tasks, but not on optimisation tasks, because the state of art models focus on approximating the full GP by using sparse sets (generally known as inducing inputs). It works for the regression tasks where overall accuracy is important, however, it may fail for Bayesian optimisation, as it may require more accurate models in important regions (high function values or high uncertainty). The method proposed in [12] is the only work that tried to model a sparse GP more suited to BO. Their method use weighted KL divergence to maintain fixed size inducing points in an online fashion by removing the worst performing points. By using function values as an weight, they ensure that more inducing points are selected in the high function value regions, compared to the low function value regions. However, the search of Bayesian optimisation is also guided by the variance, so sometimes when the function value regions are low but it is under explored, then it can be a good location to try. Unfortunately, their method does not consider such cases. Hence, developing a sparse GP framework which is tuned to the goal of Bayesian optimisation is still open.
In this paper, we propose a novel modification to the sparse spectrum Gaussian process approach to make it more suitable for Bayesian optimisation. We utilise the enhanced SSGP marginal likelihood to select our sparse set. This regulariser allows us to minimise the differences between two distributions that are global maxima of the full GP and that of the sparse model. The form of the regulariser is selected as the entropy of the maxima distribution (in the BO community, it is often known as p(x * )). Since without this regulariser the original formulation results in a more confident GP, we maximise the entropy of p(x * ). We also derive the formula of using Thompson sampling in SSGP to estimate the global maxima of the sparse model. However, computing p(x * ) can be very expensive. Hence, as an alternative we propose to use Expected Improvement acquisition function, which often follows the shape as p(x * ), as it's proxy. We demonstrate our method on four synthetic functions and applications on hyperparameter tuning for machine learning algorithms as well as alloy optimisation. In all experiments our method empirically demonstrates the superior performance of optimisation over standard Sparse spectrum methods.

Background
In this paper we consider global optimisation problems in the form (1) where f : x → R and X is a compact subspace in the d dimensional real space R d .

Bayesian Optimisation
There are two main components in Bayesian optimisation. The first component is to model the function using Gaussian process as a prior. The other one is to search the best location to perform the new experiment. Acquisition function is the surrogate utility function that guides the search for the new point.
Gaussian Process Gaussian process (GP) [13] is a distribution over functions and it is specified by its mean function µ(x) and covariance function k(x, x ′ ), where x ∈ R d is the input space. A sample from a Gaussian process is a function Without loss of generality, the prior mean function can be assumed as zero function, so that the Gaussian process can be fully defined by the k(x, x ′ ). The squared exponential (SE) function is a popular choice of kernel function where ρ l is the kernel length scale, and σ 2 f is the signal variance. Given a set of observations , the joint probability distribution of observations D and a new point {x t+1 , f t+1 } is given as , and K is the kernel matrix presented as The predictive distribution of f t+1 follows a normal distribution with the predicted mean and variance derived as If the observations are contaminated with zero-mean Gaussian noise , the the predicted mean and variance can be represented as GP is prohibitive for large data sets because computational complexity of prediction is O(t 3 ) due to the need for the inversion of the covariance matrix, where t is the number of function observations. Acquisition Functions Acquisition function known as its effective capability of trading off exploration and exploitation in Bayesian optimisation. Exploration means sampling in the areas where the epistemic uncertainty about the function values are high. Exploitation means sampling in the areas where the mean prediction for function values are high. In this paper, we use Expected Improvement (EI) [14] as the criteria for our acquisition function. Assume that the current maximum is f (x + ). The improvement function I(x) is written as: The analytic form of E(I(x)) can be obtained as . Φ(z) and φ(z) are the CDF and PDF of standard normal distribution.

Proposed Framework
Directly Applying sparse GP to optimization problems might not receive promised performance. This is because that the challenge of using sparse GP for Bayesian optimisation is different from GP regression. More specifically, the most of existing sparse models focus on approximating the full GP as closely as possible by using minimum resources . This motivation is desirable for regression tasks, however, has significantly drawbacks in optimisation tasks as the limited resources of the sparse models may be allocated to closely model regions of parameter space that perform poorly (over-fitting) and are therefore less important for optimisation. As a result, we propose a novel sparse spectrum Gaussian process model tailoring Bayesian optimization. Our approach to sparse spectrum GP is similar to that of [SSGP], but utilises a regularisation term in marginal likelihood to select our optimum features (sparse set). This regularisation term allows us to minimise the difference between full GP p(x * ) distribution and SSGP p e (x * ) distribution, so that the model can better allocate the limited modeling capacity of the sparse GP to further the goal of optimisation. We provide two different forms of regularisation: one with the entropy of the maxima distribution, p e (x * ), and another slightly approximate but computationally efficient version by treating expected improvement acquisition function as a proxy for p e (x * ). We start with the formulation of naive sparse spectrum GP and then develop the formulation for enhanced sparse spectrum GP.

Bayesian Optimisation using Naive Sparse Spectrum GP
We have mentioned sparse spectrum GP (SSGP) in the first section and claimed that we are using it as our baseline. We will introduce details of it in this part. Sparse spectrum GP is a spare approximation that sparsifies power spectral density of a generic GP with stationary covariance function. For a stationary GP, the power is equal to the prior variance k(x, x) = k(0) = σ 2 f . The frequency vector s has the same length D as the input vector x. The d-th element of s can be interpreted as the frequency associated to the d-th input dimension. The Wiener-Khintchine theorem states that the power spectrum and the auto correlation of the random process constitute a Fourier pair. In our case, given that f (·) is drawn from a stationary Gaussian process, the autocorrelation function is equal to the stationary covariance function, and we have: where k(x i , x j ) denotes the generic covariance function for a stationary Gaussian process in the time domain, while S(s) denotes the power spectral density of a stationary random process expresses how the power is distributed over the frequency domain and τ = x i − x j is a stationary covariance function.
In Bochner's theorem, it has stated that any stationary covariance function k(τ ) can be represented as the Fourier transform of a positive finite measure []. This means that the power spectrum S(s) is a positive finite measure, and in particular that it is proportional to a probability measure, S(s) ∝ ps(s). By evaluating the covariance function above at τ = 0 we can obtain the proportionality constant directly as: Since S(s) is proportional to a multivariate probability density in s, so we can rewrite the covariance function as an expectation: where the superscript asterisk denotes complex conjugation. To be noted, this expectation is an exact expansion of the generic covariance function in the form of distribution over frequencies. By using simple Monte Carlo we can approximate this integral by obtaining an average of a few samples corresponding to a finite set of frequencies, which also known as spectral points. A valid Monte Carlo procedure is to sample frequencies always as a pair {s r , −s r }as the fact that the power spectrum is symmetric around zero. This also helps of preserving the property of the exact expansion of covariance function, and now the expansion is : [e 2πis T r xi (e 2πis T r xj ) * + (e 2πis T r xi ) * e 2πis T r xj ] where m is the number of spectral frequencies (features),s r is drawn from ps(s) and Re[·] denotes the real part of a complex number. Until now the only detail left of spare spectrum GP is how do we select the m number of optimal spectral frequencies. By noticing k( we can then rewrite the approximated covariance function as: where we define the column vector of length 2m containing the evaluation of the m pairs of trigonometric functions at x: Then the log marginal likelihood can be given as All hyperparameters include m are denoted by Θ. We can then learn the optimal features m as well as all kernel hyperparameters by maxmising the log marginal likelihood. Once the optimal features have been obtained, the posterior mean and variance can be estimated and then BO can be performed. The posterior mean and variance for sparse spectrum GP can write as: We have discussed before, the SSGP is designed for the regression task. For the goal of optimisation, we expect the knowledge about the global maxima from the sparse model is as close to that from the full GP as possible. An intuitive representation for the knowledge about the global optimum is the posterior distribution of the global optimum, denoted as p(x * ). We illustrate the difference of p s (x * ) from the SSGP (we talk how to obtain it later) and p(x * ) from the full GP by using the example in Fig 1a-1b. We can see that p s (x * ) is more narrow than p(x * ). In other words, p s (x * ) has small entropy than p(x * ) . The entropy calculation of both distributions also numerically proved the visualised results and we show them at right side of each plot (middle). Hence, we develop a new sparse GP model which can resemble more closely to the knowledge of the p(x * ) as that of the full GP. , EI function (bottom) and global maxima distributions (middle) of (a) Full GP (p(x * )) (b) SSGP with 30 optimal features (p s (x * )) (c) Enhanced SSGP with 30 optimal features (p e (x * )).

Bayesian Optimisation using Enhanced Sparse Spectrum GP
We denote the posterior distribution of the global maxima from our method as p e (x * ). Following the discussion above, a straightforward idea is to minimise the difference between two distributions p e (x * ) and p(x * ) and KL divergence is an usual way to measure the difference, or D KL (p e (x * ) || p(x * )). However, we cannot access p(x * ) due to the unavailability of the full GP so we cannot compute the KL divergence above. Note that p(x * ) is fixed given the full GP. Therefore, we can simplify our task to the problem of making the entropy of p e (x * ) close to that of p(x * ). A subsequent question is whether we should maxmise or minimise the entropy of p e (x * ). It has been demonstrated in the original work that the SSGP tends to overfit. This overfitting property of SSGP indicates that the entropy for the posterior distribution of the global maxima in the sparse spectrum GP is less than that of the full GP. Therefore, we will maximise the entropy of p e (x * ) along with the fitting term for finding the optimal frequencies (sparse set) and the hyperparameters of the sparse spectrum GP. Based on the idea above, we propose an enhanced sparse spectrum Gaussian process (Enhanced SSGP). The loss function is defined as where the first term is the log marginal likelihood in SSGP, the second term is a regularisation term consistently the entropy for the posterior distribution of the global maxima p e (x * ) and λ is the trade-off parameter. We obtain the optimal features w by maximising the Eq.(2) as w = argmax(logp(y 1:t |Θ) + λlogH(p e (x * ))) Finally we expand the Eq.(2) and maxising the following loss function Similar to the SSGP, the log marginal likelihood of our enhanced SSGP is a function of w. The p e (x * ) is dependent on w through the GP. Next, we propose a method for Thompson sampling in the sparse spectrum GP to compute p e (x * ).
Since SSGP framework provides the formulation using random features φ(x) so it is amenable to Thompson sampling [15]. Now we can consider a linear model f (x) = φ(x) T θ where θ ∽ N (0, I) has a standard Gaussian distribution and observations . Then the posterior of θ given (D n , φ(x)) is a normal N (A −1 Φ T y, A −1 σ 2 noise ), where A and Φ have already been defined in sparse spectrum GP.
We then let φ (i) and θ (i) be a random set of features and corresponding posterior weights sampled both according to the generative process given above and they can be used to construct the function g (i) (x) = φ (i) (x) T θ (i) , which is an approximate posterior sample of f . We then maximise this function to obtain x (i) * = argmax x∈X g (i) (x), which is approximately distributed according to p e (x * | D n ). Once we have generated sufficient samples, we can use kernel density estimation (KDE) method to estimate the probability density function of samples which is a representation of p e (x * ) distribution. We visualise our p e (x * ) in Fig 1c. Expected Improvement Acquisition Function as a Proxy Now we know that changing w can result in different p e (x * ). However, p e (x * ) has no analytic form so that optimising the loss function presents heavy computation burden. To reduce the computation, we propose to use EI as a proxy for p e (x * ). This choice is reasonable since they both measure the belief about the location of the global optimum. These can also be observed from Fig 1 and we visualise the result of using EI as a proxy in Fig 2where the entropy of EI is calculated through kernel density estimation. To end up, we can obtain the optimal features by maximising the loss function below

Experiments
We evaluate our method on applications of Bayesian optimisation. We compare the proposed method enhanced SSGP with the following baselines for optimisation: -Full Gaussian process (Full GP) -Sparse spectrum Gaussian process (SSGP) We use four benchmark functions and two real world problems, namely alloy optimisation using a thermodynamic simulator and hyperparemeter tuning as our test bed.

Experiment with Synthetic Functions
The benchmark functions we compare are: In all settings, we use EI as acquisition function and use the DIRECT [16] optimiser to optimise EI. We run each method 50 trials with different initialisation and report the simple regret along with standard errors. Simple regret is defined as is the global maxima and f (x + ) = max x∈{x1:t} f (x) is the best value so far. Fig 3 plots the simple regret vs iterations for optimising 2D and 5D Ackley functions, 3D Rosenbrock function and 6D Hartmann function. In Fig 3a 2D Ackley function , we use 20 observations and the number of optimal features for both SSGP and our method is set to 20. In terms of kernel parameters, we use the isotropic length scale σ l = 0.5, signal variance σ 2 f = 2 and noise variance σ 2 noise = (0.01) 2 for all three methods.   Bayesian optimisation with Full GP with SE kernel performs the best for 2D Ackley function, and it is easy to explain since all observations have been used in the algorithm. Our approach performs better than SSGP as it reaches closer to the global optimum. In Fig 3b 5D Ackley example, we use 50 observations and 20 optimal features. The isotropic length scale, signal variance and noise variance are set as σ l = 0.1, σ 2 f = 2, σ 2 noise = (0.01) 2 respectively. We can see that our method performs the best at the end. It is explainable since we only use 20 optimal features so the inversion of the covariance matrix is more accurate than the inversion of Full GP's. We also receive similar results for 3D Rosenbrock function and 6D Hartmann function in Fig 3c and Fig 3d respectively.

Alloy Optimisation
In this experiment we aim to design an alloy through a thermodynamic simulator called ThermoCalc [17]. Given a composition of an alloy the simulator can compute thermodynamic equillibrium and predict the micro-structure of the resultant alloy using CALPHAD [18] methodology. In a joint project with our metallurgist collaborators we aim to design an alloy with a micro-structure that contains as much fraction of FCC phases as possible. The search space was defined in a 15 dimensional space of the following elements: Fe, Ni, Cr, Ti, Co, Al, Mn, Cu, Si, Nb, Mo, W, Ta, C, N. For each composition we measure the amount of FCC in fraction, and we want it to maximise it. Since each composition takes around 10 minutes of computation to measure the amount of FCC, it fits perfectly in our definition of semi-expensive functions. In our experiment we use 50 initial points and 20 sparse features. The results in Fig 4c shows convergence of our method compared to the vanilla sparse spectrum Gaussian process based approach. The difference is more pronounced with higher number of observations. We were able to run only once due to time limitation.

Hyperparameter Tuning
We experiment with two real world datasets for tuning hyperparameters of Support vector machines (SVM) classifier. We optimise two hyperparameters in SVM which are the cost parameter (C) and the width of the RBF kernel (γ). The search bounds for the two hyperparameters are C =10 λ where λ ∈ [−3, 3] and γ = 10 ω where ω ∈ [−3, 0] correspondingly. To make our search bounds manageable, we optimise for λ and ω. The two datasets LiverDisorders and BreastCancer are available from UCI data repository [19]. We run each method with 50 trials different initialisation, and the results are plotted in

Conclusion
In this paper we proposed an enhancement to the sparse spectrum Gaussian process method to make it more suitable for Bayesian optimisation application. The original formulation results in an over-confident Gaussian process. Bayesian optimisation using such an over-confident Gaussian process may fare poorly as