Bayesian optimization in ab initio nuclear physics

TL;DR

This paper investigates the use of Bayesian optimization to efficiently determine coupling constants in complex nuclear physics models, especially when traditional methods are computationally expensive.

Contribution

It demonstrates that Bayesian optimization is effective for low-dimensional parameter spaces in nuclear physics models, aiding in the calibration of coupling constants.

Findings

Bayesian optimization performs well with low-dimensional parameters.

It can be useful for optimizing a small set of coupling constants.

Potential application in determining three-nucleon forces.

Abstract

Theoretical models of the strong nuclear interaction contain unknown coupling constants (parameters) that must be determined using a pool of calibration data. In cases where the models are complex, leading to time consuming calculations, it is particularly challenging to systematically search the corresponding parameter domain for the best fit to the data. In this paper, we explore the prospect of applying Bayesian optimization to constrain the coupling constants in chiral effective field theory descriptions of the nuclear interaction. We find that Bayesian optimization performs rather well with low-dimensional parameter domains and foresee that it can be particularly useful for optimization of a smaller set of coupling constants. A specific example could be the determination of leading three-nucleon forces using data from finite nuclei or three-nucleon scattering experiments.

Tables

Table 1. Table 1: Three different parameter domains (small, medium, large) of the 12 parameters (LECs) in the nuclear physics objective function we study here. The LECs govern the short-range contact potential ( C ~ ⋅ subscript ~ 𝐶 ⋅ \tilde{C}_{\cdot} ) at leading order (unit: 10 4  GeV -2 ), the short-range contact potential ( C ⋅ subscript 𝐶 ⋅ C_{\cdot} ) at next-to-leading order (unit: 10 4  GeV -4 ), and the sub-leading long-ranged pion potential ( c ⋅ subscript 𝑐 ⋅ c_{\cdot} ) at next-to-next-to-leading order (unit: GeV -1 ).

	Domains
parameter	small	medium	large
${\tilde{C}}_{{}_{}{}^{1}S_0^{}}^{\left(np\right)}$	(-0.2,-0.1)	(-5,+5)	(-5,+5)
${\tilde{C}}_{{}_{}{}^{3}S_1^{}}$	(+2,+3)	(-5,+5)	(-5,+5)
$C_{{}_{}{}^{1}S_0^{}}$	(-0.2,-0.1)	(-5,+5)	(-5,+5)
$C_{{}_{}{}^{3}S_1^{}}$	(-1,+1)	(-5,+5)	(-5,+5)
$C_{{}_{}{}^{3}P_0^{}}$	(-1,+1)	(-5,+5)	(-5,+5)
$C_{{}_{}{}^{3}P_1^{}}$	(-1,+1)	(-5,+5)	(-5,+5)
$C_{{}_{}{}^{3}P_2^{}}$	(-1,+1)	(-5,+5)	(-5,+5)
$C_{{}_{}{}^{1}P_1^{}}$	(-1,+1)	(-5,+5)	(-5,+5)
$C_{E_1}$	(-1,+1)	(-5,+5)	(-5,+5)
$c_1$	(-0.76,-0.72)	(-0.76,-0.72)	(-5,+5)
$c_3$	(-3.66,-3.56)	(-3.66,-3.56)	(-5,+5)
$c_4$	(+2.41,+2.47)	(+2.41,+2.47)	(-5,+5)

Equations

$${\mathbf{x}}_{\star }=\underset{\mathbf{x}\in \mathcal{X}}{\text{arg min}},f(\mathbf{x}),$$

$$\left[\begin{array}[]{c}y_{1}\\ \vdots\\ y_{n}\end{array}\right]\sim\mathcal{N}\left(\left[\begin{array}[]{c}\mu_{0}(\mathbf{x}_{1})\\ \vdots\\ \mu_{0}(\mathbf{x}_{n})\end{array}\right],\mathbf{K}=\left[\begin{array}[]{ccc}k(\mathbf{x}_{1},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{1},\mathbf{x}_{n})\\ \vdots&\ddots&\vdots\\ k(\mathbf{x}_{n},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{n},\mathbf{x}_{n})\end{array}\right]\right),$$

$$\displaystyle\begin{split}{}&\mu(\mathbf{x})=\mu_{0}(\mathbf{x})+\left[\begin{array}[]{c}k(\mathbf{x_{1}},\mathbf{x})\\ \vdots\\ k(\mathbf{x_{n}},\mathbf{x})\end{array}\right]^{T}\left[\begin{array}[]{ccc}k(\mathbf{x}_{1},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{1},\mathbf{x}_{n})\\ \vdots&\ddots&\vdots\\ k(\mathbf{x}_{n},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{n},\mathbf{x}_{n})\end{array}\right]^{-1}\left[\begin{array}[]{c}y_{1}-\mu_{0}(\mathbf{x}_{1})\\ \vdots\\ y_{n}-\mu_{0}({\mathbf{x}_{n}})\end{array}\right]\\ {}&\sigma(\mathbf{x})^{2}=k(\mathbf{x},\mathbf{x})-\left[\begin{array}[]{c}k(\mathbf{x_{1}},\mathbf{x})\\ \vdots\\ k(\mathbf{x_{n}},\mathbf{x})\end{array}\right]^{T}\left[\begin{array}[]{ccc}k(\mathbf{x}_{1},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{1},\mathbf{x}_{n})\\ \vdots&\ddots&\vdots\\ k(\mathbf{x}_{n},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{n},\mathbf{x}_{n})\end{array}\right]^{-1}\left[\begin{array}[]{c}k(\mathbf{x_{1}},\mathbf{x})\\ \vdots\\ k(\mathbf{x_{n}},\mathbf{x})\end{array}\right].\end{split}$$

$$\begin{array}{rl}{\mathcal{A}}_{\mathrm{EI}}(\mathbf{x})=& ⟨\max (0,f_{\min }-f(\mathbf{x}))⟩={\int }_{-\infty }^{\infty }\max (0,f_{min}-f)\mathcal{N}(f(\mathbf{x});\mu (\mathbf{x}),\sigma (\mathbf{x}{)}^2)df(\mathbf{x})=\\ & {\int }_{-\infty }^{f_{\min }}(f_{\min }-f)\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{f-\mu }{2{\sigma }^2}\right]df=\\ & (f_{\min }-\mu )\Phi \left(\frac{f_{\min }-\mu }{\sigma }\right)+\sigma \varphi \left(\frac{f_{\min }-\mu }{\sigma }\right)=\sigma \left[z\Phi (z)+\varphi (z)\right],\end{array}$$

$$\mathcal{A}(\mathbf{x}{)}_{\mathrm{LCB}}=\beta \sigma (\mathbf{x})-\mu (\mathbf{x}).$$

$$\begin{array}{rl}f(x,y)=& -\sum _{i=1}^5\left\{\frac{c_i\cos \left[\pi ((x-a_i{)}^2+(y-b_i{)}^2)\right]}{\exp \left[\frac{1}{\pi }((x-a_i{)}^2+(y-b_i{)}^2)\right]}\right\},\\ & \mathbf{a}=(3,5,2,1,7),\\ & \mathbf{b}=(5,2,1,4,9),\\ & \mathbf{c}=(1,2,5,2,3).\end{array}$$

$$d_s(\alpha )=\frac{1}{|\mathcal{P}|}\mathrm{size}\left\{p\in \mathcal{P}:\frac{t_{s,p}}{D_p+1}\le \alpha \right\}.$$

$$f(x_0)-f(x)\ge (1-\tau )(f(x_0)-f_L).$$

$$P(N,D)=2^{-N+1}\sum_{k=0}^{D-1}\left(\begin{array}[]{c}N-1\\ k\end{array}\right).$$

$$f(\mathbf{x})=\sum _{g=1}^{N_G}\underset{{\nu }_g}{\min }\left\{\sum _{i=1}^{N_{g,d}}{\left(\frac{{\nu }_g{\mathcal{O}}_{g,i}^{\mathrm{model}}(\mathbf{x})-{\mathcal{O}}_{g,i}^{\mathrm{experiment}}}{{\sigma }_{g,i}}\right)}^2+{\left(\frac{1-{\nu }_g}{{\sigma }_{g,0}}\right)}^2\right\}.$$

$$f(\mathbf{x})=\frac{1}{2}\sum _{i=1}^p{R}_i(\mathbf{x}{)}^2.$$

$$m_k({\mathbf{x}}_k+\mathbf{s})=f({\mathbf{x}}_k)+{\mathbf{g}}_k^T\mathbf{s}+\frac{1}{2}{\mathbf{s}}^T{\mathbf{H}}_k\mathbf{s}$$

$$q_k^{(i)}(\mathbf{x})=c^{(i)}+(\mathbf{x}-{\mathbf{x}}^k{)}^T{\mathbf{g}}^{(i)}+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^k{)}^T{\mathbf{H}}^{(i)}(\mathbf{x}-{\mathbf{x}}^k)$$

$$M_k({\mathbf{x}}_k+\mathbf{s})=f({\mathbf{x}}_k)+{\mathbf{s}}^T\sum _{i=1}^p{R}_i({\mathbf{x}}_k){\mathbf{g}}^{(i)}+\frac{1}{2}{\mathbf{s}}^T\sum _{i=1}^p\left[{\mathbf{g}}^{(i)}{\left({\mathbf{g}}^{(i)}\right)}^T+R_i({\mathbf{x}}_k){\mathbf{H}}^{(i)}\right]\mathbf{s}.$$

$$B_k=\left\{\mathbf{x}\in {\mathbb{R}}^D:||\mathbf{x}-{\mathbf{x}}_k||\le {\delta }_k\right\}.$$

$$\begin{array}{rl}& f(\mathbf{x})=-A\exp \left[\frac{B}{d}\sum _{i=1}^d{x}_i^2\right]-\exp \left[\frac{1}{d}\sum _{i=1}^d\cos (C{x}_i)\right]+A+e,\\ & \text{where }A=20,B=\frac{1}{2},\text{ and }C=2\pi .\\ & \text{Domain: }{x}_i\in [-30,+30],\text{ for }i=1,\dots ,d.\\ & \text{Global minimum: }{\mathbf{x}}_{\star }=(0,\dots ,0),\text{ and }f({\mathbf{x}}_{\star })=0.\end{array}$$

$$\begin{array}{rl}& f(\mathbf{x})=-{\left[\frac{1}{d}\sum _{i=1}^d{g}_i(x_i)\right]}^2,\\ & \text{where }{g}_i(x_i)=\{\begin{array}{ll}-\frac{x_i}{{\alpha }_i}+\frac{4}{5}& \text{if }0\le x_i<\frac{4}{5}{\alpha }_i\\ +5\frac{x_i}{{\alpha }_i}-4& \text{if }\frac{4}{5}{\alpha }_i<x_i\le {\alpha }_i\\ +5\frac{x_i-{\alpha }_i}{{\alpha }_i-1}+1& \text{if }{\alpha }_i<x_i\le \frac{1+4{\alpha }_i}{5}\\ +\frac{x_i-1}{1-{\alpha }_i}+\frac{4}{5}& \text{if }\frac{1+4{\alpha }_i}{5}<x_i\le 1\end{array}\\ & \text{and }{\alpha }_i=\frac{i}{d+1}.\\ & \text{Domain: }{x}_i\in [0,1],\text{ for }i=1,\dots ,d.\\ & \text{Global minimum: }{\mathbf{x}}_{\star }=({\alpha }_1,\dots ,{\alpha }_d),\text{ and }f({\mathbf{x}}_{\star })=-1.\end{array}$$

$$\begin{array}{rl}& f(\mathbf{x})=10n+\sum _{i=1}^d\left[x_i^2-10\cos (2\pi x_i)\right].\\ & \text{Domain: }{x}_i\in [-5.12,5.12],\text{ for }i=1,\dots ,d.\\ & \text{Global minimum: }{\mathbf{x}}_{\star }=(0,\dots ,0),\text{ and }f({\mathbf{x}}_{\star })=0.\end{array}$$

$$\begin{array}{rl}& f(\mathbf{x})=\sum _{i=1}^{d-1}\left[100(x_{i+1}-x_i^2{)}^2+(1-x_i{)}^2\right].\\ & \text{Domain: }{x}_i\in [-2.048,2.048],\text{ for }i=1,\dots ,d.\\ & \text{Global minimum: }{\mathbf{x}}_{\star }=(1,\dots ,1),\text{ and }f({\mathbf{x}}_{\star })=0.\end{array}$$

$$\begin{array}{rl}& f(\mathbf{x})=418.9829d-\sum _{i=1}^d\left[x_i-\sin \left(\sqrt{|x_i|}\right)\right].\\ & \text{Domain: }{x}_i\in [-500,500],\text{ for }i=1,\dots ,d.\\ & \text{Global minimum: }{\mathbf{x}}_{\star }=(420.9687,\dots ,420.9687),\text{ and }f({\mathbf{x}}_{\star })\approx 0.\end{array}$$

$$\begin{array}{rl}& f(\mathbf{x})=\sum _{i=1}^d{x}_i^2.\\ & \text{Domain: }{x}_i\in [-5.12,5.12],\text{ for }i=1,\dots ,d.\\ & \text{Global minimum: }{\mathbf{x}}_{\star }=(0,\dots ,0),\text{ and }f({\mathbf{x}}_{\star })=0.\end{array}$$

$$V({\mathbf{p}}_{\mathbf{f}},{\mathbf{p}}_{\mathbf{i}})=V_{\mathrm{long}-\mathrm{range}}({\mathbf{p}}_{\mathbf{f}},{\mathbf{p}}_{\mathbf{i}})+V_{\mathrm{contact}}({\mathbf{p}}_{\mathbf{f}},{\mathbf{p}}_{\mathbf{i}}),$$

$$V_{1\pi}^{(\rm LO)}(\mathbf{p_{f}},\mathbf{p_{i}})=\frac{g_{A}^{2}}{4f_{\pi}^{2}}\>\mbox{\boldmath$\tau$}_{1}\cdot\mbox{\boldmath$\tau$}_{2}\>\frac{\mbox{\boldmath$\sigma$}_{1}\cdot\mathbf{q}\,\,\mbox{\boldmath$\sigma$}_{2}\cdot\mathbf{q}}{q^{2}+m_{\pi}^{2}}\,,$$

$$\displaystyle\begin{split}{}&V_{2\pi}^{({\rm NLO})}(\mathbf{p_{f}},\mathbf{p_{i}})+V_{2\pi}^{({\rm NNLO})}(\mathbf{p_{f}},\mathbf{p_{i}})=\\ {}&L(q;\tilde{\Lambda})\left[(\mbox{\boldmath$\tau$}_{1}\cdot\mbox{\boldmath$\tau$}_{2})\left\{4m_{\pi}^{2}(1+4g_{A}^{2}-5g_{A}^{4})+q^{2}(1+10g_{A}^{2}-23g_{A}^{4})-\frac{48g_{A}^{4}m_{\pi}^{4}}{w^{2}}\right\}\right.-\\ {}&\left.18g_{A}^{4}((\mbox{\boldmath$\sigma$}_{1}\cdot\mathbf{q}\,\,\mbox{\boldmath$\sigma$}_{2}\cdot\mathbf{q})-q^{2}(\mbox{\boldmath$\sigma$}_{1}\cdot\mbox{\boldmath$\sigma$}_{2}))\right]+A(q;\tilde{\Lambda})\left[(2m_{\pi}^{2}+q^{2})(2m_{\pi}^{2}(c_{3}-2c_{1})+c_{3}q^{2})\right.-\\ {}&\frac{1}{6}c_{4}w^{2}((\mbox{\boldmath$\tau$}_{1}\cdot\mbox{\boldmath$\tau$}_{2})\,(\mbox{\boldmath$\sigma$}_{1}\cdot\mathbf{q}\,\,\mbox{\boldmath$\sigma$}_{2}\cdot\mathbf{q})-\left.q^{2}(\mbox{\boldmath$\tau$}_{1}\cdot\mbox{\boldmath$\tau$}_{2})(\mbox{\boldmath$\sigma$}_{1}\cdot\mbox{\boldmath$\sigma$}_{2}))\right]\end{split}$$

$$\begin{array}{rl}& L(q;\tilde{\Lambda })=\frac{w}{768{\pi }^2{f}_{\pi }^{2}q}\ln \frac{{\tilde{\Lambda }}^2(2m_{\pi }^2+q^2)-2m_{\pi }^2{q}^2+\tilde{\Lambda }\sqrt{{\tilde{\Lambda }}^2-4m_{\pi }^2}qw}{2m_{\pi }^2({\tilde{\Lambda }}^2+q^2)}\\ & A(q;\tilde{\Lambda })=\frac{3g_A^2}{32\pi f_{\pi }^{4}q}\arctan \frac{q(\tilde{\Lambda }-2m_{\pi })}{q^2+2\tilde{\Lambda }{m}_{\pi }}\end{array}$$

$$\displaystyle\begin{split}V_{\rm contact}(\mathbf{p_{f}},\mathbf{p_{i}})={}&C_{S}+C_{T}(\mbox{\boldmath$\sigma$}_{1}\cdot\mbox{\boldmath$\sigma$}_{2})+C_{1}q^{2}+C_{2}k^{2}+(C_{3}q^{2}+C_{4}k^{2})(\mbox{\boldmath$\sigma$}_{1}\cdot\mbox{\boldmath$\sigma$}_{2})+\\ {}&C_{5}(-i\mathbf{S}\cdot(\mathbf{q}\times\mathbf{k}))+C_{6}(\mbox{\boldmath$\sigma$}_{1}\cdot\mathbf{q})(\mbox{\boldmath$\sigma$}_{2}\cdot\mathbf{q})+C_{7}(\mbox{\boldmath$\sigma$}_{1}\cdot\mathbf{k})(\mbox{\boldmath$\sigma$}_{2}\cdot\mathbf{k})\end{split}$$

Full text

1 Introduction

Mathematical optimization plays a central role in natural science. Indeed, most theoretical predictions are preceded by a calibration stage whereby the parameters of the model are optimized to reproduce a selected set of calibration data. In nuclear physics, the coupling constants of any theory of the strong interaction between protons and neutrons (nucleons) must be determined from experimental data before one can attempt to solve the Schrödinger equation to make quantitative predictions of the properties of atomic nuclei.

Typically, measured low-energy nucleon-nucleon ($NN$) cross sections and the properties of light nuclei with mass number $A\leq 4$ have been used for calibrating the $NN$ and three-nucleon ($NNN$) interaction sectors of the nuclear force, see e.g. Refs. [1, 2, 3] and references therein. However, it is still an open question—with several parallel and ongoing research efforts [4, 5, 6, 7]—how to best constrain the unknown coupling constants in theoretical descriptions of the nuclear interaction and how to incorporate the covariance structure of the experimental data and the model discrepancy [5, 8, 9, 10]. Recent theoretical studies [7, 11, 12, 13] indicate that $NN$ scattering data and the properties of very light nuclei do not contain enough information to constrain all directions in the parameter domain at sufficient level. Instead, it has been proposed that the pool of fit data need to be augmented with complimentary data types, such as $NNN$ scattering cross sections, and/or the ground-state properties of nuclei with $A>4$, or even empirical properties of infinite nuclear matter. However, modeling of such observables is typically much more complex and requires substantial computational effort ranging from hours to days for just one model evaluation, even on a supercomputer. Consequently, the optimization of the underlying model parameters will be difficult. The main focus of the present work is to investigate a possible strategy for mitigating this computational challenge.

Inspired by recent progress in the optimization of hyperparameters of deep neural networks [14], we explore several Bayesian optimization (BayesOpt) strategies111we employ the BaysOpt implementation provided through the Python package GPyOpt [15]. for maximizing the likelihood of objective functions based on complex models in nuclear physics. BayesOpt originated more than 50 years ago [16], it was popularized in the 1990s, see e.g. [17, 18], and has since then been applied in various disciplines; from selecting experiments in material and drug design [19] to tuning event-generators in particle physics [20]. The central idea in BayesOpt is to construct a probabilistic surrogate model, usually a Gaussian process ($\mathcal{GP}$), to capture our prior beliefs and knowledge about the objective function, $f(\mathbf{x})$, and iteratively confront the surrogate with actual data samples from $f(\mathbf{x})$ and thereby refine our posterior of this function. The main advantage of BayesOpt is that we can incorporate prior beliefs in a straightforward way. The disadvantage lies in the arbitrariness and uncertainty of a priori information.

In the following we will be dealing with complex models in nuclear physics. The origin of the underlying physics model and its parameters is briefly introduced in this section, while more details and relevant references are provided in B. Nucleons are made of quarks and gluons, and it is well known that the strong interaction between these fundamental particles is described in detail by quantum chromodynamics (QCD), which is part of the standard model of particle physics. It is equally well-known that QCD is not perturbative in the low-energy region relevant for nuclear structure physics. This prevents straightforward application of, e.g., perturbation theory to compute atomic nuclei starting from QCD. Instead, chiral effective field theory ($\chi$EFT) [21, 22, 23] is constructed as a low-energy approximation to QCD. This framework shows promising signs of being an operational approach to analyzing atomic nuclei while maintaining a firm link with the more fundamental theory. In $\chi$EFT, the long-ranged part of the nuclear interaction is described in terms of pion exchanges, while the short-ranged part is parameterized by contact interactions. The unknown parameters of this description are known as low-energy constants (LECs). It is important to realize that any realistic description of the nuclear interaction has to introduce unknown coefficients, but it has been found that $\chi$EFT is able to capture the physics with a relatively small number of parameters. Depending on the level of details included, there are typically $\sim 10$–$30$ LECs.

Clearly, if model predictions of physical observables, such as scattering cross sections, are computationally expensive, we face the challenge of optimization with a limited budget of evaluations of the objective function. For the overwhelming majority of well-developed research problems there does not exist a universal optimization strategy that guarantees to arrive at a global optimum of an objective function in a finite number of iterations. Instead, any information regarding the mathematical or computational structure of the objective function, perhaps guided by the physical nature of the underlying problem, should play an important role in the choice or design of the optimization algorithm.

In this paper, we will systematically study the application of BayesOpt to optimization problems of increasing degree of complexity. The BayesOpt algorithm is presented in Sec. 2 with its main ingredients: the $\mathcal{GP}$ kernel and the acquisition function. The performance of different optimization algorithms can be compared using a data profile. This measure, as used in the present work, is introduced in Sec. 3. In order to benchmark the performance of BayesOpt with various settings in controlled problem conditions we will employ a selected set of six test functions in $D=2,4,8$ dimensional parameter domains. This study is presented in Sec. 4, while the test functions themselves are listed in A. The main focus of this work is found in Sec. 5 with the application of BayesOpt to a real nuclear physics problem. We will use BayesOpt to optimize the 12 LECs appearing at next-to-next-to-leading order (NNLO) in $\chi$EFT, using the proton-neutron scattering data in the Granada database [24] with laboratory energies of the incident proton beam below 75 MeV. This case is complex enough to constitute a non-trivial problem from a physics as well as an optimization perspective. However, it is still computationally straightforward such that we can easily afford a detailed analysis of 12 different BayesOpt algorithms. Indeed, each evaluation of the objective function at a specific point in the parameter domain only takes a couple of seconds on a standard desktop computer. Still, evaluating the 4,096 corners in the hypercube of the corresponding parameter domain will take a couple of hours, so the premise for BayesOpt is well justified. We will compare the optimization performance of BayesOpt with the POUNDERs algorithm [25]. This is a simulation-based optimization algorithm that has been successfully applied in non-linear least-squares optimization in nuclear physics before [26, 27]. Finally, we end with a summary and outlook in Sec. 6.

2 Bayesian optimization

Without any loss of generality, we will frame the determination of the LECs in $\chi$EFT as a minimization problem. Global minimization of a function $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$, with input parameters $\mathbf{x}$ that are perhaps subject to some constraints $c(\mathbf{x})\leq 0$ and typically belong to a compact input domain $\mathcal{X}\subset\mathbb{R}^{D}$, is a long-standing and central problem in science. Here, we also specialize the formalism to scalar valued functions $f$. Mathematically, we are trying to find a global minimizer:

[TABLE]

i.e. a point that fulfill $f(\mathbf{x}_{\star})\leq f(\mathbf{x})$ for all $\mathbf{x}\in\mathcal{X}$. In general, this is an intractable problem unless we have detailed information about $\mathbf{x}$ or that the parameter domain only contains a finite number of points. In reality, we are trying to find local minimizers to $f(\mathbf{x})$, i.e. points $\mathbf{x}_{\star}$ for which $f(\mathbf{x}_{\star})\leq f(\mathbf{x})$ for all $\mathbf{x}\in\mathcal{X}$ close to $\mathbf{x}_{\star}$. A substantial complication arises if $f(\mathbf{x})$ is computationally expensive to evaluate. Then it becomes even more important to adopt a well-founded strategy, based on all present knowledge about $f(\mathbf{x})$, for carefully selecting a dataset of $n$ sequential function queries $\mathcal{D}_{1:n}=\{\mathbf{x}_{i},y_{i}\}_{i=1}^{n}$, where $y_{i}=f(\mathbf{x}_{i})$, such that we increase our chances of a rapid convergence towards $\mathbf{x}_{\star}$. A BayesOpt framework can provide this. Moreover, at each iteration BayesOpt will consider all insofar collected data points and thereby take full advantage of the history of the optimization run. Note that we refer to a set of function evaluations as data. This should not be confused with experimental data. There are two main components in a BayesOpt algorithm;

• •

A prior probabilistic belief $p(f|\mathcal{D})$ for the function $f$ given some data $\mathcal{D}$. The prior is often a $\mathcal{GP}$. This is updated in every iteration.

• •

An acquisition function $\mathcal{A}(\mathbf{x}|\mathcal{D})$ given some data $\mathcal{D}$, i.e. a heuristic that balances exploration against exploitation and determines where to evaluate the objective function $f(\mathbf{x})$ next.

The next iterate, $\mathbf{x}_{i+1}$, is selected where we expect the minimizer $\mathbf{x}_{\star}$, based on some utility function. Below, we will define two different acquisition functions $\mathcal{A}(\mathbf{x}|\mathcal{D})$, and show how to embed them in an iterative context for selecting sample points $\mathbf{x}_{i}$. In the following we will also drop the explicit data dependence in the notation for the acquisition function and only write $\mathcal{A}(\mathbf{x})$. A pseudo-code for BayesOpt is listed in Algorithm 1 and a pictorial exposition of a handful of BayesOpt iterations of a simple univariate function is provided in Fig. 1.

2.1 The prior: a Gaussian process

To model the prior $p(f|\mathcal{D})$ for the objective function we use a Gaussian process $\mathcal{GP}(\mathbf{x})$ with mean function $\mu(\mathbf{x})$ and covariance matrix $\mathbf{K}$ with entries $k_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$. The mean and covariance functions fulfill expected relations; $\mu(\mathbf{x})=\mathbb{E}[\mathbf{x}]$ and $k(\mathbf{x},\mathbf{x}^{\prime})=\mathbb{E}[(\mathbf{x}-\mu(\mathbf{x}))(\mathbf{x}^{\prime}-\mu(\mathbf{x}^{\prime}))]$ for all $(\mathbf{x},\mathbf{x}^{\prime})\in\mathcal{X}$. Any real-valued function $\mu(\cdot)$ is permissible, but for $k(\cdot,\cdot)$ the corresponding covariance matrix $\mathbf{K}$ must be positive semidefinite. A $\mathcal{GP}$ is one example of a stochastic process that is very useful in statistical modeling [28]. In brief, it is a collection of function evaluations $y_{1:n}$ at $\mathbf{x}_{1:n}$, with mean $\mu_{0}$ (often shifted to zero), that are jointly Gaussian, i.e.

[TABLE]

where $y\sim\mathcal{N}(\mu,\sigma^{2})$ denotes a normally distributed random variable $y$ with mean $\mu$ and covariance $\sigma^{2}$. After conditioning this prior with some data $\mathcal{D}_{n}=\left\{(\mathbf{x}_{1},y_{1}),\ldots,(\mathbf{x}_{n},y_{n})\right\}$, with mean $\mu_{0}$, we obtain the mean and covariance of the $\mathcal{GP}$ model for the prior according to

[TABLE]

These explicit expressions follow from the fact that the marginal distribution of a multivariate Gaussian is also Gaussian. The mean and variance of the prior, i.e. the $\mathcal{GP}$, for the schematic example in Fig. 1 are indicated with a green line and a green shaded band, respectively.

In this work we will consider three common types of covariance functions, also referred to as kernels:

• •

Squared exponential: $k(\mathbf{x},\mathbf{x}^{\prime})=\theta_{0}^{2}\exp\left(-\frac{1}{2}r^{2}\right)$

• •

Matern 3/2: $k(\mathbf{x},\mathbf{x}^{\prime})=\theta_{0}^{2}(1+\sqrt{3}r)\exp(-\sqrt{3}r)$

• •

Matern 5/2: $k(\mathbf{x},\mathbf{x}^{\prime})=\theta_{0}^{2}(1+\sqrt{5}r+\frac{5}{3}r^{2})\exp\left(-\sqrt{5}r\right)$

where $r^{2}\equiv\sum_{i=1}^{D}\frac{1}{\theta_{i}^{2}}(x_{i}-x_{i}^{\prime})^{2}$ and $\left\{\theta_{i}\right\}_{i=0}^{D}$ is a set of hyperparameters for the kernel. The correlation length(s) $\ell=(\theta_{1},\ldots,\theta_{D})$ indicate how far you need to move (along a particular direction in the parameter domain) for the function values to become uncorrelated. With automatic relevance determination (ARD), we optimize the vector of correlation lengths separately in each direction of the function domain $\mathcal{X}\subset\mathbb{R}^{D}$. Without ARD, the kernel hyperparameters are isotropic. In this work, we extract the hyperparameters using maximum likelihood estimation.

The characteristic features of the $\mathcal{GP}$s based on the three kernels listed above are illustrated in Fig. 2. Their smoothness, equivalent to the typical correlation length, is one of the key differences between them. This feature of the $\mathcal{GP}$ kernel affects the prior modeling of the objective function, and as such the ensuing performance of BayesOpt. We will see this clearly in the analyses in Secs. 4-5.

2.2 The acquisition function

The acquisition function determines the most likely improvement to the currently best minimizer in the parameter domain. In Fig. 1 the mean of our posterior belief (green line) of the unknown values of the objective function $f$ (black line) is sequentially augmented with one new data point (black dot) in each iteration. The best candidate for further minimizing $f$ at iteration $n$ is the parameter that maximizes the acquisition function (red curve). Although the acquisition function is also optimized in $\mathcal{X}\subset\mathbb{R}^{D}$ it is significantly faster to evaluate, compared to the underlying objective function $f(\mathbf{x})$, since it only relies on draws from the prior $\mathcal{GP}$. Still, the complexity of maximizing $\mathcal{A}(\mathbf{x})$ increases as we increase the dimensionality of the parameter domain $\mathcal{X}$. This aspect should not be underestimated, and it is in fact one of the main challenges with BayesOpt. Another challenge, although unlikely, could emerge if the set of collected data points $\mathcal{D}_{1:n}$ becomes very large. Indeed, in each iteration the evaluation of the $\mathcal{GP}$ requires an inversion of an $n\times n$ matrix, and the complexity of that operation naively scales as $\mathcal{O}(n^{3})$. Cholesky decomposition reduces this cost somewhat to $\mathcal{O}(n^{3}/6)$. In reality, however, this is rarely a limiting factor since we typically resort to BayesOpt when only a small number of function evaluations can be carried out in the first place.

A very appealing feature of BayesOpt is its ability to select a new point $\mathbf{x}_{n}$ in a region of $\mathcal{X}$ where the prior model of $f$ is exhibiting a large uncertainty. This means that the algorithm can be rather explorative and therefore escape a local minimum of the objective function $f$. Depending on the details of the acquisition function, the explorative nature is balanced with a certain degree of exploitation, i.e. to evaluate points in the parameter domain where the prior model for $f$ is exhibiting a low mean value. To study the exploration-exploitation balance we will consider two of the most common acquisition functions; the expected improvement (EI) and the lower confidence-bound (LCB). In the following, we denote by $f_{\rm min}$ the insofar lowest recorded value of $f(\mathbf{x})$.

The expected improvement acquisition function is defined by the expectation value of the rectifier ${\rm max}(0,f_{\rm min}-f(\mathbf{x}))$, i.e. we reward any expected reduction of $f$ in proportion to the reduction $f_{\rm min}-f(\mathbf{x})$. This can be evaluated analytically

[TABLE]

where we dropped the explicit dependence on $\mathbf{x}$ in the third step, and the cumulative distribution function and standard normal distribution are denoted $\Phi$ and $\phi$, respectively. In the last step we write the result in the standard normal variable $z=\frac{f_{\rm min}-\mu}{\sigma}$. BayesOpt will exploit regions of expected improvement when the term $z\Phi(z)$ dominates, while new, unknown regions will be explored when the second term $\phi(z)$ dominates. For the expected improvement acquisition function, the exploration-exploitation balance is entirely determined by the set of observed data $\mathcal{D}_{1:n}$ and the $\mathcal{GP}$ kernel.

The lower confidence-bound acquisition function introduces an additional parameter $\beta$ that explicitly sets the level of exploration

[TABLE]

The maximum of this acquisition function will occur for the maximum of the $\beta$-enlarged confidence envelope of the $\mathcal{GP}$. We use $\beta=2$, which is a very common setting. Larger values of $\beta$ leads to even more explorative BayesOpt algorithms.

Besides being derivative-free—although derivatives can in fact be incorporated [29]—it is in many ways the explorative nature of BayesOpt that is most attractive. This feature is most easily observed in the optimization of a rather complex two-dimensional function with several local minima, such as the Langermann function

[TABLE]

In Fig. 3 we show the sequence ${(x_{i},y_{i})}_{i=0}^{51}$ of 50 evaluation points following two initial points of BayesOpt with an Matern 3/2 kernel, expected improvement acquisition function without ARD.

It should be made clear that, given a limited computational budget, the success of BayesOpt hinges on the choice of kernel and acquisition function. In the example above with a Langermann function, the non-smooth nature of the Matern 3/2 kernel is advantageous compared to, e.g., the squared exponential kernel. The exploration-exploitation balance also leverages the success ratio of BayesOpt. To learn more about BayesOpt, it is obviously instructive to benchmark and compare different BayesOpt algorithms using well-known test functions. For this purpose we first need to select a performance measure for optimization algorithms.

3 Measuring the performance of optimization

algorithms

To analyze the performance of derivative-free optimization algorithms we follow Ref. [30] and define a data profile

[TABLE]

The data profile enables direct comparison between a set of optimization algorithms $\mathcal{S}$, all of which are applied to a set of well-defined optimization problems $\mathcal{P}$. For each $(s,p)\in\mathcal{S}\times\mathcal{P}$, the performance measure $t_{s,p}>0$ denotes the number of function evaluations that are required for optimization algorithm $s$ applied to a problem $p$ to satisfy some convergence criterion. Thus, $d_{s}(\alpha)$ is the fraction of problems that can be solved within $\alpha$ function calls. The performance measure can be further normalized to $D_{p}+1$, where $D_{p}$ denotes the dimensionality of the parameter domain in problem $p$. This is an attempt to account for some of the complexity due to a larger number of parameters in the problem. This dimensional scaling is only approximate and motivated by the $(D_{p}+1)$ function evaluations required to compute a $D$-dimensional simplex, i.e. a $D$-dimensional triangle of function evaluations.

Each combination of starting point and objective function (plus other possible constraints) constitutes a separate problem $p$. The data profile is a monotonically increasing function between zero and one, and a large value of $d_{s}(\alpha)$ for small values of $\alpha$ is better. In line with Ref. [30] we employ a convergence criterion

[TABLE]

This is fulfilled for any $x$ where the initial function value $f(x_{0})$ is reduced $1-\tau$ times the best possible reduction $f(x_{0})-f_{L}$. We will set $f_{L}$ to be the lowest value of $f$ achieved by any solver $s\in\mathcal{S}$. Although we will come close to the true solution $x_{\star}$ in a few cases, it is highly unlikely that any derivative-free solver will arrive at $f(x_{\star})$. For that, one would typically have to resort to gradient-based optimization algorithms. With BayesOpt, we will consider a 90% reduction in the function value as a sign of convergence, i.e. set $\tau=0.1$. For the purpose of comparing BayesOpt with other derivative-free optimization algorithms we will occasionally use $\tau=0.01$. When the objective functions are represented by known test functions (described in Sec. 4 and A), we will set $f_{L}=f(x_{\star})$. Throughout this study we will only allow a maximum of $\alpha=$250 function calls per solver and problem, and set $t_{p,s}=\infty$ in case the convergence criterion is not fulfilled. This upper limit on $\alpha$ should cover most scenarios that would call for the use of BayesOpt.

3.1 Selecting parameter starting points in $\mathbb{R}^{D}$

A multi-dimensional parameter domain $\mathcal{X}\subset\mathbb{R}^{D}$ of the objective function leads to the increasingly likely existence of multiple points corresponding to local and/or global minimizers. In reality, the choice of initial point, i.e. where we start the optimization run, will determine the local minimum that the optimizer converges to. When benchmarking, it will be necessary to start each solver $s$ from $N>1$ starting points in $\mathcal{X}$. In this work, we use an identical set of $N=2D$ starting points for all optimizers in a $D-$dimensional domain. This is somewhat motivated by Wendel’s theorem [31] that states that for $N$ random points on the $D-$sphere, the probability $P(N,D)$ that they all lie on some hemisphere is given by

[TABLE]

Thus for random $N=2D$ points on the $D-$sphere, $P(2D,D)=0.5$. Although we do not distribute points on a hypersphere, we use it to motivate a rule of thumb. It should be pointed out that some authors argue for a slightly more expensive rule of thumb to select $N=10D$ starting points when initializing computer simulations in $\mathbb{R}^{D}$ [32].

A priori, we do not distinguish between different parts of the parameter domain, so we select the $N$ starting points using a space-filling algorithm in the form of a quasi-random Sobol sequence. We employ a Python implementation of this algorithm. The mathematical underpinnings of the Sobol sequence are provided in the original paper [33] and a discussion related to the numerical implementation is given in e.g. Ref. [34]. This is a so-called low-discrepancy sequence, which in fact is the opposite of a random sequence. It is designed to generate each successive sample point as far away as possible from all existing sample points. This tends to sample the space more uniformly than pseudo-random numbers, at least for lower-dimensional domains. We will not go beyond $D=12$ in any part of this work. Although the Sobol sequence can exhibit gaps in multi-dimensional spaces, it has several advantages. In addition to fast generation, a Sobol sequence is straightforward to augment with additional sample points. We remind the reader that Latin hypercube sampling (LHS) [35], which is a different kind of algorithm for generating space-filling samples, is in most cases not easy to augment while preserving the latin hypercube structure. The space-filling differences between Sobol, LHS, and conventional pseudo-random numbers in $\mathbb{R}^{2}$ are illustrated in Fig. 4.

4 Analyzing a set of test functions

Before we tackle an optimization problem in nuclear physics, we will explore and benchmark the performance of BayesOpt on a set of test functions. To this end, we have selected a set of six multivariate and continuous functions $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$, each defined on some domain $\mathcal{X}\subset\mathbb{R}^{D}$. The functions are defined for any $D>0$, but we will only consider $D=2,4,8$. The set of test functions reflects an average of various spatial characteristics. Two-dimensional graphical representations are shown in Figs. 5 and 6, with explicit expressions given in A.

A comparison of two or more optimization algorithms on a finite set of test functions is neither fair nor conclusive. Indeed, although one algorithm (or class of methods) appears to be more successful in finding optima, it is clear from the no free lunch theorem in optimization that the averaged performance of any two algorithms on the set of all possible functions will be the same. Here, we merely set out to compare and analyze the performance of BayesOpt on a limited set of characteristically different continuous test functions.

We now turn our attention to the resulting data profiles when applying BayesOpt for finding the minimizer in each one of these test functions with parameter domains in $D=2,4,8$ dimensions. In total, we will analyze 12 BayesOpt algorithms composed from combining three kernels, two acquisition functions, and with or without ARD. We use a Sobol sequence to initiate each BayesOpt algorithm at $N=2D$ different starting points in each parameter domain in $\mathbb{R}^{D}$. Remember that we refer to each specific combination of starting point and test function as a problem. Combined, with the different solvers (optimization algorithm with specific settings) this is a rather large dataset and we analyze it from several different angles. The data profiles for different versions of BayesOpt applied to all six test functions in $D=2$, see Fig. 7 (top row), indicate that $d(50)\sim 0.7$, i.e. that $\sim 70$% of the problems are converged at the $\tau=0.1$ level within 50 function evaluations. The performances of the expected improvement and the lower confidence-bound acquisition functions are very similar. Although, expected improvement exhibits a slightly better and more coherent performance across all $\mathcal{GP}$ kernels, and at 150 function evaluation more than 80% of the test functions in $D=2$ are converged.

ARD, i.e. to learn each individual length scale hyperparameter in the $\mathcal{GP}$ kernel, is often an efficient way of pruning irrelevant features. We find that exploiting ARD increases the performance of BayesOpt once data from a sufficient number of function evaluations has informed the $\mathcal{GP}$ kernel on the possible anisotropic structure of the objective function. This becomes even clearer if we enforce a stronger convergence criterion with $\tau=0.01$, see Fig. 7 (bottom row). There seems to be a slight advantage of using ARD with the expected improvement acquisition function.

In $D=2$, and for this set of functions, the Matern kernels perform slightly better than the squared exponential. In general, the Matern kernels are better tailored to non-smooth objectives. This becomes even clearer when we study the performance of all BayesOpt algorithms on the Ackley function (one hole on a semi-flat surface with several periodic shallow local minima) in $\mathbb{R}^{D=2,4,8}$, see Fig. 8.

From the characteristics of the data profiles shown in Fig. 8, it is obvious that the exploratory nature of the lower confidence-bound acquisition function is highly advantageous for finding the global minimum hiding on a surface covered with local minima. Clearly, it is important to tailor the BayesOpt acquisition function and underlying $\mathcal{GP}$ kernels to the spatial structure of the objective function. This result also reflects the two fundamental and competing aspects of BayesOpt. If we incorporate prior beliefs, or even knowledge, about the objective function, then BayesOpt will perform rather well. On the other hand, the usefulness of BayesOpt will consequently be limited by the arbitrariness and uncertainty of a priori information. This is further complicated by the fact that we typically resort to BayesOpt when we know very little about the objective function in the first place, since it is computationally expensive to evaluate.

Data profiles for the test functions as we increase the dimensionality of the parameter domain from $D=2$ (in Fig. 7) to $D=4$ and $D=8$ are shown in Fig. 9.

We can conclude that having a larger number of parameters provides a significant challenge to BayesOpt, as for all optimization algorithms. Indeed, in $D=4$ dimensions it takes more than $\alpha=150$ function evaluations to reach a data profile value of $d\sim 0.5$. The more exploratory lower confidence-bound acquisition function exhibits a slightly larger performance spread with respect to different $\mathcal{GP}$ kernels. This becomes even more prominent as we increase the dimensionality of the parameter domain to $D=8$. For the set of test functions that we have employed it is marginally advantageous to be more exploratory for higher dimensional objective functions. We also note that the potential benefit of ARD requires data from more than $\alpha=250$ function evaluations. This is natural since ARD introduces more hyperparameters that need to be determined. For all problems in $D=4,8$ dimensions, BayesOpt with Matern kernels converge faster than BayesOpt with a squared exponential kernel. As above, the squared exponential kernel does not capture the high-frequency modes that are present in some of the functions.

5 Bayesian optimization of the nucleon-nucleon interaction

An EFT of the nuclear interaction essentially corresponds to a low-energy parameterization of QCD in a fashion that is consistent with the symmetries of the more fundamental theory. To devise a true EFT description of the nuclear interaction is a major intellectual challenge in ab initio nuclear theory. Several avenues are currently being explored. The physical underpinnings of the objective function we seek to minimize are described in somewhat more detail in B. Regarding the optimization of the LECs, the inclusion of calibration data from more complex atomic nuclei and low-energy nuclear processes, e.g. $NNN$ scattering, in the objective function are currently being considered in the community. Such extensions of the calibration data clearly increases the information content—but does so at the expense of an increased complexity in the computer modeling.

A specific aim of this work is to analyze the performance of BayesOpt for determining the parameters $\mathbf{x}$ in an EFT model of the $NN$ interaction. For this, we use the proton-neutron sector of an EFT at NNLO with 12 parameters and try to find the vector of model parameters $\mathbf{x}$ that are in agreement with existing experimental data on proton-neutron scattering cross sections. This type of observable is a measure of the probability of an incident particle (a proton) with a given kinetic energy to scatter into a certain solid angle element due to mutual interaction with the target particle (a neutron222Since free neutrons decay within $\sim 10$ minutes, the neutron target is a composite material containing neutrons.). We have deliberately chosen this class of calibration data since it does not render particularly challenging model evaluations. One evaluation of the full objective function $f(\mathbf{x})$ takes a merely couple of seconds on a standard desktop. Still, the complexity of the physical model provides a non-trivial testing ground for assessing nascent applications of BayesOpt in ab initio nuclear physics.

The experimental dataset is composed of $N_{G}$ groups of measurements where the $g$-th group consists of $N_{g,d}$ measured cross sections, with associated random measurement uncertainty, $\mathcal{O}_{g,i}^{\rm experiment}\pm\sigma_{g,i}$, for $i=1,\ldots,d$, with a common normalization constant $\nu_{g}$ and corresponding experimental systematic error $\sigma_{g,0}$. We employ the measurement errors as reported by each experimenter. Each group of data originates from a specific experiment. We restrict the inclusion of experimental data in the present case to laboratory scattering energies below 75 MeV. This ensures a rather fast (seconds) evaluation of the objective function and therefore enables a more detailed analysis of BayesOpt. Experimental errors across measurement groups are considered independent and identically distributed. The normalization constant, together with its uncertainty, represents the systematic uncertainty of the measurement. For an absolute measurement, the normalization is given by $\nu_{g}=1\pm 0$. Usually this means that the statistical and systematic errors have been combined with $\sigma_{g,d}$. Certain experiments are not normalized at all. Instead, only the angular- or energy-dependence of the cross section was determined. For these groups of data, $\nu_{g}$ is solved for in the optimization by minimizing the discrepancy between the model prediction $\mathcal{O}_{g,d}^{\rm model}(\mathbf{x})$ and the experimental data points $\mathcal{O}_{g,d}^{\rm experiment}$. For such freely renormalized data, the optimal $\nu_{g}$ is easily obtained in closed form. For practical purposes, the normalization error can be considered infinite in these cases, and will therefore not contribute to the objective function.

In summary, we seek to find the parameter vector $\mathbf{x}$ that minimizes the deviation between the model and the experimental data, as measured by the objective function $f$ defined below:

[TABLE]

This type of (non-linear) least-squares objective function, where we assume a normal model for the parameters, appears naturally in most parameter estimation problems. In a setting defined by $\chi$EFT it would be natural to also incorporate a theoretical model discrepancy term motivated by the low-energy EFT expansion itself [5, 8, 9, 36, 37]. However, in this paper we will focus on the challenges of mathematical optimization that are associated with a numerically costly objective function, and for simplicity therefore only incorporate the experimental errors of the data.

We define three different parameter domains for minimizing the objective function; small, medium, and large, see Tab. 1. They differ by the level of included prior knowledge regarding the range of plausible parameter values.

The limits of the large parameter domain is based on the most naive estimate without any significant prior information. In contrast, the limits of the medium domain are partly informed by prior data. Specifically, the three parameters (LECs $c_{1},c_{3},c_{4}$) associated with the long-range part of the nuclear interaction are constrained to ranges determined from a separate analysis of pion-nucleon scattering data [38, 39]. The small domain is further informed by previous experience of typical values for the LECs in the short-ranged contact potential.

The data profiles of the BayesOpt algorithms applied to all domains of the nuclear physics objective function in Eq. 10 are plotted in Fig. 10. Clearly, the expected improvement acquisition function performs slightly better than the lower confidence-bound acquisition function. It is difficult to draw any conclusions regarding an optimal choice of kernel. Compared to the test functions, the squared exponential kernel seems to work a little bit better than the Matern 5/2 kernel. This result is also somewhat surprising since the performance of the squared exponential kernel is nearly identical to the Matern 3/2 kernel.

When we analyze the performance of BayesOpt in each of the three parameter domains separately; small, medium, large, see Fig. 11, we note that the medium and large domains stand out. In the former, the expected improvement acquisition with the Matern 3/2 kernel and ARD algorithm achieves a 50% performance already after 50 iterations. This is the top performing BayesOpt algorithm with the nuclear physics objective function. In the small and large domains the Matern 3/2 kernel performs best without ARD. Since BayesOpt will work best with a sensible prior, this suggests rather short correlation lengths. In the large domain, the lower confidence-bound acquisition shows a tendency to perform better than expected improvement. This is perhaps not too surprising since the parameter domain is large enough to benefit from substantial exploration. Clearly, the ARD kernels require much more data in larger spaces in order to prune irrelevant features of the objective function.

5.1 Comparison with the POUNDERs algorithm

To facilitate a benchmark and further analysis of the strengths and weaknesses of BayesOpt, we have also optimized the nuclear physics objective function using the POUNDERs (Practical Optimization Using No Derivatives for sums of Squares) optimization algorithm [25]. This is a well-tested derivative-free algorithm for optimizing non-linear least squares problems. Indeed, it has been applied successfully in basic nuclear physics research since almost a decade [26, 27]. The key feature of POUNDERs is that it assumes a least-squares structure of the objective function, i.e. that it consists of a sum of squared residuals $R_{i}(x)$ written as

[TABLE]

Tailoring an optimization algorithm to exploit this mathematical structure, i.e. that each term is a squared function of the parameters $x$, is very fruitful. A quadratic model,

[TABLE]

for an objective function $f$, at the current iterate $\mathbf{x}_{k}$, with $\mathbf{g}_{k}=\nabla f(\mathbf{x}_{k})$ and $\mathbf{H}_{k}=\nabla^{2}f(\mathbf{x}_{k})$, is a very common choice. If the corresponding derivatives are known the subproblem of minimizing $m_{k}$ can be solved quite efficiently. However, derivatives $\nabla f$ and $\nabla^{2}f$ are considered unavailable for the present problems and only a set of function values $f(\mathbf{y}^{(j)})$, and residual values $R_{i}(\mathbf{y}^{(j)})$, for some set $\mathbf{Y}=\left\{\mathbf{y}^{(1)},\mathbf{y}^{(2)},\ldots,\mathbf{y}^{(n)}\right\}$ can be accessed. POUNDERs sets up a quadratic model for each residual $i$ in Eq. 11

[TABLE]

centered around the current iterate $\mathbf{x}_{k}\in\mathcal{X}\in\mathbb{R}^{D}$. The model for each residual $i$ is defined by the parameters $c^{(i)}\in\mathbb{R}$, $\mathbf{g}^{(i)}\in\mathbb{R}^{D}$, and $\mathbf{H}^{(i)}\in\mathbb{R}^{D\times D}$. These model parameters are determined by solving an interpolation problem in $\mathbb{R}^{D}$, see Ref.  [25] and references therein. Once the model parameters are obtained, they can be used to approximate the derivatives of the objective function. In principle, the first- and second-order derivatives of $f(\mathbf{x})$ with respect to $\mathbf{x}$ are given by $\sum_{i}\nabla R_{i}(\mathbf{x})R_{i}(\mathbf{x})$ and $\sum_{i}(\nabla R_{i}(\mathbf{x})\nabla R_{i}(\mathbf{x})+R_{i}(\mathbf{x})\nabla^{2}R_{i}(\mathbf{x}))$, respectively. Consequently, POUNDERs sets up a master model $M_{k}$ for the objective function

[TABLE]

It should be noted that the master model itself does not interpolate the objective over the interpolation set $\mathbf{Y}$.

POUNDERs is a trust-region method. The master model $M_{k}$ is believed to approximate $f$ in a neighborhood $B_{k}$ of the current iterate $\mathbf{x}_{k}$, where

[TABLE]

The master model is therefore minimized over some $B_{k}$ with radius $\delta_{k}>0$. The solution $\mathbf{s}_{k}$ to ${\rm min}\left\{M_{k}(\mathbf{x}_{k}+\mathbf{s}):||\mathbf{s}||\leq\delta_{k}\right\}$ determines the next iterate $\mathbf{x}_{k+1}$ and a new trust region radius $\delta_{k+1}$ is determined based on how good the model prediction $M_{k}(\mathbf{x}_{k}+\mathbf{s}_{k})$ was, see Ref [25] for a full specification of the algorithm. We are running POUNDERs in the default mode, meaning that we only have to input $\mathbf{x}_{0}\in\mathbb{R}^{D}$, $M_{0}$, some initial trust-region radius $0<\delta_{k}<0.5$, and lower and upper bounds $l_{i},u_{i}$ of the domain $\mathcal{X}$. We also scale the problem such that the bounds correspond to a unit hypercube $[0,1]^{D}$. This scaling is not performed in BayesOpt.

In Fig. 12 we present the data profiles for POUNDERs applied to the physics objective function in the small,medium, and large domains and compare with BayesOpt. The results are only compared with the expected improvement acquisition function as it was shown in Fig. 11 to perform significantly better than the lower confidence-bound acquisition function for this optimization problem. We remind the reader that all algorithms are initiated from an identical set of 24 starting points $\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{24}\}$ for each parameter domain.

As seen clearly in Fig. 12, the choice of initial trust-region radius $\delta_{0}$ determines the performance of POUNDERs. Setting the initial radius too small ($\delta_{0}\lesssim 0.15$) hampers the POUNDERs algorithm by trapping it in a shallow local minimum. This is not an issue with BayesOpt which continues to improve as more and more function values extends the data vector $\mathcal{D}$. As we have already noted several times, the overall performance of BayesOpt depends crucially on the prior. We also note that when POUNDERs performs well, it does so within rather few function evaluations, but halts once it is trapped in a local minimum. The advantages of BayesOpt are most prominent when optimizing over the medium domain. Regardless of kernel, BayesOpt performs rather well even with few function evaluations. In the large domain, the good performance of POUNDERs clearly indicates the usefulness of encoding prior knowledge about the mathematical structure of the objective function. There is likely some large scale structure in the objective function that the BayesOpt kernel would benefit to learn about. Therefore, it would probably be advantageous to amend the BayesOpt prior with a polynomial regression term fitted to the first few evaluations of the objective function.

BayesOpt is not intended for pinpointing the exact location of an optimum. In the neighborhood of an optimum most objective functions can be well approximated by a quadratic polynomial. For this reason, POUNDERs will always outperform BayesOpt when decreasing $\tau$ in the convergence criterion given in Eq. 8. For $\tau=0.01$, i.e. a convergence criterion corresponding to a 99% reduction of the objective function, BayesOpt will only reach a performance of $d(\alpha<250)\approx 0.15$, whereas POUNDERs can approach $d(\alpha)\approx 0.35$ for an optimal choice of the initial trust-region radius $\delta_{0}$, see Fig. 13.

6 Summary and Outlook

Some of the most interesting optimization problems in nuclear physics, as well as other fields of science, typically render computationally expensive objective functions defined on multi-dimensional parameter domains. Moreover, derivatives with respect to those parameters are usually not accessible.

In this paper we explore the potential benefits of BayesOpt (Sec. 2) for efficiently exploring the parameter space of a $\chi$EFT (B) in computationally complex circumstances. A local minimum, with realistic physical properties, in this parameter domain allows for numerical simulations of atomic nuclei and therefore improves our understanding of the origin, evolution, and structure of all matter in the universe. The underlying optimization problem is therefore well known in the nuclear physics research community and several classes of numerical optimization algorithms have already been employed; derivative-based [5] as well as derivative-free approaches [27].

BayesOpt presupposes a prior on the objective function, usually in the form of a $\mathcal{GP}$. The original optimization challenge is transformed to a design problem that boils down to choosing an appropriate acquisition function to facilitate an exploration-exploitation balance. This choice is encoded in a utility function that decides where to collect training data for the $\mathcal{GP}$. Several choices of kernels and utility functions exist. Our initial studies of BayesOpt applied to a set of six test functions with parameter domains in $\mathbb{R}^{2,4,8}$ clearly demonstrate the importance of a sensible prior assumption of the objective function, see Fig. 7. From this analysis it is also clear that BayesOpt performs rather well in low-dimensional ($\mathbb{R}^{2}$ and $\mathbb{R}^{4}$) parameter domains. It turns out that the choice of acquisition function is even more important than the choice of $\mathcal{GP}$-kernel, see Fig. 8. This is something we see also when we study the data profiles of BayesOpt applied to the nuclear physics objective, see Fig. 11.

Our main findings and conclusions can be summarized as follows:

• •

In general, BayesOpt will never find a narrow minimum nor be useful for extracting the exact location of any optimum. For that to work, in anything but a trivial case, it is necessary to have detailed information about the objective function and successfully encode this prior knowledge into the algorithm. This is a situation that we typically do not have, since BayesOpt is applied to computationally expensive objective functions. Instead, BayesOpt might find use as a first stage in a hierarchical optimization protocol to identify an interesting region of the parameter domain. It might also be advantageous to design an acquisition function that is more explorative during the first iterations, and then switch to an acquisition function that exploits more than it explores.

• •

When we compare with the POUNDERs algorithm, Sec. 5.1—a derivative-free optimization algorithm that successfully incorporates the squared-sum structure of the objective function—we find that BayesOpt in ab initio nuclear physics would probably benefit from a prior with a polynomial regression term to efficiently capture the large scale structure of the objective function.

• •

We find that the choice of acquisition function is more important than the specific form of the $\mathcal{GP}$-kernel. For the present case, the expected improvement acquisition function performed slightly better than the lower confidence-bound in smaller parameter domains, while more exploration as achieved with the lower confidence-bound acquisition function, was shown to be beneficial in larger domains.

• •

The $\mathcal{GP}$-kernel can be improved with ARD tuning of the hyperparameters. However, this feature is only useful if a minimum number of iterations can be afforded. In fact, the ARD kernels requires significantly more data in larger spaces in order to prune irrelevant features of the objective function.

• •

Although we employ an objective function consisting of independent scattering cross sections, and all of them with similar and low computational cost, a multi-fidelity scenario is equally probable. Consider e.g. to calibrate an EFT model of the nuclear interaction using scattering cross sections and data on bound states in multi-nucleon systems such as isotopes of oxygen and calcium. In such scenarios, where the computational cost of solving the Schrödinger equation for bound-states of a nucleus with $A$ nucleons naively grow exponentially with $A$, it would be interesting to study the benefits of existing BayesOpt frameworks that can maximize information gathering across multiple functions with varying degrees of accuracy and computational cost. See e.g. Ref. [40] and references therein.

Acknowledgments

The research leading to these result received financial support from the BigData@Chalmers initiative. This research also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 758027), the Swedish Research Council under Grants No. 2015-00225 and 2017-04234, and the Marie Sklodowska Curie Actions, Cofund, Project INCA 600398. Computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC, Linköping.

Appendix A Test functions

In this appendix we provide the expressions for the six test functions that we used for initial analysis of BayesOpt.

• •

Ackley: one hole on a semi-flat surface with shallow local minima.

[TABLE]

• •

Deceptive: very challenging multivariate test function for which the total size of the region with local minima is $5^{D}-1$ times larger than the region with the global minimum. This function has $3^{D}-1$ local minima in $\mathbb{R}^{D}$.

[TABLE]

• •

Rastrigin: spherical function with cosine modulation to generate frequent local minima.

[TABLE]

• •

Rosenbrock: classic test function with minimum located in very shallow valley.

[TABLE]

• •

Schwefel: smooth surface with several local minima and a global minimum located far away in a corner, which in turn is far away from the second best local minimum.

[TABLE]

• •

Sphere: in many ways the simplest possible test function. It is convex and unimodal.

[TABLE]

Appendix B Physics model: $\chi$EFT for neutron-proton scattering at NNLO

In the interest of keeping this paper somewhat self-contained, we briefly review the formalism that underpins the physics model we are optimizing. To that end, we include a selected set of expressions that are representative for the full theoretical framework for computing neutron-proton ($np$) scattering cross section starting from a potential description in $\chi$EFT. Exhaustive reviews of EFT are given in [21, 22, 23]. The model evaluation $\mathcal{O}^{\rm model}(\mathbf{x})$ of an $np$ scattering cross section to be confronted with an experimentally determined value $\mathcal{O}^{\rm experiment}$ at a given on-shell momentum $k$ proceeds in three steps:

compute the momentum-space proton-neutron interaction potential. 2. 2.

use the potential to compute the quantum-mechanical scattering amplitude by solving the non-relativistic Lippmann-Schwinger equation. 3. 3.

use the amplitude to compute a model value for the scattering cross section by evaluating the spin-scattering matrix.

B.1 The $np$ interaction potential at NNLO in $\chi$EFT

The interaction potential in $\chi$EFT consists of non-polynomial terms that describe the long-range part corresponding to pion exchanges, and contact-type interactions given by polynomial expressions corresponding to a short-range part:

[TABLE]

where $\mathbf{p_{f}}$ and $\mathbf{p_{i}}$ denote the final and initial nucleon momenta in the center-of-mass system (CMS), and the special case $|\mathbf{p_{f}}\,|=|\mathbf{p_{i}}|=k$ corresponds to an on-shell momentum. The potential can be written as a systematic expansion with high-order terms being less important than low-order ones. In this work we employ a potential including terms up to NNLO in $\chi$EFT. This means that there are terms also at leading-order (LO) and next-to-leading order (NLO). In general, at higher orders there are more pion exchanges and higher-order polynomial terms in momenta that flow through the contact diagrams. At LO, the one-pion exchange potential is given by

[TABLE]

where $\mathbf{q}\equiv\mathbf{p_{f}}-\mathbf{p_{i}}$ is the momentum transfer, \mbox{\boldmath\sigma}_{1,2} and \mbox{\boldmath\tau}_{1,2} are the spin and isospin operators of nucleon 1 and 2, $g_{A}$, $f_{\pi}$, and $m_{\pi}$ denote the axial-vector coupling constant, the pion decay constant, and the pion mass, respectively. We use $f_{\pi}=92.4$ MeV and $g_{A}=1.29$ throughout this work. Higher order corrections to one-pion exchange renormalize the coupling constants $g_{A}$ and $f_{\pi}$, and the LO long-ranged part is considered parameter-free in this work. Up to NNLO, leading and sub-leading two-pion exchange enters and the long-ranged part of the interaction is given by

[TABLE]

where $w^{2}=(4m_{\pi}^{2}+q^{2})$ and the loop functions are here written as

[TABLE]

and the so-called spectral function cutoff is set to $\tilde{\Lambda}=700$ MeV. The long-ranged part contains three of 12 unknown LECs that we seek to constrain using BayesOpt, denoted $c_{1}$, $c_{3}$, and $c_{4}$ in Eq. (24). The remaining nine LECs control the short-ranged part of the NNLO potential, which can be written as a linear combination of terms polynomial in the initial and final momenta

[TABLE]

B.2 Proton-neutron scattering observables

Proton-neutron elastic scattering observables are calculated from the spin-scattering matrix $M$ [41, 42]. This is a $4\times 4$ matrix in spin-space that operates on the initial state to give the scattered part of the final state. $M$ is related to the conventional scattering matrix $S$ by $M=\frac{2\pi}{ik}(S-1)$, where $k$ is the relative momentum between the nucleons. The $S$-matrix for the scattering channel with angular momentum $J$ can be parameterized by the Stapp phase shifts $S=e^{2i\delta_{J}}$ [43]. The Stapp phase shifts are calculated from the potential $V(\mathbf{p_{f}},\mathbf{p_{i}})$ by solving the Lippmann-Schwinger equation. This equation describes quantum-mechanical scattering, and is an integral equation of the Fredholm type that can be solved as a matrix equation. In our application, and for each value of the on-shell momentum $k$, this amounts to inverting a 200-by-200 matrix followed by a matrix-vector multiplication. The matrix inversion prevents linearizing this particular EFT model in its parameters. Although, the matrix inverse is not particularly time-consuming in the present case, it should be pointed out that the complexity of the corresponding quantum-mechanical equations for describing scattering states, or bound states, of more than two nucleons typically scale exponentially with the particle number.