Stable Bayesian Optimisation via Direct Stability Quantification

TL;DR

This paper introduces a Bayesian optimization method that quantifies and optimizes for stability in expensive functions, ensuring solutions are robust to input variations in physical and industrial processes.

Contribution

It develops a novel approach using gradient Gaussian Processes to estimate stability probabilities, guiding optimization toward stable maxima and filtering results for robustness.

Findings

Successfully finds stable maxima in synthetic problems

Effectively applies to real-world industrial scenarios

Outperforms traditional optimization methods in stability criteria

Abstract

In this paper we consider the problem of finding stable maxima of expensive (to evaluate) functions. We are motivated by the optimisation of physical and industrial processes where, for some input ranges, small and unavoidable variations in inputs lead to unacceptably large variation in outputs. Our approach uses multiple gradient Gaussian Process models to estimate the probability that worst-case output variation for specified input perturbation exceeded the desired maxima, and these probabilities are then used to (a) guide the optimisation process toward solutions satisfying our stability criteria and (b) post-filter results to find the best stable solution. We exhibit our algorithm on synthetic and real-world problems and demonstrate that it is able to effectively find stable maxima.

Figures

Tables

Table 1. Table 1: Relevant standard isotropic kernels. In this table K ​ ( 𝐱 , 𝐱 ′ ) = κ ​ ( 1 2 ​ ‖ 𝐱 − 𝐱 ′ ‖ 2 2 ) 𝐾 𝐱 superscript 𝐱 ′ 𝜅 1 2 superscript subscript norm 𝐱 superscript 𝐱 ′ 2 2 K({\bf x},{\bf x}^{\prime})=\kappa(\frac{1}{2}\|{\bf x}-{\bf x}^{\prime}\|_{2}^{2}) , s 𝑠 s is the differentiability of κ 𝜅 \kappa , κ ( q ) ​ ( r ) = ∂ q ∂ r q ​ κ ​ ( r ) superscript 𝜅 𝑞 𝑟 superscript 𝑞 superscript 𝑟 𝑞 𝜅 𝑟 \kappa^{(q)}(r)=\frac{\partial^{q}}{\partial r^{q}}\kappa(r) is the q th superscript 𝑞 th q^{\rm th} derivative ( q ∈ ℤ s + 1 𝑞 subscript ℤ 𝑠 1 q\in{\mathbb{Z}}_{s+1} ) of κ 𝜅 \kappa , L ↑ , L ↓ superscript 𝐿 ↑ superscript 𝐿 ↓ L^{\uparrow},L^{\downarrow} relate to the effective length-scales, and Δ r subscript Δ 𝑟 \Delta_{r} , Δ Δ \Delta are the Taylor remainder bounds (bounds labelled ∗ may be calculated numerically using ( 5 ) and ( 6 )).

Kernel	Derivatives	$s$
RBF	$\begin{array}{c}{\kappa }_{\left(\gamma \right)}\left(r\right)=e^{-\frac{1}{{\gamma }^2}r}\\ {\kappa }_{\left(\gamma \right)}^{\left(q\right)}\left(r\right)={\left(-\frac{1}{{\gamma }^2}\right)}^q{e}^{-\frac{1}{{\gamma }^2}r}\end{array}$	$\infty$	$\begin{array}{c}{L}_{\left(\gamma \right)}^{\uparrow }=L_{\left(\gamma \right)}^{\downarrow }=\frac{1}{{\gamma }^2}\\ {\Delta }_{\left(\gamma \right)r}\left(\delta r\right)={\Delta }_{\left(\gamma \right)}\left(\delta r\right)=0\end{array}$
$\frac{1}{2}$-Matern	$\begin{array}{c}{\kappa }_{(\frac{1}{2},\rho )}\left(r\right)=e^{-\frac{\sqrt{2r}}{\rho }}\end{array}$	$0$	$\begin{array}{c}{L}_{(\frac{1}{2},\rho )}^{\uparrow }=\sqrt{\frac{1}{2}}\frac{1}{2\rho }\\ L_{(\frac{1}{2},\rho )}^{\downarrow }=0.3764\sqrt{\frac{1}{2}}\frac{1}{2\rho }\\ {\Delta }_{(\frac{1}{2},\rho )r}\left(\delta r\right),{\Delta }_{(\frac{1}{2},\rho )}\left(\delta r\right)={}^{\ast }\end{array}$
$\frac{3}{2}$-Matern	$\begin{array}{c}{\kappa }_{(\frac{3}{2},\rho )}\left(r\right)=\left(1+\frac{\sqrt{6r}}{\rho }\right)e^{-\frac{\sqrt{6r}}{\rho }}\\ {\kappa }_{(\frac{3}{2},\rho )}^{\left(1\right)}\left(r\right)=-\frac{\sqrt{3}}{\sqrt{2}\rho }{\kappa }_{(\frac{1}{2},\rho )}\left(r\right)\end{array}$	$1$	$\begin{array}{c}{L}_{(\frac{3}{2},\rho )}^{\uparrow }=\underset{c\in \{0,1\}}{max}\{1,\frac{{\kappa }_{d+\frac{1}{2}-c}\left(\frac{1}{2}{M}^2\right)}{{\kappa }_{d+\frac{1}{2}}\left(\frac{1}{2}{M}^2\right)}\}\sqrt{\frac{3}{2}}\frac{1}{2\rho }\\ L_{(\frac{3}{2},\rho )}^{\downarrow }=0.7528\sqrt{\frac{3}{2}}\frac{1}{2\rho }\\ {\Delta }_{(\frac{3}{2},\rho )r}\left(\delta r\right),{\Delta }_{(\frac{3}{2},\rho )}\left(\delta r\right)={}^{\ast }\end{array}$
$\frac{5}{2}$-Matern	$\begin{array}{c}{\kappa }_{(\frac{5}{2},\rho )}\left(r\right)=\left(1+\frac{\sqrt{10r}}{\rho }+\frac{10r}{3{\rho }^2}\right)e^{-\frac{\sqrt{10r}}{\rho }}\\ {\kappa }_{(\frac{5}{2},\rho )}^{\left(1\right)}\left(r\right)=-\frac{\sqrt{5}}{3\sqrt{2}\rho }{\kappa }_{(\frac{3}{2},\rho )}\left(r\right)\\ {\kappa }_{(\frac{5}{2},\rho )}^{\left(2\right)}\left(r\right)=\frac{5}{6{\rho }^2}{\kappa }_{(\frac{1}{2},\rho )}\left(r\right)\end{array}$	$2$	$\begin{array}{c}{L}_{(\frac{5}{2},\rho )}^{\uparrow }=\underset{c\in \{0,1,2\}}{max}\{1,\frac{{\kappa }_{d+\frac{1}{2}-c}\left(\frac{1}{2}{M}^2\right)}{{\kappa }_{d+\frac{1}{2}}\left(\frac{1}{2}{M}^2\right)}\}\sqrt{\frac{5}{2}}\frac{1}{2\rho }\\ L_{(\frac{5}{2},\rho )}^{\downarrow }=0.5018\sqrt{\frac{5}{2}}\frac{1}{2\rho }\\ {\Delta }_{(\frac{5}{2},\rho )r}\left(\delta r\right),{\Delta }_{(\frac{5}{2},\rho )}\left(\delta r\right)={}^{\ast }\end{array}$

Table 2. Table 2: Lower bounds β d subscript 𝛽 𝑑 \beta_{d} on Matérn kernel ratios R d ​ ( r ) = κ ( d + 1 2 , ρ ) ​ ( r ) κ ( d + 3 2 , ρ ) ​ ( r ) subscript 𝑅 𝑑 𝑟 subscript 𝜅 𝑑 1 2 𝜌 𝑟 subscript 𝜅 𝑑 3 2 𝜌 𝑟 R_{d}(r)=\frac{\kappa_{(d+\frac{1}{2},\rho)}(r)}{\kappa_{(d+\frac{3}{2},\rho)}(r)} .

$d$	${\beta }_d$	$d$	${\beta }_d$	$d$	${\beta }_d$
0	0.752871	22	0.999107	44	0.99976
1	0.92244	23	0.999178	45	0.99977
2	0.96113	24	0.999241	46	0.99978
3	0.976487	25	0.999297	47	0.999789
4	0.984203	26	0.999347	48	0.999797
5	0.988643	27	0.999391	49	0.999805
6	0.991437	28	0.999432	50	0.999813
7	0.993311	29	0.999468	51	0.99982
8	0.99463	30	0.999501	52	0.999826
9	0.995593	31	0.999531	53	0.999833
10	0.996318	32	0.999559	54	0.999839
11	0.996878	33	0.999584	55	0.999844
12	0.997319	34	0.999607	56	0.99985
13	0.997673	35	0.999628	57	0.999855
14	0.997961	36	0.999648	58	0.99986
15	0.998198	37	0.999666	59	0.999864
16	0.998397	38	0.999682	60	0.999869
17	0.998564	39	0.999698	61	0.999873
18	0.998706	40	0.999712	62	0.999877
19	0.998829	41	0.999725	63	0.99988
20	0.998934	42	0.999738	64	0.999884
21	0.999026	43	0.999749	65	0.999887

Equations

$$\begin{array}[]{rl}m_{\mathbb{D}}\!\left({\bf x}\right)&\!\!\!\!={\bf k}^{{\rm T}}\!\left({\bf x}\right)\left({\bf K}+{\sigma}^{2}{\bf I}\right)^{-1}\!{\bf y}\\ \!\!\!\!\!\lambda_{\mathbb{D}}\!\left({\bf x},{\bf x}^{\prime}\right)&\!\!\!\!=K\!\left({\bf x},{\bf x}^{\prime}\right)-{\bf k}^{{\rm T}}\!\left({\bf x}\right)\left({\bf K}+{\sigma}^{2}{\bf I}\right)^{-1}\!{\bf k}\left({\bf x}^{\prime}\right)\!\!\!\!\!\!\!\end{array}$$

$$\begin{array}[]{rl}{\bf m}_{\mathbb{D}}^{\!(q)}\!\left({\bf x}\right)&\!\!\!\!\!=\!\left({\nabla}_{\bf x}^{\otimes q}{\bf k}^{{\rm T}}\!\left({\bf x}\right)\right)\left({\bf K}+{\sigma}^{2}{\bf I}\right)^{-1}{\bf y}\\ \!\!\!\!\!\!{\mbox{\boldmath$\Lambda$}}_{\mathbb{D}}^{\!(q)}\!\left({\bf x},\!{\bf x}^{\prime}\right)&\!\!\!\!\!=\!\nabla_{\bf x}^{\otimes q}\nabla_{{\bf x}^{\prime}}^{\otimes q}{}^{{\rm T}}\!K\!\left({\bf x},\!{\bf x}^{\prime}\right)\!-\!\left({\nabla}_{\bf x}^{\otimes q}{\bf k}^{{\rm T}}\!\!\left({\bf x}\right)\right)\!({\bf K}+\!\!\!\!\!\!\\ &\!\!\!\ldots+{\sigma}^{2}{\bf I})^{-1}\!\left({\nabla}_{{\bf x}^{\prime}}^{\otimes q}{\bf k}^{{\rm T}}\!\left({\bf x}^{\prime}\right)\right)^{{\rm T}}\end{array}$$

$$\begin{array}[]{l}{\rm vec}\left(\nabla_{{\bf x}}^{\otimes q}\nabla_{{\bf x}^{\prime}}^{\otimes q{\rm T}}K\left({\bf x},{\bf x}^{\prime}\right)\right)=\left(\nabla_{{\bf x}^{\prime}}^{\otimes q}\otimes\nabla_{{\bf x}}^{\otimes q}\right)K\left({\bf x},{\bf x}^{\prime}\right)\end{array}$$

$$\begin{array}[]{l}K\left({\bf x},{\bf x}^{\prime}\right)=\kappa\left(\frac{1}{2}\left\|{\bf x}-{\bf x}^{\prime}\right\|_{2}^{2}\right)\end{array}$$

$$\begin{array}[]{l}{\nabla_{{\bf x}^{\ldots}}^{\otimes q}K\left({\bf x},{\bf x}^{\prime}\right)=\left(-1\right)^{\alpha}\mathop{\sum}\limits_{i=0}^{\left\lfloor\frac{q}{2}\right\rfloor}{\bf a}_{(i,q)}\left({\bf x}^{\prime}-{\bf x}\right)\kappa^{(q-i)}\left(\frac{1}{2}\left\|{\bf x}-{\bf x}^{\prime}\right\|_{2}^{2}\right)}\end{array}$$

$$\begin{array}[]{rl}{\bf a}_{(i,q)}\left({\bf d}\right)&\!\!\!=\mathop{\sum}\limits_{{\bf j}\in{\mathbb{J}}_{(i,q)}}\;\;\mathop{\otimes}\limits_{k=0}^{q-1}\left\{\begin{array}[]{ll}{{\bf d}}&{{\tt{if}}\;j_{k}=0}\\ {{\mbox{\boldmath$\delta$}}_{j_{k}}}&{\tt otherwise}\\ \end{array}\right.\end{array}$$

$$\begin{array}[]{l}\scriptstyle{{\mathbb{J}}_{(i,q)}=\{\left.{\bf j}\in{\mathbb{Z}}^{q}\right|\left\{j_{0},j_{1},\ldots,j_{q-1}\right\}=\left\{0,\ldots,0,-1,-1,-2,-2,\ldots,-i,-i\right\}\wedge\ldots}\\ \scriptstyle{\mathop{\operatorname{argmin}}\left\{\left.j_{k}\right|j_{k}=-1\right\}\leq\mathop{\operatorname{argmin}}\left\{\left.j_{k}\right|j_{k}=-2\right\}\leq\ldots\leq\mathop{\operatorname{argmin}}\left\{\left.j_{k}\right|j_{k}=-i\right\}\}}\end{array}$$

$$\begin{array}[]{l}{(\overbrace{{\overset{{a\;{\tt{terms}}}}{\ldots}\otimes{\mbox{\boldmath$\delta$}}_{l}\otimes\ldots}}^{{b\;{\tt{terms}}}}\otimes{\mbox{\boldmath$\delta$}}_{l}\otimes\ldots)_{\bf i}=(\delta_{i_{a},i_{b}}(\overbrace{{\overset{{a\;{\tt{terms}}}}{\ldots}\otimes{\bf 1}\otimes\ldots}}^{{b\;{\tt{terms}}}}\otimes{\bf 1}\otimes\ldots))_{\bf i}}\end{array}$$

$$\begin{array}[]{rl}{\bf x}^{*}=\mathop{\operatorname{argmax}}\limits_{{\bf x}\in{\mathbb{S}}}f\left({\bf x}\right)\end{array}$$

$$\begin{array}[]{l}L^{\downarrow q}\kappa\left(r\right)\leq\left|\kappa^{(q)}\left(r\right)\right|\leq L^{\uparrow q}\kappa\left(r\right)\;\forall q\in{\mathbb{Z}}_{s+1}\\ \left|\kappa\left(r+\delta r\right)-\mathop{\sum}\limits_{q\in{\mathbb{Z}}_{s+1}}\frac{1}{q!}\delta r^{q}\kappa^{(q)}\left(r\right)\right|\leq\Delta_{r}\left(\delta r\right)\end{array}$$

$$\begin{array}[]{l}\Delta\left(\delta r\right)=\mathop{\sup}\limits_{r\in\left[0,\frac{1}{2}M^{2}\right)}\frac{\Delta_{r}\left(\delta r\right)}{\kappa\left(r\right)}\end{array}$$

$$\begin{array}[]{r}{\Delta_{r}\left(\delta r\right)\approx\max\{\left.E_{r}\left(\delta r\right),E_{r}\left(s_{i}\right)\right|s_{0},s_{1},\ldots s_{R_{A}-1}\sim{\mathcal{U}}\left(0,\delta r\right)\}}\end{array}$$

$$\begin{array}[]{l}E_{r}\left(\delta r\right)=\left|\kappa\left(r+\delta r\right)-\mathop{\sum}\limits_{q\in{\mathbb{Z}}_{s+1}}\frac{1}{q!}\delta r^{q}\kappa^{(q)}\left(r\right)\right|\end{array}$$

$$\begin{array}[]{r}{\Delta\left(\delta r\right)\approx\max\{\left.\frac{\Delta_{r}\left(\delta r\right)}{\kappa\left(r\right)}\right|r_{0},r_{1},\ldots,r_{R_{B}-1}\sim{\mathcal{U}}\left(0,\frac{1}{2}M^{2}\right)\}}\end{array}$$

$$\begin{array}[]{l}\frac{B^{q}}{q!}\left\|\nabla_{\bf x}^{\otimes q}f\left({\bf x}\right)\right\|_{2}\leq{\mu}\end{array}$$

$$\begin{array}[]{l}{\mathbb{P}}=\left\{p\in{\mathbb{Z}}_{s}+1\left|U_{p}\left(B\right)\leq A\right.\right\}\\ \end{array}$$

$$\begin{array}[]{l}{\mathbb{S}}_{{{\mu}^{-}_{1:p}}}\subseteq{\mathbb{S}}_{{(A,B)}}\subseteq{\mathbb{S}}_{{{\mu}^{+}_{1:p}}}\end{array}$$

$$\begin{array}[]{l}F=\kappa\left(0\right)\sqrt{\frac{1}{\Gamma\left(\frac{n}{2}+1\right)}\left(\frac{\sqrt{\pi}M}{2}\right)^{n}}G\end{array}$$

$$\begin{array}[]{r}{\left|R_{q:{\bf x}}\left(\delta{\bf x}\right)\right|\leq\frac{D}{\sqrt{(q+1)!}}\frac{\left(\sqrt{2L^{\uparrow}}B\right)^{q+1}\!\!\!-\left(\sqrt{2L^{\uparrow}}B\right)^{s+1}}{1-\sqrt{2L^{\uparrow}}B}F}\\ \end{array}$$

$$\begin{array}[]{l}D=0.816\pi^{\frac{1}{4}}e^{\frac{1}{2}\left(\sqrt{L^{\uparrow}}M\right)^{2}}+\ldots\\ \frac{L^{\uparrow}-L^{\downarrow}}{L^{\uparrow}}\mathop{\sum}\limits_{i=0}^{s}\mathop{\sum}\limits_{i=0}^{\left\lfloor\frac{1}{2}\left\lfloor\frac{q}{2}\right\rfloor\right\rfloor}\frac{\sqrt{q!}}{2^{2i}(2i)!(q-4i)!}\left(\sqrt{L^{\uparrow}-L^{\downarrow}}M\right)^{q-4i}\\ \end{array}$$

$$\begin{array}[]{r}{p_{\rm min}=\max\Big{\{}1,\Big{\lceil}\left(\sqrt{2L^{\uparrow}}B\right)^{2}\exp\Big{(}1+W_{0}\Big{(}\frac{2}{e\left(\sqrt{2L^{\uparrow}}B\right)^{2}}\log\ldots}\\ {\ldots\left(\frac{1}{\sqrt{\sqrt{2\pi}}}\frac{DF}{A}\frac{1}{1-\sqrt{2L^{\uparrow}}B}\right)\Big{)}\Big{)}-1\Big{\rceil}\Big{\}}}\\ \end{array}$$

$$\begin{array}[]{l}\!\!\!\!\!\!\!\!\!\begin{array}[]{r}{p^{\rm rec}=\max\max\Big{\{}1,\Big{\lceil}\left(\sqrt{2L^{\uparrow}}B\right)^{2}\exp\Big{(}1+W_{0}\Big{(}\frac{2}{e\left(\sqrt{2L^{\uparrow}}B\right)^{2}}\log\ldots}\\ {\ldots\left(\frac{1}{\sqrt{\sqrt{2\pi}}}\frac{DF}{A}\frac{1}{1-\sqrt{2L^{\uparrow}}B}\right)\Big{)}\Big{)}-1\Big{\rceil}\Big{\}}}\\ \end{array}\!\!\!\!\!\!\!\!\!\\ \!\!\!\!\!\!\!\!\!\begin{array}[]{r}{{{{\mu}^{\pm}}}=A\pm\frac{D}{\sqrt{(q+1)!}}\frac{\left(\sqrt{2L^{\uparrow}}B\right)^{q+1}\!\!\!-\left(\sqrt{2L^{\uparrow}}B\right)^{s+1}}{1-\sqrt{2L^{\uparrow}}B}F}\\ \end{array}\!\!\!\!\!\!\!\!\!\end{array}$$

$$\begin{array}[]{rl}\left.f\left({\bf x}\right)\right|{\mathbb{D}}&\!\!\!\sim{\mathcal{N}}\left(m_{\mathbb{D}}\left({\bf x}\right),\lambda_{{\mathbb{D}}}\left({\bf x},{\bf x}^{\prime}\right)\right)\\ \left.{\nabla}_{\bf x}^{\otimes p}f\left({\bf x}\right)\right|{\mathbb{D}}&\!\!\!\sim{\mathcal{N}}\left({\bf m}_{\mathbb{D}}^{(p)}\left({\bf x}\right),{{\mbox{\boldmath$\Lambda$}}}_{{\mathbb{D}}}^{(p)}\left({\bf x},{\bf x}^{\prime}\right)\right)\\ \end{array}$$

$$\begin{array}[]{rll}s_{{{\mu}_{q}}}\left(\left.{\bf x}\right|{\mathbb{D}}\right)&\!\!\!\triangleq\Pr\left(\left.{\bf x}\in{\mathbb{S}}_{{{\mu}_{q}}}\right|{\mathbb{D}}\right)&\!\!\!=\Pr\left(\left\|{\bf v}_{(q)}\right\|_{2}\leq{\mu}\right)\\ \end{array}$$

$$\begin{array}[]{rll}s_{{{\mu}_{1:p}}}\left(\left.{\bf x}\right|{\mathbb{D}}\right)&\!\!\!\triangleq\Pr\left(\left.{\bf x}\in{\mathbb{S}}_{{{\mu}_{1:p}}}\right|{\mathbb{D}}\right)&\!\!\!=\!\!\mathop{\prod}\limits_{q\in{\mathbb{Z}}_{p}+1}\!s_{{{\mu}_{q}}}\left(\left.{\bf x}\right|{\mathbb{D}}\right)\\ \end{array}$$

$$\begin{array}[]{l}{\bf v}\sim{\mathcal{N}}\left(\frac{B^{q}}{q!}{\bf m}_{{\mathbb{D}}}^{(q)}\left({\bf x}\right),\left(\frac{B^{q}}{q!}\right)^{2}{\mbox{\boldmath$\Lambda$}}_{{\mathbb{D}}}^{(q)}\left({\bf x},{\bf x}\right)\right)\end{array}$$

$$\begin{array}[]{l}\left\|\left.f\right|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}\right\|_{W_{2}^{p,\infty}}\triangleq\mathop{\sup}\limits_{q\in{\mathbb{Z}}_{p}+1}\left\|D^{(q)}f\right\|_{L_{2}^{\infty}({\mathbb{S}}_{{{\mu}_{1:p}}})}\leq{\mu}\end{array}$$

$$\begin{array}[]{rl}g\left({\mathbb{F}}\right)&\!\!\!=y^{+}_{{\mathbb{F}}}-\chi\end{array}$$

$$\begin{array}[]{l}g\left.\left(\left\{\left({\bf x},f\left({\bf x}\right)\right)\right\}\right)\right|{\mathbb{F}}\sim{\mathcal{N}}_{\left(\chi,\infty\right)}\left(m_{g;{\mathbb{F}}}\left({\bf x}\right),\lambda_{g;{\mathbb{F}}}\left({\bf x},{\bf x}\right)\right)\end{array}$$

$$\begin{array}[]{rl}m_{g:{\mathbb{F}}}\left({\bf x}\right)&\!\!\!=m_{{\mathbb{F}}}\left({\bf x}\right)-\chi+\frac{\phi\left(\chi\right)}{\Phi\left(\chi\right)-1}\lambda_{{\mathbb{F}}}\left({\bf x},{\bf x}\right)\\ \lambda_{g:{\mathbb{F}}}\left({\bf x},{\bf x}\right)&\!\!\!=\lambda_{{\mathbb{F}}}\left({\bf x},{\bf x}\right)\left(1+\frac{\chi\phi\left(\chi\right)}{\Phi\left(\chi\right)-1}-\left(\frac{\phi\left(\chi\right)}{\Phi\left(\chi\right)-1}\right)^{2}\right)\\ \end{array}$$

Full text

1 Introduction

A canonical application of Bayesian optimisation is experimental design. Typically one aims to find the optimal experimental parameters - ratios of chemicals, temperatures etc - that maximise some form of experimental yield or return. Implicit in this task is the assumption of repeatability, specifically that if we run the same experiment twice we will obtain the same result. However in all physical experiments there are limitations (both practical and financial) on how precisely one can control the experimental conditions such as ingredient quality (eg type and quantity of any impurities) or oven temperature, and this intrinsic imprecision will manifest in variability in experimental outcomes. If this variability is small then it may be acceptable, but when it is significant it may represent the difference between a good outcome (for example an alloy that is strong and lightweight for aircraft design) or an unacceptable one.

Similar problems also arise outside of the industrial and experimental setting. [11, 10] observes that when tuning hyperparameters we may see the phenomena of false maxima, which are sharp peaks in the performance surface that may be present when the testing set is small that disappear altogether when the size of the testing set increases. Subsequently a simple Bayesian otimisation for hyper-parameter selection may recommend “optimal” hyper-parameters that refer to “optima” that have no objective reality, being a figment of the (small) training set.

Our aim in this paper is twofold. First we show how (in)stability may be characterised and detected using Gaussian Process models, and secondly we show how Bayesian optimisation may be steered to avoid unstable regions and only report stable optima. We begin by characterising instability in terms of maximal output variation bounds given specified (bound) input perturbations: we call this ${(A,B)}$-stability. We then demonstrate how gradient bounds on the first $p$ derivatives (which we call ${{\mu}_{1:p}}$-stability) may be used as a surrogate for ${(A,B)}$-stability, and how the probability of a function being ${{\mu}_{1:p}}$-stable at a point may be calculated using gradient Gaussian process models. Finally we present two modified acquisition function that may be used in Bayesian optimisation to steer the procedure away from unstable regions and toward stable ones.

1.1 Notation

Sets are written ${\mathbb{A}},{\mathbb{B}},\ldots$; where ${\mathbb{R}}_{+}$ is the positive reals, ${\mathbb{Z}}_{+}=\{1,2,\ldots\}$, ${\mathbb{Z}}_{n}=\{0,1,\ldots,n-1\}$, and $\bar{{\mathbb{R}}}_{+}$ is the non-negative reals. $|{\mathbb{A}}|$ is the cardinality of ${\mathbb{A}}$. Column vectors are bold lower case ${\bf a},{\bf b},\ldots$. Matrices are bold upper case ${\bf A},{\bf B},\ldots$. Element $i$ of vector ${\bf a}$ is $a_{i}$. Element $i,j$ of matrix ${\bf W}$ is $W_{i,j}$. ${\bf a}^{{\rm T}}$ is the transpose, ${\bf a}\otimes{\bf b}$ the Kronecker product, and ${\bf a}^{\otimes p}={\bf a}\otimes\overset{p\;{\rm terms}}{\ldots}\otimes{\bf a}$ the Kronecker power. ${\bf 1}$ a vector of $1$s, ${\bf 0}$ a vector of [math]s, and ${\bf I}$ the identity matrix. ${\nabla}_{\bf x}=[\frac{\partial}{\partial x_{0}}\,\frac{\partial}{\partial x_{1}}\,\ldots\,\frac{\partial}{\partial x_{n-1}}]^{{\rm T}}$. The indicator function is denoted $\mathbbold{1}\left({\tt Q}\right)$ and is $1$ if boolean ${\tt Q}$ is true, [math] otherwise. Logical conjunction is indicated with $\wedge$. Logical disjunction is indicated with $\vee$. The principle branch of the Lambert $W$-function is denoted $W_{0}$. The PDF and CDF of the standard normal distribution are denoted $\phi$ and $\Phi$, respectively.

2 Background

Bayesian optimisation [3, 7, 18, 6] is an optimisation technique designed for optimising expensive (in terms of economic cost, time etc) functions $f$ in the fewest evaluations possible. A Bayesian optimiser maintains a model of $f$ (usually a Gaussian process, as described shortly). At each iteration $t$ the optimiser selects a sample ${\bf x}_{t}\in{\mathbb{X}}$ to maximise an acquisition function $a_{t}:{\mathbb{X}}\to{\mathbb{R}}$ based on this model. This point is evaluated (often noisily) to obtain $y_{t}=f({\bf x}_{t})+\epsilon_{i}$, the model updated, and the process repeated. Acquisition functions are designed to trade-off exploitation of known-good regions and exploration of unknown ones. Typical acquisition functions include expected improvement (EI) [7], GP-UCB [18] and Predictive Entropy Search (PES) [6].

2.1 Gaussian Processes and Derivatives

A gaussian process $\operatorname{{\mathcal{GP}}}(\mu,K)$ is a distribution on a space of functions with mean $\mu:{\mathbb{R}}^{n}\to{\mathbb{R}}$ and covariance $K:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{R}}$. Assume $f:{\mathbb{X}}\subseteq{\mathbb{R}}^{n}\to{\mathbb{R}}\sim\operatorname{{\mathcal{GP}}}(0,K({\bf x},{\bf x}^{\prime}))$ is a draw from an unbiased Gaussian process [8, 15]. The posterior of $f$ given ${\mathbb{D}}=\{({\bf x}_{i},y_{i})|y_{i}=f({\bf x}_{i})+\epsilon_{i},\epsilon_{i}\sim{\mathcal{N}}(0,\sigma^{2})\}$ is $f({\bf x})|{\mathbb{D}}\sim{\mathcal{N}}(m_{\mathbb{D}}({\bf x}),\lambda_{\mathbb{D}}({\bf x},{\bf x}))$, where:

[TABLE]

${\bf y},{\bf k}({\bf x})\in{\mathbb{R}}^{|{\mathbb{D}}|}$, ${\bf K}\in{\mathbb{R}}^{|{\mathbb{D}}|\times|{\mathbb{D}}|}$, $k({\bf x})_{i}=K({\bf x},{\bf x}_{i})$, and $K_{i,j}=K({\bf x}_{i},{\bf x}_{j})$.

The gradient of a Gaussian process is an (independent [19]) Gaussian process [13, 14, 17] if the kernel is differentiable, and so on too are higher order gradients of Gaussian processes. In vectorised form, denoting the Kronecker power ${\bf a}^{\otimes q}={\bf a}\otimes\overset{q\;{\rm terms}}{\ldots}\otimes{\bf a}$, the posterior of ${\nabla}_{\bf x}^{\otimes q}f$ given ${\mathbb{D}}$ is {\nabla}_{\bf x}^{\otimes q}f({\bf x})|{\mathbb{D}}\sim{\mathcal{N}}({\bf m}_{\mathbb{D}}^{(q)}({\bf x}),{{\mbox{\boldmath\Lambda}}}_{{\mathbb{D}}}^{(q)}({\bf x},{\bf x})), where:

[TABLE]

and we note that:

[TABLE]

Relevant gradient calculations for standard $K$ functions can be found in [9]. Alternatively for the isotropic kernels:

[TABLE]

assuming $\kappa$ is differentiable in closed form the following result, along with table 1, may be used to calculate the required derivatives:

Theorem 1

Let $K({\bf x},{\bf x}^{\prime})=\kappa(\frac{1}{2}\|{\bf x}-{\bf x}^{\prime}\|_{2}^{2})$ be an isotropic kernel, where $\kappa$ is $s$-times differentiable. Denote by $\nabla_{{\bf x}^{\ldots}}^{\otimes q}$ a mixed Kronecker derivative of order $q$ (e.g. $\nabla_{{\bf x}^{\ldots}}^{\otimes 2}$ may be $\nabla_{{\bf x}}\otimes\nabla_{{\bf x}}$, $\nabla_{{\bf x}^{\prime}}\otimes\nabla_{{\bf x}^{\prime}}$, $\nabla_{{\bf x}}\otimes\nabla_{{\bf x}^{\prime}}$ or $\nabla_{{\bf x}^{\prime}}\otimes\nabla_{{\bf x}}$), where $\alpha$ is the number of times $\nabla_{{\bf x}^{\prime}}$ appears in $\nabla_{{\bf x}^{\ldots}}^{\otimes q}$. Then $\forall q\in{\mathbb{Z}}_{s+1}$:

[TABLE]

where $\kappa^{(c)}(x)=\frac{\partial^{c}}{\partial x^{c}}\kappa(x)$;

[TABLE]

[TABLE]

and we have used the symbolic notation (where ${\bf i}\in{\mathbb{Z}}_{n}^{q}$ is a multi-index, noting that {\mbox{\boldmath\delta}}_{l}’s appear in pairs in ${\bf a}_{(i,q)}$ $\forall l=-1,-2,\ldots,-i$):

[TABLE]

Proof:

The complete proof of this theorem is presented in the appendix. The proof begins by assuming that $\alpha=0$ (that is, $\nabla_{{\bf x}^{\prime}}$ does not appear in the Kronecker gradient, so $\nabla_{{\bf x}^{\ldots}}^{\otimes q}K({\bf x},{\bf x}^{\prime})=\nabla_{{\bf x}}^{\otimes q}K({\bf x},{\bf x}^{\prime})$) and proving the special case inductively. The general case $\alpha\geq 0$ follows by observing the sign anti-symmetry of the gradients with respect to ${\bf x}$ and ${\bf x}^{\prime}$. $\square$

3 Problem Statement

Let $f:{\mathbb{X}}\to{\mathbb{R}}_{+}$. We assume that $f$ may be evaluated (with noise and significant expense) but that its derivatives may not. Our aim is to find the stable maxima:

[TABLE]

where ${\mathbb{S}}\subseteq{\mathbb{X}}$ is the stable subset of ${\mathbb{X}}$. To achieve this we must (a) quantify what we mean by stability in practical terms, and (b) incorporate this into the acquisition function used by the Bayesian optimiser.

3.1 Assumptions

For the purposes of this paper we assume:

${\mathbb{X}}\subset{\mathbb{R}}^{n}$ compact, $\|{\bf x}-{\bf x}^{\prime}\|_{2}\leq M$ $\forall{\bf x},{\bf x}^{\prime}\in{\mathbb{X}}$. 2. 2.

$f:{\mathbb{X}}\subseteq{\mathbb{R}}^{n}\to{\mathbb{R}}_{+}\sim\operatorname{{\mathcal{GP}}}(0,K({\bf x},{\bf x}^{\prime}))$. 3. 3.

$\|f\|_{{\mathbb{H}}_{K}}\leq G$, where $\|\cdot\|_{{\mathbb{H}}_{K}}$ is the reproducing kernel Hilbert space norm. 4. 4.

$K({\bf x},{\bf x}^{\prime})=\kappa(\frac{1}{2}\|{\bf x}-{\bf x}^{\prime}\|_{2}^{2})$ is isotropic kernel (covariance), $\kappa$ is completely monotone, positive, $s$-times differentiable, and there exist $L^{\uparrow}\geq L^{\downarrow}\in{\mathbb{R}}_{+}$, $\Delta_{r}:{\mathbb{R}}_{+}\to{\mathbb{R}}_{+}$ non-decreasing such that:

[TABLE]

and we define the overall Taylor bound for $\kappa$ as:

[TABLE]

Of these assumptions only assumption 4 is the only non-trivial. We have considered only isotropic kernels as these represent the most common kernels in the Gaussian process literature, and restricted our choice to positive (valued) kernels (excluding for example the wave kernel) rather than Bernstein to enable us to construct various bounds on the remainder of $f$. The parameters $L^{\uparrow},L^{\downarrow}$ (and their existance and finiteness) is required to allow us to bound the Taylor expansion of $f$, which forms the basis of our defintion of stability, while the non-decreasing (in $\delta r$) bound on the remainder of the Taylor expansion is a convenience factor allowing us to use a richer range of (non-infinitely-differentiable) kernels. Examples of kernels satisfying the conditions of this assumption are presented in table 1.

On a technical point, we note that the remainder bounds $\Delta_{r}$, $\Delta$ can be difficult to calculate in closed form. As discussed in the appendix, if a (tight) closed-form bound is not available then these terms may be approximated using Monte-Carlo simulation [4]. Specifically:

[TABLE]

where:

[TABLE]

is a tight bound on the absolute remainder of the Taylor expansion of $\kappa$, and samples are drawn to ensure $\Delta_{r}(\delta r)$ is increasing with respect to $\delta r$. Obviously more samples $R_{A}$ will give a more accurate bound, while fewer samples will be faster to evaluate. Likewise:

[TABLE]

where the total number of samples required for this approximation is $R_{A}R_{B}$. We note that this need only be calculated twice in our algorithm, so it is feasible to use a larger number of samples to ensure accuracy. See appendix for further discussion and relevant derivations.

3.2 Related Work

The works most closely related to the present work are unscented Bayesian optimisation [12] and stable Bayesian optimisation [11, 10]. Both of these works attempt to find stability in terms of input noise by translating it to output (target) noise. [12] does this using the unscented transformation, while [11, 10] constructs a new acquisition function combining the effects of epistemic variance (“standard” variance in the output due to limited samples and noisy measurements) and aleatoric variance due to input perturbations translated into output through the objective function. Thus unstable regions of the objective function become regions of high uncertainty, which the algorithm may subsequently avoid. However there is no guarantee that such approaches will avoid unstable regions, particularly those that combine instability and particularly high (relative) return, so variability of results may still be a problem.

4 Stability - Definition and Quantification

In this section we present two definitions of stability, ${(A,B)}$-stability and ${{\mu}_{1:p}}$-stability. ${(A,B)}$-stability is defined in terms of the sensitivity of the output to variation in the input - the smaller $|f({\bf x})-f({\bf x}+\delta{\bf x})|$ is for bounded $\delta{\bf x}$, the more stable $f$ is at ${\bf x}\in{\mathbb{X}}$. This is a practical definition for the experimenter, but is difficult to quantify in practice. Alternatively, ${{\mu}_{1:p}}$-stability defines stability in terms of gradients (to order $p$). This is not as useful for the experimenter, but, as we will show, may be readily quantified using gradient Gaussian processes. In this section we will relate these two definitions and demonstrate that ${{\mu}_{1:p}}$-stability may be used as a surrogate for ${(A,B)}$-stability, allowing the experimenter to specify stability constraints in the more practical ${(A,B)}$-stable form and then enforce them in terms of the more practical ${{\mu}_{1:p}}$-stability form.

4.1 Defining Stability

${(A,B)}$-stability is defined as follows:

Definition 1** (${(A,B)}$-stability)**

Let $A,B\in{\mathbb{R}}_{+}$. We say that $f$ is ${(A,B)}$-stable at point ${\bf x}\in{\mathbb{X}}$ if $|f({\bf x}+\delta{\bf x})-f({\bf x})|\leq A$ $\forall\delta{\bf x}:\|\delta{\bf x}\|_{2}\leq B$. The set of all ${(A,B)}$-stable points for $f$ is denoted ${\mathbb{S}}_{{(A,B)}}$.

Intuitively a function $f$ is ${(A,B)}$-stable at ${\bf x}$ if input perturbation of magnitude less than $B$ leads to output variation of magnitude less than $A$.

Alternatively, stability may be defined by bounding the derivatives of $f$ up to some order $p$. This is motivated by the observation that, if the derivative $\nabla_{\bf x}f({\bf x})$ is large then small changes in ${\bf x}$ will lead to large changes in $f({\bf x})$; and if the vectorised Hessian $\nabla_{\bf x}^{\otimes 2}f({\bf x})$ is large then, even if the gradient is small at ${\bf x}$, small (finite) changes in ${\bf x}$ may nevertheless cause us to “fall off” the sharp (unstable) peak at this point. Thus we would like to label regions with large derivatives $\nabla_{\bf x}f({\bf x})$ or large vectorised Hessian $\nabla_{\bf x}^{\otimes 2}f({\bf x})$ as unstable; hence, generalising to arbitrary order, we define ${{\mu}_{1:p}}$-stability by:

Definition 2** (${{\mu}_{1:p}}$-stability)**

Let ${\mu},B\in{\mathbb{R}}_{+}$ and $p\in{\mathbb{Z}}_{+}$. We say that $f$ is ${{\mu}_{1:p}}$-stable at point ${\bf x}$ if $\forall q\in{\mathbb{Z}}_{p}+1$:

[TABLE]

The set of all ${{\mu}_{1:p}}$-stable points for $f$ is denoted ${\mathbb{S}}_{{{\mu}_{1:p}}}$. We also say that $f$ is ${{\mu}_{q}}$-stable at ${\bf x}$ for a given $q\in{\mathbb{Z}}_{+}$ if the gradient bound is met for the $q$ specified.

4.2 Connection Between ${(A,B)}$- and ${{\mu}_{1:p}}$-Stability

The forms of stability we have defined (${(A,B)}$-stability and ${{\mu}_{1:p}}$-stability) are related through the following key result, which (a) shows that ${(A,B)}$-stability is equivalent to ${{\mu}_{1:p}}$-stability in the limit $p\to\infty$ for appropriately conditions on $f$ and selected $\mu$ and (b) suggests how the paremeters $p$ and ${\mu}$ may be selected given $A,B\in{\mathbb{R}}_{+}$ and the specifics of the kernel $K$ ($s$, $L_{(q)}$, $\Delta_{r}(\delta r)$ and $\Delta(\delta r)$, as per section 3.1).

Theorem 2

Let $A,B\in{\mathbb{R}}_{+}$, $s\in{\mathbb{Z}}_{+}$. Under the default assumptions, suppose the remainders of $f({\bf x}+\delta{\bf x})$ Taylor expanded about ${\bf x}\in{\mathbb{X}}$ to order $q$ satisfy the bound $|R_{q:{\bf x}}(\delta{\bf x})|\leq U_{q}(B)$ $\forall\delta{\bf x},\|\delta{\bf x}\|_{2}\leq B$. Define:

[TABLE]

If ${\mathbb{P}}\neq\emptyset$, $p\in{\mathbb{P}}$, and ${{{\mu}^{\pm}}}=(A\pm U_{p}(B))$ then, using ${{\mu}^{+}_{1:p}}$-stability and ${{\mu}^{-}_{1:p}}$-stability to denote ${{\mu}_{1:p}}$-stability with, respectively, ${\mu}={{\mu}^{+}}$ and ${\mu}={{\mu}^{-}}$, we have:

[TABLE]

Proof:

This follows from the definitions in a straightforward manner applying standard inequalitites. See appendix for details. $\square$

This theorem suggests that we may use ${{\mu}_{1:p}}$-stability as a proxy for ${(A,B)}$-stability, and suggests a range ${\mu}\in[{{\mu}^{-}},{{\mu}^{+}}]$ of choices for ${\mu}$ to approximate ${(A,B)}$-stability given $A,B\in{\mathbb{R}}_{+}$, as shown for example in figure 1. This is desirable because the derivatives of a Gaussian Process are Gaussian Processes (to order $s$, see section 2.1), which will allow us to directly calculate the probability that $f$ is ${{\mu}_{1:p}}$-stable at a point ${\bf x}$ given observations ${\mathbb{D}}$, which allows us to quantify the expected gain for a particular recommendation and thus construct a sensible acquisition function for our Bayesian optimiser. Note that:

• •

Smaller ${\mu}$ (e.g. ${\mu}={{\mu}^{-}}$) defines a conservative approximation excluding marginally stable points, while larger ${\mu}$ (e.g. ${\mu}={{\mu}^{+}}$) defines a more liberal approximation possibly including marginally unstable points.

• •

If the Taylor expansion of $\kappa$ converges (so $U_{q}(B)$ decreases with $q$) then larger $p$ values will result in better approximation of ${(A,B)}$-stability. However this must be balanced against the computational cost of calculating means and variances of $n^{p}$-dimensional ($p^{\rm th}$-order derivative) Gaussian processes. In practice we found this to be of little concern as $p\leq 2$ typically suffices, which bounds the gradient and Hessian, where bounding the gradient excludes unstable maxima on the boundaries of ${\mathbb{X}}$,111Other maxima will have zero gradient by first-order optimality conditions. and bounding the Hessian excludes unstable, quadratic-type maxima.

The convergence rate of the Taylor expansion of $f$ depends on the isotropic kernel function $K$ of the Gaussian process from which $f$ was drawn, as quantified by the following theorem (proven in the appendix), where for clarity we consider the simplified case $s=\infty$, $\Delta(r)=0$ (the more general case is presented in the appendix):

Theorem 3

Under the default assumptions $|f({\bf x})|\leq F$ $\forall{\bf x}\in{\mathbb{X}}$, where:

[TABLE]

and the remainders of $f({\bf x}+\delta{\bf x})$ Taylor expanded around ${\bf x}\in{\mathbb{X}}$ to order $q\in{\mathbb{Z}}_{s}+1$ are bounded by:

[TABLE]

$\forall\delta{\bf x}:\|\delta{\bf x}\|_{2}\leq B<\frac{1}{\sqrt{2L^{\uparrow}}}$, where:

[TABLE]

Moreover $|R_{p:{\bf x}}(\delta{\bf x})|\leq A$ $\forall p\geq p_{\rm min}$, where:

[TABLE]

where $W_{0}$ is the principle branch of the Lambert $W$-function.

Proof:

A proof is given in the appendix. Several steps are required. As preliminary, we show that the number of terms in the gradients $\nabla_{\bf x}^{\otimes q}K({\bf x},{\bf x}^{\prime})$ is equal to the number of terms in the Hermite polynomial $H_{q}$. This is leveraged to construct a bound on the remainders of the Taylor expansion of $\kappa$. Noting that $f$ is in a reproducing kernel Hilbert space, the bound on the remainder of $\kappa$ is used to bound the remainder of the Taylor expansion of $f$. Finally Stirlings approximation is used to find $p_{\rm min}$. $\square$

This theorem provides the details required to use ${{\mu}_{1:p}}$-stability as a proxy for ${(A,B)}$-stability, as suggested by theorem 2. In particular, it suggests that we choose $p=p^{\rm rec}$ and ${\mu}\in[{{\mu}^{-}},{{\mu}^{+}}]$, where, using the constants in the theorem, and provided $B<\frac{1}{\sqrt{2L}}$:

[TABLE]

where $F=\kappa(0)\sqrt{\frac{1}{\Gamma(\frac{n}{2}+1)}(\frac{\sqrt{\pi}M}{2})^{n}}G$. Note that the restriction $B<{1}/{\sqrt{2L^{\uparrow}}}$ on input variation is actually a requirement that the input variation be less than an amount proportional to the (effective) length-scale of the kernel $K$.

Finally we note that in practice we have observed that it is almost never necessary to test ${{\mu}_{1:p}}$-stability past $3^{\rm rd}$-order (or even $2^{\rm nd}$ order) in most cases when using an RBF kernel. This appears to be due to two factors:

• •

The scaled gradients $\frac{B^{q}}{q!}\nabla_{\bf x}^{\otimes q}f({\bf x})$ taper off much more quickly than the bounds in theorem 3 may suggest, presumably due to the number of approximations (upper bounds) required to obtain the said bounds. For example in figure 1 we see that by $3^{\rm rd}$-order the scaled gradients fall well within the bounds of ${{\mu}_{1:p}}$-stability.

• •

Even if a higher-order derivative fails to meet the bound requirement $\|\frac{B^{q}}{q!}\nabla_{\bf x}^{\otimes q}f({\bf x})\|_{2}\leq{\mu}$, usually a lower-order derivative will also fail to meet this bound, rendering the (more computationally expensive) higher-order test superfluous.

Next we consider how the stability of a point may be quantified when the derivatives are approximated using the derivatives of the Gaussian process model of $f$.

4.3 Quantifying Stability

We now show how the derivatives of the Gaussian process model of $f$ may be used to calculate the posterior probability that $f$ is ${{\mu}_{1:p}}$-stable at ${\bf x}\in{\mathbb{X}}$. Using the notation of section 2.1, given ${\mathbb{D}}=\{({\bf x}_{i},y_{i})|y_{i}=f({\bf x}_{i})+\epsilon_{i},\epsilon_{i}\sim{\mathcal{N}}(0,\sigma^{2})\}$:

[TABLE]

where means and variances are given by (1) and (2). This allows us to calculate the posterior probabilities of ${{\mu}_{q}}$-stability and ${{\mu}_{1:p}}$-stability, specifically:

Theorem 4

The posterior probability of $f$ being ${{\mu}_{q}}$-stable at ${\bf x}$ given ${\mathbb{D}}$ is:

[TABLE]

where {\bf v}_{(q)}\sim{\mathcal{N}}\big{(}\frac{B^{q}}{q!}{\bf m}_{{\mathbb{D}}}^{(q)}({\bf x}),(\frac{B^{q}}{q!})^{2}{\mbox{\boldmath\Lambda}}_{{\mathbb{D}}}^{(q)}({\bf x},{\bf x})\big{)}, and posterior probability of $f$ being ${{\mu}_{1:p}}$-stable at ${\bf x}$ is:

[TABLE]

Proof:

The first result follows from the properties of the Gaussian process model of $f$, and the second from the fact that $\nabla_{\bf x}f$, $\nabla_{\bf x}^{\otimes 2}f$, $\ldots$ are independent. $\square$

We call $s_{{{\mu}_{1:p}}}({\bf x}|{\mathbb{D}})$ the stability score of ${\bf x}$ given ${\mathbb{D}}$. These stability scores form the basis for our proposed acquisition functions in subsequent sections. Stability scores may be calculated by Monte-Carlo estimation [4]. That is, generate a set of random vectors:

[TABLE]

and test what fraction satisfy $\|{\bf v}\|_{2}\leq{\mu}$. Note that $\|{\bf v}\|_{2}$ is $1$-dimensional, so the number of samples required to achieve a given accuracy does not depend on the dimension $n$ or the order $q$.

4.4 Connection to Sobolev Norms

As an aside, it is interesting to note the connection between ${{\mu}_{1:p}}$-stability and Sobolev norms. If we let $D^{(q)}=\frac{B^{q}}{q!}\nabla_{\bf x}^{\otimes q}$ be a (scaled) derivative operator and denote by $f|_{{\mathbb{S}}}$ the restriction of $f$ to ${\mathbb{S}}\subset{\mathbb{X}}$, we see that ${\mathbb{S}}_{{{\mu}_{1:p}}}$ is the largest subset of ${\mathbb{X}}$ such that the Sobolev-type seminorm222To make this a Sobolev norm $f$ would also need to be bounded. Without this additional requirement it may be seen that $\|(f+g)|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}\|_{W_{2}^{p,\infty}}\leq\|f|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}\|_{W_{2}^{p,\infty}}+\|g|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}\|_{W_{2}^{p,\infty}}$ and $\|af|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}\|_{W_{2}^{p,\infty}}=|a|\|f|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}\|_{W_{2}^{p,\infty}}$, but $\|f|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}\|_{W_{2}^{p,\infty}}=0$ for all non-varying $f$, so this is a seminorm rather than a norm. of $f|_{{\mathbb{S}}_{{{\mu}_{1:p}}}}$ ($f$ restricted to ${\mathbb{S}}_{{{\mu}_{1:p}}}$) satisfies:

[TABLE]

where $\|{\bf g}\|_{L_{2}^{\infty}({\mathbb{S}})}\triangleq\sup_{{\bf x}\in{\mathbb{S}}}\|{\bf g}({\bf x})\|_{2}$.

5 Stable Bayesian Optimisation

Having established preliminary results we now move on to define our stable Bayesian optimisation algorithm. We do this in two parts: first we construct stable forms of the expected improvement (EI) [3, 7] and GP upper confidence bound (GP-UCB) [18] acquisition functions, then we present the complete stable Bayesian optimisation algorithm.

5.1 Gain, Stable Gain and Acquisition Functions

We begin by introducing the concept of gain:

Definition 3** (Gain)**

Let $\chi\in{\mathbb{R}}$ be a lower bound on $f$, and let ${\mathbb{F}}=\{(\tilde{\bf x}_{i},\tilde{y}_{i})|\tilde{y}_{i}=f(\tilde{\bf x}_{i})+\tilde{\epsilon}_{i}\}$ be a set of observations of $f$. The gain of ${\mathbb{F}}$ is the maximum improvement over $\chi$ for any observation in ${\mathbb{F}}$:

[TABLE]

where $y^{+}_{{\mathbb{F}}}={\max}\{\chi,\tilde{y}_{i}|(\tilde{\bf x}_{i},\tilde{y}_{i})\in{\mathbb{F}}\}$.

Recall that the posterior $f({\bf x})|{\mathbb{F}}\sim{\mathcal{N}}(m_{\mathbb{F}}({\bf x}),\lambda_{\mathbb{F}}({\bf x},{\bf x}))$ is normally distributed under the default assumptions. It follows that:

[TABLE]

follows a truncated normal distribution with:

[TABLE]

where $\phi$ and $\Phi$ are the PDF and CDF of the standard normal distribution. Note that we may write the EI [3, 7] and GP-UCB [18] acquisition functions in terms of the gain:

[TABLE]

where $z({\bf x})=\frac{m_{\mathbb{D}}({\bf x})-y_{{\mathbb{D}}}^{+}}{\lambda_{\mathbb{D}}^{1/2}({\bf x},{\bf x})}$.

We wish to reformulate these acquisition functions so that only points at which $f$ is ${{\mu}_{1:p}}$-stable contribute to the result. Our approach is to re-write these in terms of the ${{\mu}_{1:p}}$-stable gain, which we define to be the gain due to the subset of ${{\mu}_{1:p}}$-stable points in the set of observations ${\mathbb{F}}$ - that is:

Definition 4** (Stable Gain)**

Let $\chi\in{\mathbb{R}}$ be a lower bound on $f$, and let ${\mathbb{F}}=\{(\tilde{\bf x}_{i},\tilde{y}_{i})|\tilde{y}_{i}=f(\tilde{\bf x}_{i})+\tilde{\epsilon}_{i}\}$ be a set of observations of $f$. The ${{\mu}_{1:p}}$-stable gain* of ${\mathbb{F}}$ is the maximum improvement over $\chi$ for any ${{\mu}_{1:p}}$-stable observation in ${\mathbb{F}}$:*

[TABLE]

where $y^{+}_{{\mathbb{F}}{{\mu}_{1:p}}}={\max}\{\chi,\tilde{y}_{i}|(\tilde{\bf x}_{i},\tilde{y}_{i})\in{\mathbb{F}}\wedge\tilde{\bf x}_{i}\in{\mathbb{S}}_{{{\mu}_{1:p}}}\}$.

As usual, under the default assumptions the posterior $f({\bf x})|{\mathbb{F}}\sim{\mathcal{N}}(m_{\mathbb{F}}({\bf x}),\lambda_{\mathbb{F}}({\bf x},{\bf x}))$ is normally distributed. It is readily seen that:

[TABLE]

follows a truncated normal distribution with:

[TABLE]

By analogy with the (standard) EI and GP-UCB acquisition functions we define the expected improvement in stable gain (EISG) and stable GP-UCB (UCBSG) acquisition functions:

Definition 5** (EISG Acquisition Function)**

The expected improvement in stable gain (EISG) acquisition function is:

[TABLE]

Definition 6** (UCBSG Acquisition Function)**

The GP-UCB in stable gain (UCBSG) acquisition function is:

[TABLE]

These may be calculated with the help of the theorems:

Theorem 5

Let ${\mathbb{D}}=\{({\bf x}_{i},{y}_{i})|{y}_{i}=f({\bf x}_{i})+\epsilon_{i}\}$. Assume without loss of generality that ${y}_{0}\leq{y}_{1}\leq\ldots$ and define $y_{-1}=\chi$, $y_{|{\mathbb{D}}|}=\infty$. Under the usual assumptions the EISG acquisition function reduces to:

[TABLE]

where $z_{i}({\bf x})=\frac{m_{\mathbb{D}}({\bf x})-y_{i}}{\lambda^{1/2}_{\mathbb{D}}({\bf x},{\bf x})}$, $\Delta\hat{y}_{i}({\bf x})=\frac{{y}_{i}-{y}_{i-1}}{\lambda^{1/2}_{{\mathbb{D}}}({\bf x},{\bf x})}$ and:

[TABLE]

so $\Delta\Phi_{|{\mathbb{D}}|}({\bf x})=\Phi(z_{|{\mathbb{D}}|-1}({\bf x}))$ and $\Delta\phi_{|{\mathbb{D}}|}({\bf x})=\phi(z_{|{\mathbb{D}}|-1}({\bf x}))$. The weights $\omega_{0}$, $\omega_{1}$, $\ldots$, $\omega_{|{\mathbb{D}}|}$ are given by:

[TABLE]

where $\omega_{|{\mathbb{D}}|}=1$.

Proof:

The complete proof is technical and can be found in the appendix. $\square$

Theorem 6

Let ${\mathbb{D}}=\{({\bf x}_{i},{y}_{i})|{y}_{i}=f({\bf x}_{i})+\epsilon_{i}\}$. Under the usual assumptions the EISG acquisition function reduces to:

[TABLE]

Proof:

This follows from definition 12 using (10). $\square$

Note that, in the absense of stability constraints or in the limit ${\mu}\to\infty$ the stability scores $s_{{{\mu}_{1:p}}}({\bf x}|{\mathbb{D}})\to 1$ $\forall{\bf x}\in{\mathbb{X}}$, so $\omega_{i}\to 0$ $\forall i\in{\mathbb{Z}}_{|{\mathbb{D}}|}$ and $\omega_{|{\mathbb{D}}|}=1$, sp the EISG and UCBSG acquisition functions reduce to the standard (non stability constrained) forms.

5.2 Stable Bayesian Optimisation via Direct Stability Quantification

Our Stable Bayesian optimisation via Direct Stability Quantification algorithm is presented in algorithm 1. Once the operating parameters ${\mu}$ and $p$ have been selected the algorithm proceeds as per standard Bayesian optimisation, excepting that the final recommendation is selected to maximise expected ${{\mu}_{1:p}}$-stable gain. Note that:

• •

The parameters $A,B$ control the stability constraints applied to the solution as per definition 1.

• •

The policy control parameter $\gamma\in[0,1]$ controls whether the approximation of ${(A,B)}$-stability with ${{\mu}_{1:p}}$-stability is conservative ($\gamma=0$), which may exclude some marginally stable points from the search, or liberal ($\gamma=1$), which may include marginally unstable points. Unless otherwise stated we have used a maximally conservative ($\gamma=0$) policy.

• •

The pragmatic limit parameter $p^{\rm max}$ controls the maximum order to which the stability scores are approximated. This is based on the observation that the $p$ value selected from the theory is almost always overly large, leading to excessive computational cost. Experimentally we have observed that $p^{\rm max}=3$ suffices in most cases, so this may be assumed unless otherwise stated.

• •

Based on our experimental results we recommend that the GP-UCB in stable gain acquisition function be used at all times.

6 Experimental Results

6.1 Simulated Experiments

In our first experiment we consider the simulated objective:

[TABLE]

where ${\mathbb{X}}=[0,1]$ and $\gamma=0.03535$, with stability parameters $A=0.2$, $B=0.0125$, as shown in figure 1. This function has an unstable maxima at $x=0.5$ and a stable maxima at $x=\frac{4}{5}$, as well as stable local (but not global) maxima at $x=\frac{1}{8},\frac{3}{8},\frac{1}{2},\frac{5}{8}$. It was chosen because the distinction between ${(A,B)}$-stable regions and ${(A,B)}$-unstable regions is not immediately obvious on inspection.

We have compared EISG (expected improvement in stable gain) and UCBSG (GP-UCB in stable gain) acquisition functions as well as unscented Bayesian optimisation and the stable Bayesian optimisation of [11, 10], with results shown in table 2. All experiments were repeated $10$ times. Note that neither unscented Bayesian optimisation nor [11, 10] are directly designed for this task and required some tweaking (in particular significantly increasing the variance of the input noise over that suggested by $B=0.0125$ to avoid always converging to the global maxima). Even after tweaking these algorithms still occasionally converged to the unstable maxima, so to ensure a fair comparison we have filtered out such cases.

The UCBSG acquisition function outperformed all other algorithms for this experiment. The reason for this is clear from figure 3, which shows $f({\bf x}_{t})$ and associated stable gains for recommendations over time. The EISG acquisition function tends to become “stuck” exploring the unstable global maxima, testing the same point over and over again. This provides no additional information for the gradient GP (and thus no additional information to update stability scores), as gradients are informed by the spread of samples around a point, so no additional information is gained and the process repeats. By contrast the explicit exploration term in the UCBSG acquisition function ensures a better spread of samples, so gradients (and thus stability scores) are correctly learnt and samples increasingly focus on the stable maxima.

7 Conclusions

In this paper we have studied the problem of finding stable maxima for expensive functions using a Gradient-based constraint as a surrogate for ${(A,B)}$-stability. We have also presented some theoretical analysis of the commection between ${(A,B)}$-stability and its surrogate to obtain bounds on the various parameters required. Our optimisation method is based on Bayesian optimisation. Using the novel concept of stable gain we have presented two acquisition function designed to avoid unstable regions in favour of stable solutions, namely expected improvement in stable gain (EISG) and GP upper confidence bound in stable gain (UCBSG), and experimentally we have compared these and also unscented Bayesian optimisation and stable Bayesian optimistion. Experimental results indicate that UCBSQ outperforms the alternative methods both in terms of reliability (likelihood that it will find a stable maxima) and convergence.

Appendix A Derivatives of Isotropic Kernels

Isotropic kernels [5] are kernels of the form:

[TABLE]

Assume $K$ is $s$-times differentiable. In this section we consider the calculation of derivatives up to order $s$. Assuming that:

[TABLE]

can be calculated in closed form, we have the following results:

Theorem 1

Let $K({\bf x},{\bf x}^{\prime})=\kappa(\frac{1}{2}\|{\bf x}-{\bf x}^{\prime}\|_{2}^{2})$ be an isotropic kernel, where $\kappa$ is $s$-times differentiable. Denote by $\nabla_{{\bf x}^{\ldots}}^{\otimes q}$ a mixed Kronecker derivative of order $q$ (e.g. $\nabla_{{\bf x}^{\ldots}}^{\otimes 2}$ may be $\nabla_{{\bf x}}\otimes\nabla_{{\bf x}}$, $\nabla_{{\bf x}^{\prime}}\otimes\nabla_{{\bf x}^{\prime}}$, $\nabla_{{\bf x}}\otimes\nabla_{{\bf x}^{\prime}}$ or $\nabla_{{\bf x}^{\prime}}\otimes\nabla_{{\bf x}}$), where $\alpha$ is the number of times $\nabla_{{\bf x}^{\prime}}$ appears in $\nabla_{{\bf x}^{\ldots}}^{\otimes q}$. Then $\forall q\in{\mathbb{Z}}_{s+1}$:

[TABLE]

where $\kappa^{(c)}(x)=\frac{\partial^{c}}{\partial x^{c}}\kappa(x)$;

[TABLE]

[TABLE]

and we have used the symbolic notation (where ${\bf i}\in{\mathbb{Z}}_{n}^{q}$ is a multi-index, noting that {\mbox{\boldmath\delta}}_{l}’s appear in pairs in ${\bf a}_{(i,q)}$ $\forall l=-1,-2,\ldots,-i$):

[TABLE]

Proof:

We begin by assuming $\alpha=0$. Defining ${\bf d}={\bf x}^{\prime}-{\bf x}$ we see that:

[TABLE]

which confirms the first expression in the theorem for $q\leq 4$ when $\alpha=0$. More generally, suppose that for some $q>0$:

[TABLE]

Then:

[TABLE]

and it follows by index rearrangement that:

[TABLE]

Therefore by induction $\forall q\in{\mathbb{Z}}_{+}$:

[TABLE]

The final result ($\alpha\geq 0$) follows by observing the sign-anti-symmetry of ${\bf x}$ and ${\bf x}^{\prime}$ in all expressions. $\square$

Appendix B Proof of Theorems 2 and 3

In this section we prove theorems 1 and 2 from the body of the paper. Before proceeding with this we first establish some preliminary results. Finally, we consider some examples of kernels and derive the relevant constants relating to the theorems. Throughout this section we use the shorthand:

[TABLE]

We will also be using the Hermite polynomials $H_{q}$ and the normalised Hermite polynomials (Hermite functions) $h_{q}$ [1]:

[TABLE]

where:

[TABLE]

As per the paper, it is assumed throughout that:

${\mathbb{X}}\subset{\mathbb{R}}^{n}$ compact, $\|{\bf x}-{\bf x}^{\prime}\|_{2}\leq M$ $\forall{\bf x},{\bf x}^{\prime}\in{\mathbb{X}}$. 2. 2.

$f:{\mathbb{X}}\subseteq{\mathbb{R}}^{n}\to{\mathbb{R}}_{+}\sim\operatorname{{\mathcal{GP}}}(0,K({\bf x},{\bf x}^{\prime}))$. 3. 3.

$\|f\|_{{\mathbb{H}}_{K}}\leq G$, where $\|\cdot\|_{{\mathbb{H}}_{K}}$ is the reproducing kernel Hilbert space norm. 4. 4.

$K({\bf x},{\bf x}^{\prime})=\kappa(\frac{1}{2}\|{\bf x}-{\bf x}^{\prime}\|_{2}^{2})$ is isotropic kernel (covariance), $\kappa$ is completely monotone, positive, $s$-times differentiable, and there exist $L^{\uparrow}\geq L^{\downarrow}\in{\mathbb{R}}_{+}$, $\Delta_{r}:{\mathbb{R}}_{+}\to{\mathbb{R}}_{+}$ non-decreasing such that:

[TABLE]

and we define the overall Taylor bound for $\kappa$ as:

[TABLE]

B.1 Preliminary Results

The following preliminary results are required:

Theorem I

The number of terms in the sum ${\bf a}_{(i,q)}$ as defined by (15) is:

[TABLE]

which are the same terms that occur in the Hermite polynomial (LABEL:eq:hermpoly).

Proof:

We aim to count the number of distinct vectors ${\bf j}\in{\mathbb{Z}}^{q}$ such that $\{j_{0},j_{1},\ldots,j_{q-1}\}=\{0,0,\ldots,0,-1,-1,-2,-2,\ldots,-i,-i\}$ and $\operatorname{argmin}\{j_{k}:j_{k}=-1\}\leq\operatorname{argmin}\{j_{k}:j_{k}=-2\}\leq\ldots\leq\operatorname{argmin}\{j_{k}:j_{k}=-i\}$.

Ignoring constraints, there are $q!$ permutations $(j_{0},j_{1},\ldots,j_{q-1})$ such that $\{j_{0},j_{1},\ldots,j_{q-1}\}=\{0,0,\ldots,0,-1,-1,-2,-2,\ldots,-i,-i\}$. Of these, $(q-2i)!$ are redundant reshuffles of [math] elements, $2$ are redundant reshuffles of $-1$ elements, $2$ are redundant reshuffles of $-2$ elements, $\ldots$, and $2$ are redundant reshuffles of $-i$ elements. Thus there are $\frac{q!}{2^{i}(q-2i)!}$ distinct vectors ${\bf j}$ such that $\{j_{0},j_{1},\ldots,j_{q-1}\}=\{0,0,\ldots,0,-1,-1,-2,-2,\ldots,-i,-i\}$. Note that only $1$ out of every $i!$ of these vectors satisfies the condition $\operatorname{argmin}\{j_{k}:j_{k}=-1\}\leq\operatorname{argmin}\{j_{k}:j_{k}=-2\}\leq\ldots\leq\operatorname{argmin}\{j_{k}:j_{k}=-i\}$, leaving a total of $n_{(i,q)}$ terms in the sum (each corresponding to a vector ${\bf j}$).

The final result follows from the definition of the Hermite polynomials ([1], table 22.3). $\square$

Theorem II

Under the default assumptions:

[TABLE]

$\forall q\in{\mathbb{Z}}_{s+1}$, where:

[TABLE]

where $D_{(q)}^{\updownarrow}=0$ if $L^{\uparrow}=L^{\downarrow}$.

Proof:

Recall the definition of ${\bf a}_{(i,q)}$ in theorem 1 and $n_{(i,q)}$ from theorem I. Using multi-index notation we see that:

[TABLE]

and so:

[TABLE]

By assumption 4 it follows that, defining $L^{\updownarrow}=L^{\uparrow}-L^{\downarrow}$ and letting $L_{[i]}\in[L^{\uparrow},L^{\downarrow}]$ $\forall i\in{\mathbb{Z}}_{q+1}$ such that the first statement in the following is true (this is always possible by the definition of $L^{\uparrow}$, $L^{\downarrow}$ in assumption 4 and the complete monotonicity of $\kappa$):

[TABLE]

So, by the definition of the normalised Hermite polynomial (LABEL:eq:hermpoly):

[TABLE]

Note from [1, 2] that $|h_{q}(x)|<0.816$, so:

[TABLE]

where:

[TABLE]

is non-decreasing for $t\in{\mathbb{R}}_{+}$, so by the assumed bounds:

[TABLE]

or, re-writing:

[TABLE]

Finally, noting that $\frac{L^{\uparrow}-L^{\downarrow}}{L^{\uparrow}}\in[0,1]$, it follows that:

[TABLE]

and the result follows. $\square$

Theorem III

Under the default assumptions, the remainders of $\kappa({\bf x}+\delta{\bf x})$ Taylor expanded around ${\bf x}\in{\mathbb{X}}$, to order $p\in{\mathbb{Z}}_{s+1}$ are bounded by:

[TABLE]

$\forall\delta{\bf x}:\|\delta{\bf x}\|_{2}\leq B<\frac{1}{\sqrt{2L^{\uparrow}}}$, where:

[TABLE]

where $D^{\updownarrow}=0$ if $L^{\uparrow}=L^{\downarrow}$.

Proof:

Taylor expanding to order $p\in{\mathbb{Z}}_{s+1}$:

[TABLE]

and so:

[TABLE]

Using theorem II we see that:

[TABLE]

where:

[TABLE]

is non-increasing with $p$. Moreover:

[TABLE]

which is well-defined by assumption $\sqrt{2L^{\uparrow}}B<1$. Therefore:

[TABLE]

and the result follows, noting that $\Delta_{r}$ is non-decreasing. $\square$

B.2 Main Proofs

Theorem 2

Let $A,B\in{\mathbb{R}}_{+}$, $s\in{\mathbb{Z}}_{+}$. Under the default assumptions, suppose the remainders of $f({\bf x}+\delta{\bf x})$ Taylor expanded about ${\bf x}\in{\mathbb{X}}$ to order $q$ satisfy the bound $|R_{q:{\bf x}}(\delta{\bf x})|\leq U_{q}(B)$ $\forall\delta{\bf x},\|\delta{\bf x}\|_{2}\leq B$. Define:

[TABLE]

If ${\mathbb{P}}\neq\emptyset$, $p\in{\mathbb{P}}$, and ${{{\mu}^{\pm}}}=(A\pm U_{p}(B))$ then, using ${{\mu}^{+}_{1:p}}$-stability and ${{\mu}^{-}_{1:p}}$-stability to denote ${{\mu}_{1:p}}$-stability with, respectively, ${\mu}={{\mu}^{+}}$ and ${\mu}={{\mu}^{-}}$, we have:

[TABLE]

Proof:

Suppose ${\bf x}\in{\mathbb{S}}_{{(A,B)}}$. As $f$ is $s$-times differentiable and ${(A,B)}$-stable at ${\bf x}$, and using the fact that $({\bf a}\otimes{\bf b})^{{\rm T}}({\bf c}\otimes{\bf d})=({\bf a}^{{\rm T}}{\bf c})({\bf b}^{{\rm T}}{\bf d})$:

[TABLE]

$\forall q\in{\mathbb{Z}}_{p}+1$, $\delta{\bf x}:\|\delta{\bf x}\|_{2}\leq B$. Hence:

[TABLE]

Substituting, we find that, under the conditions of the theorem $\forall q\in{\mathbb{Z}}_{p}+1$:

[TABLE]

and hence ${\bf x}\in{\mathbb{S}}_{{{\mu}^{+}_{1:p}}}$, ${\mathbb{S}}_{{(A,B)}}\subseteq{\mathbb{S}}_{{{\mu}^{+}_{1:p}}}$.

Next suppose ${\bf x}\in{\mathbb{S}}_{{{\mu}^{-}_{1:p}}}$. As $f$ is $s$-times differentiable:

[TABLE]

$\forall q\in{\mathbb{Z}}_{p}+1$. We want to show that this implies ${\bf x}\in{\mathbb{S}}_{{(A,B)}}$ or, equivalently, that $|f({\bf x})-f({\bf x}+\delta{\bf x})|\leq A$ $\forall\delta{\bf x}:\|\delta{\bf x}\|_{2}\leq B$. It suffices to show that for $p$, ${{\mu}^{-}}$ specified, $\forall\delta{\bf x}:\|\delta{\bf x}\|_{2}\leq B$:

[TABLE]

or, equivalently, that for $p$, ${{\mu}^{-}}$ specified, ${{\mu}^{-}}\leq A-U_{p}\left(B\right)>0$, which is true by definition of ${{\mu}^{-}}$ and $p$ in the theorem. $\square$

Theorem 3

Under the default assumptions $|f({\bf x})|<F$, where:

[TABLE]

and the remainders of $f({\bf x}+\delta{\bf x})$ Taylor expanded around ${\bf x}\in{\mathbb{X}}$ to order $q\in{\mathbb{Z}}_{s}+1$ are bounded by:

[TABLE]

$\forall\delta{\bf x}:\|\delta{\bf x}\|_{2}\leq B<\frac{1}{\sqrt{2L^{\uparrow}}}$, where:

[TABLE]

noting that $D^{\updownarrow}=0$ if $L^{\uparrow}=L^{\downarrow}$. If $\Delta(\frac{1}{2}B^{2})F<A$ then $|R_{p:{\bf x}}(\delta{\bf x})|\leq A$ $\forall p\geq p_{\rm min}$, where:

[TABLE]

where $W_{0}$ is the principle branch of the Lambert $W$-function.

Proof:

We have that ${\mathbb{X}}$ is compact and finite dimensional with maximum (Euclidean) distance between points in ${\mathbb{X}}$ being $M$. The maximum hypervolume of ${\mathbb{X}}$ satisfying these criteria is that of an $|{\mathbb{X}}|$-ball with diameter $M$ - that is, the Lebesgue measure of ${\mathbb{X}}$ is bounded as:

[TABLE]

We have that $f\in{\mathbb{H}}_{K}$, where ${\mathbb{H}}_{K}$ is the reproducing kernel Hilbert space associated with $K$, as $f$ is a draw from an unbiased Gaussian process with zero mean and kernel $K$. Hence $\exists\alpha\in L_{2}({\mathbb{X}})$ such that:

[TABLE]

where $\|f\|_{{\mathbb{H}}_{K}}=\|\alpha\|_{L_{2}({\mathbb{X}})}$ (as $f\in{\mathbb{H}}_{K}$ there exist at least one $\alpha\in L_{2}({\mathbb{X}})$ such that $f$ has the above form, and by definition $\|f\|_{{\mathbb{H}}_{K}}=\inf_{\alpha}\|\alpha\|_{H}$, so we choose the minimum norm $\alpha$). Using standard properties of $L_{p}$-norms, we also have that:

[TABLE]

Moreover, using Hölder’s inequality and the positivity and complete monotonicity of $\kappa$:

[TABLE]

Next, let $p\in{\mathbb{Z}}_{s+1}$. We know that:

[TABLE]

and the Taylor expansion of $f$ to order $p$ about ${\bf x}$ is:

[TABLE]

where $R_{p:{\bf x}}(\delta{\bf x})$ is the remainder; and hence:

[TABLE]

where:

[TABLE]

Defining $S_{p:\delta{\bf x}}({\bf x})=R_{p:{\bf x}}(\delta{\bf x})$, we see that:

[TABLE]

and hence, by Hölder’s inequality:

[TABLE]

Recall from theorem III that, defining $D=D^{\uparrow}+D^{\updownarrow}$:

[TABLE]

and hence:

[TABLE]

Finally we must prove the bound $p\geq p_{\rm min}$ so that $|R_{p:{\bf x}}(\delta{\bf x})|\leq A$ $\forall\delta{\bf x}$ satisfying relevant bounds. First we note that, trivially:

[TABLE]

Hence it suffices that $p_{\rm min}$ satisfies:

[TABLE]

By Stirling’s approximation we know that [16]:

[TABLE]

so it suffices to find $p_{\rm min}$ such that:

[TABLE]

or, equivalently, taking the natural log both sides and simplifying:

[TABLE]

Let $y=\frac{p_{\rm min}+1}{e(\sqrt{2L^{\uparrow}}\|\delta{\bf x}\|_{2})^{2}}$. Then the preceding equation reduces to finding $y$ such that:

[TABLE]

which is simply the inverse of the Lambert $W$-function (principle branch). Hence it suffices that:

[TABLE]

That is:

[TABLE]

which completes the proof. $\square$

Appendix C Properties of Standard Isotropic Kernel

In this section we consider the two kernels that are appropriate for our method and one counter-example to illustrate the limitations.

C.1 RBF Kernel

The RBF kernel is defined by:

[TABLE]

where $\gamma\in{\mathbb{R}}_{+}$. We see immediately that $\kappa$ is infinitely differentiable, where:

[TABLE]

It follows trivially that:

[TABLE]

where:

[TABLE]

Note also that the Taylor expansion of the RBF kernel is convergent, so:

[TABLE]

C.2 Matérn Kernels

The Matérn kernel is defined by:

[TABLE]

where $\nu,\rho\in{\mathbb{R}}_{+}$ and ${\rm H}_{\nu}$ is a modified Bessel function of the second kind.333We use ${\rm H}_{\nu}$ rather than ${\rm K}_{\nu}$ here to avoid confusion between the modified Bessel function and the kernel. From [1] ((9.6.28) with trivial rearrangement) we have that:

[TABLE]

and hence $\forall q\in{\mathbb{Z}}_{\lceil\nu\rceil}$:

[TABLE]

Note that, while this indicates that derivatives do exist to arbitrary order for $\nu\in{\mathbb{R}}_{+}\backslash{\mathbb{Z}}$ as ${\rm H}_{-\nu}={\rm H}_{\nu}$ (the gamma function has poles at $\nu=0,-1,-2,\ldots$, so the derivative is ill-defined if $\nu\in{\mathbb{Z}}_{+}$ and $q\geq\nu$), the result only defines a kernel for $q\leq{\lceil\nu\rceil-1}$. Thus the derivatives of a Gaussian process with a Matérn kernel are only Gaussian processes to order $q\leq{\lceil\nu\rceil-1}$. Of equal importance here, the derivatives of order $q>{\lceil\nu\rceil-1}$ have a pole at $r=0$, so the Taylor series approximation will construct in our proofs will diverge when constructed to order $p>{\lceil\nu\rceil-1}$. So in effect, for practical purposes, we say that $K$ is $\lceil\nu\rceil-1$ times differentiable. Of particular interest are the cases:

[TABLE]

where $d\in{\mathbb{Z}}_{\infty}$, where the latter is simply the RBF kernel.

We postulate the following:

Postulate IV

For all $d\in{\mathbb{Z}}_{\infty}$ the ratio function:

[TABLE]

has only three stationary points - one local maxima $R_{d}\left(0\right)=1$ and two local minima $R_{d}\left(\pm\tilde{r}_{d}\right)<1$ - and in the limits $\mathop{\lim}\limits_{r\to\pm\infty}R_{d}\left(r\right)=\infty$. Furthermore, defining $\beta_{d}=\mathop{\min}\limits_{r}R_{d}(r)=R_{d}(\pm\tilde{r}_{d})$ $\forall d\in{\mathbb{Z}}_{\infty}$:

[TABLE]

Table 2 gives $\nu_{d}$ for $d\in{\mathbb{Z}}_{66}$ (obtained by simulation).

Discussion:

We have not been able to prove this postulate. Figure 4 shows $R_{d}(r)$ for $d=0,1,2$, which conforms to the postulate, and we have simulated (and confirmed) the postulate up to $d=83$ (at which point we ran into floating point problems due to the large factorials involved). We note that this far exceeds practical requirements - most practitioners consider only $\nu\in\{\frac{1}{2},\frac{3}{2},\frac{5}{2}\}$ (i.e. $d\in\{0,1,2\}$). $\square$

Assuming the postulate is correct we have the following result:

Theorem V

For all $d\in{\mathbb{Z}}_{\infty}\backslash\{0\}$, $q\in{\mathbb{Z}}_{d}+1$, $0\leq r\leq\frac{1}{2}M^{2}$:

[TABLE]

where:

[TABLE]

Proof:

Start with (17):

[TABLE]

and hence:

[TABLE]

so by postulate IV it follows that, as $0\leq r\leq\frac{1}{2}M^{2}$:

[TABLE]

and, using postulate IV and recalling that $\Gamma(d+\frac{1}{2})\geq\Gamma(\frac{3}{2})=\frac{1}{2}\sqrt{\pi}$ $\forall d\in{\mathbb{Z}}_{\infty}\backslash\{0\}$, $q\in{\mathbb{Z}}_{d}+1$:

[TABLE]

and so:

[TABLE]

where:

[TABLE]

As $\Gamma(d+\frac{1}{2})\geq\Gamma(\frac{3}{2})=\frac{1}{2}\sqrt{\pi}$ $\forall d\in{\mathbb{Z}}_{\infty}\backslash\{0\}$ for $q\in{\mathbb{Z}}_{d}+1$ we have:

[TABLE]

and the result follows. $\square$

Finally we note that the remainders of the Taylor expansion of the Matérn kernels do not converge as the derivatives exist only to finite order (excepting the case $\nu\to\infty$, which corresponds to the RBF kernel). The remainders $\Delta_{(d+\frac{1}{2},\rho)r}(\delta r)$ and $\Delta_{(d+\frac{1}{2},\rho)}(\delta r)$ appear non-trivial to calculate, and we have been unable to find a closed-form bound. Section D discusses how these may be approximated.

C.3 A Counter-Example: the Rational Quadratic Kernel

The rational quadric kernel is defined by [5]:

[TABLE]

where $\theta\in{\mathbb{R}}_{+}$ and:

[TABLE]

We see immediately that $K\in{\mathbb{C}}_{\infty}$, and $\forall q\in{\mathbb{Z}}_{\infty}$:

[TABLE]

However when we attempt to find $L_{(\theta)}^{\uparrow}$, $L_{(\theta)}^{\downarrow}$ to satisfy assumption 4 - that is, $L_{(\theta)}^{\uparrow}$, $L_{(\theta)}^{\downarrow}$ satisfying:

[TABLE]

we immediately see that no such can exist, as the central term grows factorially with $q$, so $L_{(\theta)}^{\uparrow},L_{(\theta)}^{\downarrow}\to\infty$; and while we may artificially bound the differentiability to $s<\infty$, the resulting bound on the remainders of the Taylor expansion of $f$ (see theorem 3) will grow exponentially with $s$ as:

[TABLE]

rendering the bound useless in this particular case and making the rational quadratic kernel unsuitable for our purposes here.

Appendix D A Note on the Estimation of the Intrinsic Remainders $\Delta_{r}\left(\delta r\right)$ and $\Delta\left(r\right)$ for Non-Convergent Kernels

In the previous section the Matérn kernels discussed have non-convergent Taylor expansions and hence non-zero intrinsic remainders $\Delta_{r}(\delta r)$ and $\Delta(\delta r)$ (that is, one cannot obtain an arbitrarily accurate approximation of $\kappa(r+\delta r)$ by simply extending the Taylor series about $r$ to arbitrary order). We also noted that these remainders may be difficult to bound (tightly) in closed-form. In this section we discuss how they may be approximated using a simple Monte-Carlo approach [4].

We proceed as follows. Let $K({\bf x},{\bf x}^{\prime})=\kappa(\frac{1}{2}\|{\bf x}-{\bf x}^{\prime}\|_{2}^{2})$ be an $s$-times differentiable isotropic kernel, where $s$ is finite. Define:

[TABLE]

to be the actual (tight) remainder bound on the Taylor expansion of $\kappa$ about $r$ to maximal order $s$ (for example if we are using a Matérn kernel of order $\frac{5}{2}$ then $s=2$, so this is easily calculated).

The intrinsic remainder bound $\Delta_{r}(\delta r)$ must satisfy:

Remainder bound: $\Delta_{r}(\delta r)\geq E_{r}(\delta r)$ $\forall\delta r\geq 0$. 2. 2.

Non-decreasing: $\Delta_{r}(\delta r)\geq\Delta_{r}(\delta s)$ $\forall\delta s\in[0,\delta r]$.

It follows that we may estimate a lower bound on $\Delta_{r}(\delta r)$ by sampling:

[TABLE]

where the number of samples $R_{A}$ controls the accuracy of this estimate. Moreover we can use the same approach to approximate $\Delta(\delta r)$:

[TABLE]

where the number of meta-samples $R_{B}$ controls the accuracy of this estimate along with $R_{A}$.

The total number of samples required to approximate $\Delta(\delta r)$ in this scheme is $R_{A}R_{B}$. This may seem large if good accuracy is desired, but it is only required twice in the Bayesian optimisation algorithm - once to calculate $p_{\rm min}$, once to calculate $U_{p}(B)$ - so in practice this presents no real difficulty.

Appendix E Proof of Theorem 5

Before proceding with the proof of theorem 5 we first prove the Theorem:

Theorem VI** (Expected Stable Gain)**

Let ${\mathbb{D}}=\{({\bf x}_{i},{y}_{i})|{y}_{i}=f({\bf x}_{i})+\epsilon_{i}\}$, and ${\mathbb{F}}=\{(\tilde{\bf x}_{i},\tilde{y}_{i})|\tilde{y}_{i}=f(\tilde{\bf x}_{i})+\tilde{\epsilon}_{i}\}$. Assume without loss of generality that $\tilde{y}_{0}\leq\tilde{y}_{1}\leq\ldots$ and define $\tilde{y}_{-1}=\chi$. Given:

[TABLE]

the expected stable gain of ${\mathbb{F}}$ given ${\mathbb{D}}$ is:

[TABLE]

Proof:

Using the fact that $f$ is a draw from a Gaussian Process, defining {\bf v}_{j}^{(q)}\sim{\mathcal{N}}(\frac{B^{q}}{q!}{\bf m}_{\mathbb{D}}^{(q)}(\tilde{\bf x}_{j}),(\frac{B^{q}}{q!})^{2}{\mbox{\boldmath\Lambda}}_{\mathbb{D}}^{(q)}(\tilde{\bf x}_{j},\tilde{\bf x}_{j})) $\forall j$ and applying Theorem 2:

[TABLE]

where the range of the product arises from the assumed ordering on $\tilde{y}_{i}$. $\square$

Having established the preliminary we now prove the theorem:

Theorem 5

Let ${\mathbb{D}}=\{({\bf x}_{i},{y}_{i})|{y}_{i}=f({\bf x}_{i})+\epsilon_{i}\}$. Assume without loss of generality that ${y}_{0}\leq{y}_{1}\leq\ldots$ and define $y_{-1}=\chi$, $y_{|{\mathbb{D}}|}=\infty$. Given:

[TABLE]

the EISG acquisition function reduces to:

[TABLE]

where $z_{i}({\bf x})=\frac{m_{\mathbb{D}}({\bf x})-y_{i}}{\lambda^{1/2}_{\mathbb{D}}({\bf x},{\bf x})}$, $\Delta\hat{y}_{i}({\bf x})=\frac{{y}_{i}-{y}_{i-1}}{\lambda^{1/2}_{{\mathbb{D}}}({\bf x},{\bf x})}$ and:

[TABLE]

so $\Delta\Phi_{|{\mathbb{D}}|}({\bf x})=\Phi(z_{|{\mathbb{D}}|-1}({\bf x}))$ and $\Delta\phi_{|{\mathbb{D}}|}({\bf x})=\phi(z_{|{\mathbb{D}}|-1}({\bf x}))$. The weights $\omega_{0}$, $\omega_{1}$, $\ldots$, $\omega_{|{\mathbb{D}}|}$ are given by:

[TABLE]

Proof:

Working from the definition of EISG:

[TABLE]

where the outer expectation is with regard to $f({\bf x})\sim{\mathcal{N}}(m_{\mathbb{D}}({\bf x}),\lambda_{\mathbb{D}}({\bf x},{\bf x}))$ and the inner expectation with regard to \nabla_{\bf x}^{(q)}f({\bf x})\sim{\mathcal{N}}({\bf m}_{\mathbb{D}}^{(q)}({\bf x}),{\mbox{\boldmath\Lambda}}^{(q)}({\bf x},{\bf x})), and we have defined:

[TABLE]

Using Lemma VI:

[TABLE]

and hence $\forall k\in{\mathbb{Z}}_{|{\mathbb{D}}|+1}$:

[TABLE]

and:

[TABLE]

and the first result follows by summing over all $k$.

The recursive form of the weights $\omega_{i}$ may be deduced by inspection, using the convention that the empty product evaluates to $1$. $\square$