Quantifying Correlated Truncation Errors in Effective Field Theory

TL;DR

This paper introduces a Bayesian Gaussian process-based model to quantify correlated truncation errors in effective field theories, improving the accuracy and validation of EFT predictions, demonstrated on nucleon-nucleon scattering data.

Contribution

It develops a novel Gaussian process framework for correlated EFT truncation errors, enabling better uncertainty quantification and validation tools.

Findings

Effective Bayesian model for correlated EFT errors

Improved credible interval estimation for EFT predictions

Validated approach on nucleon-nucleon scattering data

Abstract

Effective field theories (EFTs) organize the description of complex systems into an infinite sequence of decreasing importance. Predictions are made with a finite number of terms, which induces a truncation error that is often left unquantified. We formalize the notion of EFT convergence and propose a Bayesian truncation error model for predictions that are correlated across the independent variables, e.g., energy or scattering angle. Central to our approach are Gaussian processes that encode both the naturalness and correlation structure of EFT coefficients. Our use of Gaussian processes permits efficient and accurate assessment of credible intervals, allows EFT fits to easily include correlated theory errors, and provides analytic posteriors for physical EFT-related quantities such as the expansion parameter. We demonstrate that model-checking diagnostics---applied to the case of multiple curves---are powerful tools for EFT validation. As an example, we assess a set of nucleon-nucleon scattering observables in chiral EFT. In an effort to be self contained, appendices include thorough derivations of our statistical results. Our methods are packaged in Python code, called gsum, that is available for download on GitHub.

Figures

Tables

Table 1. Table 1: Parameter estimates for the emulators in Fig.  LABEL:sub@fig:toy_coeffs_and_diagnostics_coeffs using a nugget of σ n 2 = 10 − 10 superscript subscript 𝜎 𝑛 2 superscript 10 10 \sigma_{n}^{2}=10^{-10} .

$𝜽$	Prior	MAP	True
$\mu$	$𝒩(0,0)$	0	0
$\overline{c}$	${\chi }^{--2}(1,1)$	0.97	1
$\ell$	$U(0,\infty )$	0.20	0.2
Q	$U(0,1)$	0.48	0.5

Table 2. Table 2: Notation for vectors in different spaces.

Notation	Example	Description
—	$x\in {ℝ}^d$	Input space (suppressed)
Bold	${𝐜}_n\in {ℝ}^N$	A set using inputs ${\left\{x_i\right\}}_{i=1}^N$
Arrow	${\overrightarrow{c}}_k\in {ℝ}^{n_c}$	Orders $c_0,c_1,\dots ,c_k$
Bold Italics	$𝜽$	GP parameters $\{\mu ,\overline{c},\ell \}$
—	$\mu$, $\eta$, $V$	Mean basis space
Composition	${\overrightarrow{𝐜}}_k\in {ℝ}^{n_c}\times {ℝ}^N$	Combination of above
	$𝐱\in {ℝ}^N\times {ℝ}^d$

Equations

$$y_{\text{exp}}(x)=y_{\text{th}}(x;\vec{a})+\delta y_{\text{th}}(x)+\delta y_{\text{exp}}(x),$$

$$y_n(x)=y_0(x)+\Delta y_1(x)+\cdots +\Delta y_n(x),$$

$$y_k(x)=y_{\mathrm{ref}}(x)\sum _{n=0}^k{c}_n(x)Q^n(x).$$

$$\begin{array}{rl}{y}_0(x)& \equiv y_{\mathrm{ref}}(x)c_0(x)\\ \Delta y_n(x)& \equiv y_{\mathrm{ref}}(x)c_n(x)Q^n(x).\end{array}$$

$$\delta y_k(x)=y_{\mathrm{ref}}(x)\sum _{n=k+1}^{\infty }{c}_n(x)Q^n(x).$$

$$f(x)\sim \mathrm{GP}[m(x),\kappa (x,x^{\prime })].$$

$$\mathbf{f}|\mathbf{x}\sim \mathcal{N}(\mathbf{m},K).$$

$$\kappa (x,x^{\prime };\bar{c},L)$$

$$={\bar{c}}^2{e}^{-\frac{1}{2}(x-x^{\prime }{)}^{⊺}{L}^{-1}(x-x^{\prime })},$$

$$\to {\bar{c}}^2{e}^{-(x-x^{\prime }{)}^{⊺}(x-x^{\prime })/2\ell },$$

$$\displaystyle\begin{bmatrix}\mathbf{f}_{1}\\ \mathbf{f}_{2}\end{bmatrix}\Big{|}\,{\mathbf{x}},\bm{\theta}\sim\mathcal{N}{\left(\begin{bmatrix}\mathbf{m}_{1}\\ \mathbf{m}_{2}\end{bmatrix},\begin{bmatrix}K_{11}&K_{12}\\ K_{21}&K_{22}\end{bmatrix}\right)}.$$

$${\mathbf{f}}_2|\mathbf{x},{\mathbf{f}}_1,\mathbf{\theta }\sim \mathcal{N}({\tilde{\mathbf{m}}}_2,{\tilde{K}}_{22}),$$

$${\tilde{\mathbf{m}}}_2$$

$${\tilde{K}}_{22}$$

$${\tilde{\mathbf{m}}}_2$$

$${\tilde{K}}_{22}$$

$$\mathrm{pr}(\mathbf{\theta }|\mathbf{x},\mathbf{f})\propto \mathrm{pr}(\mathbf{f}|\mathbf{x},\mathbf{\theta })\mathrm{pr}(\mathbf{\theta }).$$

$$\mathrm{pr}({\mathbf{f}}_2|\mathbf{x},{\mathbf{f}}_1)=\int \mathrm{pr}({\mathbf{f}}_2|\mathbf{x},{\mathbf{f}}_1,\mathbf{\theta })\mathrm{pr}(\mathbf{\theta }|{\mathbf{x}}_1,{\mathbf{f}}_1)\mathrm{d}\mathbf{\theta }.$$

$$c_n(x)|\mathbf{\theta }$$

$$\mathbf{\theta }$$

$$AX+BY\sim \mathcal{N}(A{\mu }_1+B{\mu }_2,A{K}_1{A}^{⊺}+B{K}_2{B}^{⊺}).$$

$$\delta y_k(x)|\mathbf{\theta },Q\sim \mathrm{GP}[m_{\delta k}(x),{\bar{c}}^2{R}_{\delta k}(x,x^{\prime };\ell )],$$

$$m_{\delta k}(x)$$

$$R_{\delta k}(x,x^{\prime };\ell )$$

$$\mu ,{\bar{c}}^2\sim \mathrm{N}{\chi }^{--2}({\eta }_0,V_0,{\nu }_0,{\tau }_0^2).$$

$$\mu |{\bar{c}}^2$$

$${\bar{c}}^2$$

$$f(z;\nu ,{\tau }^2)=\frac{(\nu {\tau }^2/2{)}^{\nu /2}}{\Gamma (\nu /2)z^{1+\nu /2}}\exp (-\frac{\nu {\tau }^2}{2z}).$$

$$\mu ,{\bar{c}}^2|{\vec{\mathbf{y}}}_k,\ell ,Q\sim \mathrm{N}{\chi }^{--2}(\eta ,V,\nu ,{\tau }^2)$$

$$E[\mu |{\vec{\mathbf{y}}}_k,\ell ,Q]$$

$$E[{\bar{c}}^2|{\vec{\mathbf{y}}}_k,\ell ,Q]$$

$$\delta y_k$$

$${\mathrm{TP}}_{\nu }\left\{b_{\delta k}(x)\eta ,{\tau }^2\left[R_{\delta k}(x,x^{\prime };\ell )+b_{\delta k}(x)V{b}_{\delta k}(x^{\prime })\right]\right\}.$$

$$m_{\text{th}}(x)$$

$${\Sigma }_{\text{th}}(x,x^{\prime };\ell )$$

Full text

I Introduction

The power counting of an effective field theory (EFT) mandates how to organize an infinite number of operators into a sequence of decreasing importance. The precision of EFT predictions is then governed by the uncertainty in the following quantities: (1) the fit of the parameters, or low energy constants (LECs), of the EFT, (2) the error due to truncation of the infinite EFT series, and (3) other approximations made in prediction. For chiral EFT calculations, recent advances in precision motivate a rigorous accounting of all these theoretical uncertainties, but such an accounting is desirable for EFTs across all domains. This work will focus on the uncertainty due to truncation, but our results have consequences for how LECs should be fit Wesolowski et al. (2019).

In previous works, we proposed a pointwise Bayesian statistical model to estimate EFT truncation errors for predicted observables $y$ Furnstahl et al. (2015); Melendez et al. (2017) (see also Ref. Cacciari and Houdeau (2011)). This model formalizes the notion of convergence in $y$, which allows one to credibly assert which observations are consistent with theory. One can incorporate expert knowledge on the convergence pattern (inherited from the EFT power counting) into prior distributions, and subsequently update these beliefs given order-by-order predictions $\{y_{n}\}$ Cacciari and Houdeau (2011). If the observable is itself a function $y({x})$, e.g., a cross section at a range of energies, Refs. Furnstahl et al. (2015); Melendez et al. (2017) simply compute its posterior in a pointwise manner. This pointwise model is tractable, but it is flawed since it ignores information about $y({x})$ as a whole, that is, the correlations of $y$ at nearby ${x}$.

Here we extend the pointwise model to functions $y({x})$, encoding the idea of curvewise convergence for observables via Gaussian processes (GPs). GPs are powerful tools for both regression and classification, and have become popular in statistics, physics, applied mathematics, machine learning, and geostatistics Sacks et al. (1989); Cressie (1992); Rasmussen and Williams (2006). Their popularity is due in part to their modeling flexibility and the mathematical convenience of Gaussian distributions. Despite the leap from points to curves, our algorithm remains analytic due to a convenient yet flexible choice of priors. The GP parameters are interpretable from an EFT convergence standpoint, and can be easily calibrated against known order-by-order predictions.

We develop the intuition behind our convergence model and provide a primer on GPs (Sec. II), and then explore ways to assess whether our assumptions are failing (Sec. III). To display the procedure in action, we show success and failure using toy EFT predictions (Sec. IV) and assess nucleon-nucleon ($NN$) scattering predictions in chiral EFT (Sec. V). We then discuss future prospects and conclude the main text (Sec. VI) before providing an extended discussion on the statistical derivations relevant to the pointwise and GP models (Appendix A).

II The Model

II.1 Setup

Predictions of experimental quantities can be inaccurate due to shortcomings of the theoretical framework and noisy measurements. This work focuses on building a model discrepancy term for EFT truncation errors. The exact manner by which specific errors propagate can be complicated, but here we assume a simple additive model of theoretical discrepancy,

[TABLE]

where ${x}\in\mathbb{R}^{d}$ is the independent variable,111The differential cross section, as a function of both energy and scattering angle, is an example where $d=2$. $\vec{a}$ are the fitted parameters of the theory prediction $y_{\textup{th}}$ (here, the LECs of an EFT222Note that the LECs $\vec{a}$ are not the same as the observable coefficients $c_{n}$ introduced in Eq. (3). These parameters $\vec{a}$ do not play an explicit role in our treatment of truncation errors, although they do appear via a brief discussion of EFT parameter fitting in Sec. II.4.), and $y_{\textup{exp}}$ are experimental measurements. The uncertainties associated with the theory and experiment are given by $\delta y_{\textup{th}}$ and $\delta y_{\textup{exp}}$. Of course, if the exact values of $\delta y_{\textup{th}}$ and $\delta y_{\textup{exp}}$ were known, then there would be no uncertainty. Instead, $\delta y_{\textup{th}}$ and $\delta y_{\textup{exp}}$ are treated as random variables, which makes Eq. (1) true in a distributional sense. The simple relationship (1) directly impacts both the fitting protocol and error propagation to observables. For example, if one assumes that $\delta y_{\textup{th}}=0$ and $\delta y_{\textup{exp}}$ is normally distributed, this leads to a least squares log likelihood function for $\vec{a}$ given $y_{\textup{exp}}$ Wesolowski et al. (2019). Here we use general properties of EFTs to develop a Bayesian model for $\delta y_{\textup{th}}({x})$ with a non-trivial correlation structure depending on ${x}$. Unless otherwise stated, we assume that the EFT has already been fit to experimental data and return to discuss the fitting procedure in Sec. II.4.

Denote an observable $y$ calculated with an $n$th-order EFT ($\mathrm{N}^{n}\mathrm{LO}$) as $y_{n}$. Suppose that the EFT expansion has only been computed up to $k$th order. Though $y_{k}$ may work well as an approximation to reality, we seek a probability distribution for the full summation333 We are aware that EFTs are often asymptotic in nature. We seek the best possible estimate up to where the sum starts diverging. Our assumption is that the divergence occurs sufficiently beyond the $k$th order that the truncation model error estimate up to order $\infty$ is a good approximation. $y\equiv y_{\infty}$ based on the convergence pattern of all known predictions $y_{0},y_{1},\dots,y_{k}$. That is, $y_{\textup{th}}=y_{k}$ and we wish to quantify $\delta y_{\textup{th}}=\delta y_{k}$. This hierarchy of predictions could instead be written as a leading-order calculation $y_{0}$ and a set of higher-order corrections $\Delta y_{n}$, with $n\in\{1,2,\cdots,k\}$. Given this change of variables, the $n\leq k$ order prediction is

[TABLE]

where each term in the series is known.

Without further knowledge of the convergence pattern in Eq. (2), it is difficult to make progress in approximating $y$. To proceed, we employ prior information about the construction of EFTs, that is, if the EFT is working as advertised, each correction should be roughly suppressed by the dimensionless expansion parameter $Q$—in accordance with the EFT power counting. Here we assume for simplicity the following power counting: the $\mathrm{N}^{n}\mathrm{LO}$ correction $\Delta y_{n}$ is suppressed by $Q^{n}$, though other prescriptions, e.g., $Q^{2n}$, can be easily incorporated into this framework. That is, up to dimensionful scales, one would generally expect $y_{0}$ is $\order{Q^{0}}$, and $\Delta y_{n}$ is $\order{Q^{n}}$.

Inspired by Eq. (2) and the power counting of the EFT, Ref. Furnstahl et al. (2015) proposed the following factorization

[TABLE]

Here, $y_{\mathrm{ref}}$ is a dimensionful quantity that sets the scale of variation, while the $c_{n}$ are dimensionless observable coefficients. Since all scales have been factored into $y_{\mathrm{ref}}$ and $Q$, the $c_{n}$ should be natural, or order 1, assuming they have not been fine tuned. Combinatorial factors, if known, should also appear in Eq. (3), else the naturalness assumption may need to be modified to account for them.

The road-map from continuous observable quantities $y_{0},\dots,y_{k}$ to coefficients is shown visually for toy EFT predictions in Fig. 1. Predictions for the four lowest orders of the EFT are plotted in Fig. 1(a), and show a steady convergence towards the final result. Suppose that the $\mathrm{N}^{3}\mathrm{LO}$ EFT is the state of the art; thus our task is to estimate the uncertainty in $y_{3}$. The change of variables described by Eq. (2) is illustrated in Fig. 1(b). Note that there is a clear hierarchy to the corrections: they tend to become smaller as order increases. Assume that $y_{\mathrm{ref}}$ is known, say, from dimensional analysis, and $Q$ is known from theory arguments. The coefficients of the expansion can then be extracted using Eq. (3), as shown in Fig. 1(c). The hierarchy has vanished, and we are left with a set of curves that appear randomly drawn from an underlying process. Reference Melendez et al. (2017) displays $c_{n}$ for a collection of $NN$ scattering observables in chiral EFT, many of which exhibit these properties.

Before moving on, we would like to correct some common misconceptions about Eq. (3). We have made no assumptions about the structure of the expansion thus far. In fact, Eq. (3) is not a formal expansion of $y$ in powers of $Q$. Rather, we have only noted that if we assume values for $y_{\mathrm{ref}}$ and $Q$ along with the factorization of Eq. (3), then there is a one-to-one correspondence between the predictions $\{y_{0},\dots,y_{k}\}$ and the coefficients $c_{n}$ given by

[TABLE]

The coefficients are not the LECs of the EFT. Because we are making predictions that take the LECs as input, the $c_{n}$ are potentially complicated functions of the LECs, $Q$, and other variables. Therefore, we do not have a general proof that the naturalness of the LECs propagates to the $c_{n}$; rather, it must be verified in each application. If some of the $c_{n}$ do not conform to our assumptions, they may need to be left out of the analysis.

Our manipulations have thus far have simply defined our terminology, but the arguments above lead to a physics-based uncertainty model that can be tested against reality. By a logical extension of Eq. (3), the truncation error $\delta y_{k}$ would consist of all terms not included in the sum, i.e.,

[TABLE]

We can now describe the details of the EFT convergence model. By induction, we assume that the properties of the unobserved $c_{n}$ for $n>k$ are the same as the $c_{n}$ for $n\leq k$. Specifically, we assume that the $c_{n}$ are independent and identically distributed (i.i.d.) random curves. Moreover, we have some prior knowledge of the curves’ properties: they should be smooth and naturally sized, i.e., have a standard deviation of order 1. Data from the lower-order predictions can then be used to refine the estimates for the sizes and shapes of all subsequent $c_{n}$ and, via Eq. (5), estimate the truncation error. This is shown visually in the Bayesian network of Fig. 2.

Because the coefficients are treated as random curves, the distribution from which they are drawn is known as a stochastic process. A stochastic process is essentially an infinite dimensional generalization of a probability distribution. Processes associate a random variable with each point in a continuous domain, hence their infinite dimensionality, and correlations between the random variables at different points can be built in. Through the appropriate assumption on this correlation structure, one can enforce desirable qualities in random draws from the process, such as continuity of the drawn functions and their derivatives. GPs are particularly useful and will be employed throughout this paper.

II.2 A Brief Introduction to Gaussian Processes

II.2.1 Definitions

GPs are a popular tool for nonparametric regression due to their flexible nature and tractable analytic forms. They are used in fields such as statistics, machine learning, and geostatistics, where GP regression is known as kriging Sacks et al. (1989); Cressie (1992); Rasmussen and Williams (2006). For appropriate priors, GPs arise as the limit of a specific type of neural network Neal (1996) and are equivalent to common interpolants, e.g., splines MacKay (1998). Here we will review the most common applications of GPs: interpolation and regression. Then we discuss the calibration of the GP parameters, which is a crucial step in our truncation error model. For more in-depth discussions of GPs, see Refs. Rasmussen and Williams (2006); MacKay (1998, 2003).

The defining quality of a GP is Rasmussen and Williams (2006)

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution.

A GP is specified by a mean function $m({x})$ and a positive semi-definite covariance function, or kernel, $\kappa({x},{x}^{\prime})$, where ${x}\in\mathbb{R}^{d}$. Assuming both the mean and covariance functions are known, a GP $f({x})$ is denoted by444The $z\sim\cdots$ notation is a common shorthand in statistical literature for “$z$ is distributed as.” Some authors use $\operatorname{pr}(z)=\cdots$ as well Gelman et al. (2013). Also “$z\,|\,I$” is read as “$z$ given $I$.”

[TABLE]

We do not operate on the infinite-dimensional object in practice; rather, the definition of GPs allows us to compute using a finite number of points. Let ${\mathbf{x}}=\{{x}_{i}\}_{i=1}^{N}$ be a set of $N$ input points and $\mathbf{f}=\{f({x}_{i})\}_{i=1}^{N}$ be the corresponding set of function values.555To specify how the GP enters in our truncation model, we must work in multiple vector spaces. The bold notation for vectors signals an element of $\mathbb{R}^{N}$ while we use arrows to indicate a vector in the space of EFT orders $\mathbb{R}^{n_{c}}$. Bold with arrows is used for elements of both spaces. The $d$ dimensional input space is not explicitly indicated. This notation is summarized in Table 2 in Appendix A. Henceforth, we refer to coefficient values computed at particular input points as “data”, because these are the points used to train the GP emulator.

Define $\mathbf{m}=m({\mathbf{x}})\in\mathbb{R}^{N}$ and $K=\kappa({\mathbf{x}},{\mathbf{x}})\in\mathbb{R}^{N\times N}$. $\mathbf{f}$ is then distributed as the multivariate normal

[TABLE]

The fact that Eq. (6) implies Eq. (7) is the definition of a GP in a mathematical form. Often the explicit conditioning on input points ${\mathbf{x}}$ from $\mathbf{f}\,|\,{\mathbf{x}}$ is dropped for notational simplicity. The $N$-dimensional Gaussian distribution of Eq. (7) can then easily be used to draw samples of $\mathbf{f}$ using the Cholesky decomposition of $K$.

The mean function is an a priori “best guess” of $f$, and should capture overall trends in the data; the GP then soaks up deviations from this mean in accordance with the chosen kernel Gelman et al. (2013). Basis function models for $m({x})$ are convenient choices due to their analytic tractability and flexibility Gelman et al. (2013). But the effects of mean functions can also be shuffled into the kernel (see Appendix A) or simply subtracted from the data in preprocessing, hence many authors set the mean to zero. The convergence model proposed here uses a constant (possibly unknown) mean $\mu$, which we leave in the mean function slot for clarity.

One of the most popular covariance functions is the squared exponential (a.k.a. radial basis function or Gaussian). It is defined by

[TABLE]

where $r({x},{x}^{\prime};L)$ is the correlation function and $L$ is a $d\times d$ diagonal matrix of correlation length parameters, allowing one to model cases where the correlation length varies by direction Rasmussen and Williams (2006). For simplicity our examples are in one dimension ($d=1$), so $L$ reduces to a single correlation length $\ell$. Functions drawn from a GP with a squared exponential kernel are infinitely differentiable. The marginal variance $\bar{c}^{2}$ controls the width of variation about the mean, while the length scale $\ell$ controls the correlation between points in the input domain, see Fig. 3. Some have argued that the squared exponential is too smooth, and thus not applicable in many situations Stein (1999). Thus, an alternative kernel is the Matérn, which provides an extra parameter $\nu$ and results in realizations that are $\lceil\nu\rceil-1$ times differentiable, reproducing the squared exponential in the limit $\nu\to\infty$.

The squared exponential kernel, along with the Matérn, are stationary, which means that they only depend on the absolute difference in input space $|{x}-{x}^{\prime}|$. Draws from a stationary process would, on average, behave similarly across the domain of interest. Stationarity is a strong assumption about the physics of the system, which can be advantageous in data-sparse regimes that could benefit from extra constraints. But its validity should be checked. We introduce diagnostics that can point to issues with the assumption of stationarity, along with the choices of mean and covariance, in Sec. III.

We assume the functional form of the mean and covariance functions (e.g., squared exponential) are known. Their parameters,666We choose to call $\mu$, $\bar{c}$, etc., the parameters of the GP to distinguish them from the hyperparameters we will introduce later. They are more regularly referred to as hyperparameters. on the other hand, can be updated based on observed data. We denote the set of all parameters in both $m$ and $\kappa$ (e.g., $\mu$, $\bar{c}^{2}$, $\ell$) as $\bm{\theta}$. We are careful to write $\cdot\,|\,\bm{\theta}$ when $\bm{\theta}$ is given (i.e., conditioned upon) and omit $\bm{\theta}$ otherwise.

II.2.2 Interpolation and Regression

To use GPs for interpolation or regression, we must be able to incorporate (noisy) observations into our predictions. Interpolation, in this case, simply refers to the procedure of finding a curve that passes directly through known data, regardless of whether the predictions are made within or outside the support of the data. Consider the partition ${\mathbf{x}}={\begin{bmatrix}{\mathbf{x}}_{1}&{\mathbf{x}}_{2}\end{bmatrix}}^{\intercal}$ and $\mathbf{f}=f({\mathbf{x}})={\begin{bmatrix}\mathbf{f}_{1}&\mathbf{f}_{2}\end{bmatrix}}^{\intercal}$ into $N_{1}$ training and $N_{2}$ test points, with corresponding partitioning of $\mathbf{m}$ and $K$. We would like to predict the unseen function values $\mathbf{f}_{2}$ given observed function values $\mathbf{f}_{1}$. By definition of GPs, their $N_{1}+N_{2}$ dimensional joint distribution is Gaussian:

[TABLE]

Then by standard Gaussian identities Rasmussen and Williams (2006),

[TABLE]

where

[TABLE]

Since we assume that $K$ can factorize into a marginal variance $\bar{c}^{2}$ and a correlation matrix $R$, i.e., $K=\bar{c}^{2}R$, this can be written as

[TABLE]

This work considers a hierarchy of noiseless computer simulations as its data, so interpolation is the primary fitting method and an example is shown in Fig. 4. Nevertheless, interpolation may require a small amount of noise called a nugget for the numerical stability of the matrix inversion Ababou et al. (1994); Andrianakis and Challenor (2012); Mohammadi et al. (2016); Pepelyshev (2010); Ranjan et al. (2011).

To make predictions of $\mathbf{f}_{2}$ assuming Gaussian white noise in the observed data $\mathbf{f}_{1}$, one need only make the replacement $K_{11}\to K_{11}+\sigma_{n}^{2}\mathds{I}_{N_{1}}$ in Eqs. (12) and (13), leaving $K_{21}$, etc., unchanged. This is what we refer to as regression. The identity matrix is denoted by $\mathds{I}$ with the rank as subscript, and here $\sigma_{n}^{2}$ is the nugget used to regularize the matrix inversion. To make predictions that also include the effects of noise, one must also make the replacement $K_{22}\to K_{22}+\sigma_{n}^{2}\mathds{I}_{N_{2}}$. An example of regression with noise is shown in Fig. 4. We include the nugget via $\bar{c}^{2}(R+\sigma_{n}^{2}\mathds{I}_{N})$ to integrate $\bar{c}^{2}$ out analytically.

II.2.3 Calibration

Until now, we have assumed that the GP parameters $\bm{\theta}$ are known and correct, but that is rarely ever the case in practice. We often want to find the $\bm{\theta}$ that produce the best fit in interpolation/regression, or $\bm{\theta}$ may be of interest in its own right. The procedure of finding the posterior for $\bm{\theta}$ values is known as calibration, or simply Bayesian inference for the model parameters Santner et al. (2018).

Given data $\mathbf{f}$ at input points ${\mathbf{x}}$, the posterior for $\bm{\theta}$ is

[TABLE]

If $f$ is modeled as a GP, then Eq. (16) is, up to a normalization constant, Eq. (7) multiplied by the priors. The posterior can be sampled via Markov chain Monte Carlo (MCMC) methods, or optimized using gradient descent to find its maximum a posteriori (MAP) value $\bm{\theta}_{\mathrm{MAP}}$. When a single value for $\bm{\theta}$ is used as a “best guess,” e.g., $\bm{\theta}_{\mathrm{MAP}}$, this is called a point estimate.

For certain cases, one can find analytic posteriors for parameters in $\bm{\theta}$ by making use of conjugate priors, whose posteriors have the same functional form as the prior Gelman et al. (2013). For our applications only a few features of the prior information are important; we are generally insensitive to details of the priors’ functional form. Our strategy is then to make use of conjugate priors where they are flexible enough to capture distributions’ relevant features, and find MAP values for the remaining parameters.

Our conjugate priors on the GP parameters $\mu$ and $\bar{c}^{2}$ are then characterized by hyperparameters that quantify the mean and variance for the prior distributions of both, as well as the heaviness of the $\bar{c}^{2}$ distribution’s tail. Our model is also easily built into existing MCMC samplers, such as PyMC3 Salvatier et al. (2016), where any functional form for the priors can be used.

To make predictions unconditional on $\bm{\theta}$, which take all plausible values into account, one must compute

[TABLE]

This can be approximated by replacing the integral over $\bm{\theta}$ by evaluation of the integrand at $\bm{\theta}_{\mathrm{MAP}}$, computed numerically with Monte Carlo methods, or evaluated analytically with the help of conjugacy. For example, Appendix A shows how $\mu$ and $\bar{c}^{2}$ are integrated out of Eq. (17) analytically via conjugacy.

This brief introduction to GPs is enough to set the stage for our model of truncation errors.

II.3 The GP Model of EFT Convergence

Recall from Sec. II.1 that after making the transformation from prediction functions $\{y_{n}\}$ to coefficient functions $\{c_{n}\}$, the coefficients appeared as naturally sized random curves. Here we will specify exactly what we mean by random and discuss how we can extract meaningful insight from this assumption.

We formalize our EFT convergence assumptions via Eqs. (3) and (5), where the coefficients $c_{n}$ are independent draws from a single underlying GP. We assume a constant mean function $m({x})=\mu$ for simplicity and a kernel $\kappa({x},{x}^{\prime};\bar{c},\ell)=\bar{c}^{2}r({x},{x}^{\prime};\ell)$ that need not be stationary. Thus,

[TABLE]

Figure 3 shows random draws from such a process and a comparison to the model described in Refs. Furnstahl et al. (2015); Melendez et al. (2017).

The parametrization of Eq. (18) can be generalized as appropriate. Non-constant mean functions can be incorporated by promoting $\mu$ to a vector of regression coefficients multiplying basis functions. For the model of the $c_{n}$’s that we use here this basis would contain only one, constant, function. But if one wants to build specific structure into the mean function for $y_{k}$ and $\delta y_{k}$, then the notion of a basis is useful for discussing the distribution of these quantities. We discuss that more general case, adopting vector and matrix notation, in Appendix A. The factorization of the kernel into a marginal variance $\bar{c}^{2}$ and a correlation function $r$ that depends on a correlation length $\ell$ may seem restrictive. On the contrary, this kernel can incorporate a wide range of functional structures since $r$ can be non-stationary in general and $\ell$ can stand in for any set of correlation parameters. Thus here we make $\mu$, $\bar{c}^{2}$, and $\ell$ the explicit parameters of our GP model: this allows us to show how they can be calibrated or integrated out of this system.

The distribution of $\delta y_{k}$, and thus the full prediction $y=y_{k}+\delta y_{k}$, can be derived using the remarkable property of Gaussian random variables: they are closed under addition and matrix multiplication. Let $X\sim\mathcal{N}(\mu_{1},K_{1})$ and $Y\sim\mathcal{N}(\mu_{2},K_{2})$ be length-$N$ independent random variables and let $A,B$ be known $M\times N$ matrices. Then

[TABLE]

Equation (5), which defines the truncation error $\delta y_{k}$, is a geometric sum over normally distributed variables. Using Eqs. (18) and (20), one can readily find the prior

[TABLE]

where

[TABLE]

We have defined the basis function $b_{\delta k}({x})$ as it will recur below. Although $R_{\delta k}$ may look like a correlation matrix, the inclusion of the $y_{\mathrm{ref}}$ and $Q$ factors mean that $R_{\delta k}({x},{x};\ell)\neq 1$, in general. The marginal variance of Eq. (21) is thus, in general, ${x}$ dependent and equal to $\bar{c}^{2}R_{\delta k}({x},{x};\ell)$.

Equations (21)–(23) define the functional form of the truncation errors, but how should one estimate their parameters in a statistically rigorous manner? The simplest (and least rigorous) way is to assume some reasonable values for $\bm{\theta}$ and $Q$ based on expert knowledge of the system. For example, $\mu=0$ and $\bar{c}\approx 1$ is a good place to start if one can assume that the coefficients are naturally sized and equally likely to be positive or negative. Then $\ell$ and $Q$ could be based on arguments specific to the observable and EFT under consideration. Again, $y_{\mathrm{ref}}$ is assumed to be known in this work, possibly based on dimensional analysis for the observable. The next best way would be to combine expert knowledge of the system with the known low-order predictions to compute point estimates of $\bm{\theta}$ and $Q$. (Remember, information only propagates from the order-by-order predictions to $\delta y_{k}$ via $\bm{\theta}$ and $Q$, as shown in Fig. 2.) Finally, the most rigorous Bayesian method is to compute the full posterior for $\bm{\theta}$ and $Q$ and marginalize (integrate) over all of their possible values. Depending on the specific application of this model and the needs of the user, different degrees of rigor may be warranted.

Assuming that we want to combine expert knowledge of the system with the order-by-order EFT predictions to estimate $\delta y_{k}$, we must start by formalizing our knowledge as Bayesian priors. We have checked the prior dependence of our pointwise model in previous works and found that truncation error estimates are fairly insensitive to their exact choice Furnstahl et al. (2015); Melendez et al. (2017). Therefore, we choose to work with conjugate priors, which yield intuitive analytic posteriors and point estimates for the GP parameters $\bm{\theta}$. We place a (scaled) normal-inverse-$\chi^{2}$ prior777 A more common and mathematically equivalent choice is the normal-inverse-gamma prior. We find the normal-inverse-$\chi^{2}$ prior more intuitive to work with in this context and thus promote its use. on $\mu$ and $\bar{c}^{2}$,

[TABLE]

This joint prior factorizes as $\operatorname{pr}(\mu,\bar{c}^{2})=\operatorname{pr}(\mu\,|\,\bar{c}^{2})\operatorname{pr}(\bar{c}^{2})$, where

[TABLE]

The inverse $\chi^{2}$ density is given by

[TABLE]

See Fig. 5 for a graphical representation of Eqs. (25)–(26). The meaning of the hyperparameters $\eta_{0},V_{0},\nu_{0},\tau_{0}^{2}$ becomes clear from Eqs. (25)–(26).

• •

The prior mean $\eta_{0}$: Our best guess for $\mu$ before seeing the data.

• •

The prior dispersion888We have not found a standard name for this quantity in the literature. The dispersion, a generic term describing the variability or spread of a distribution, is as good as any. $V_{0}$: When combined with $\bar{c}^{2}$, this is our estimated uncertainty in $\mu$. It makes sense that if the spread of the $c_{n}$ is large, then we would be less certain about their mean.

• •

The prior degrees-of-freedom $\nu_{0}$: This quantity essentially describes our uncertainty in $\bar{c}^{2}$. For $\nu=0$, Eq. (27) reduces to the scale invariant Jeffreys prior used in Ref. Melendez et al. (2017), whereas for $\nu\to\infty$, it becomes sharply peaked.

• •

The prior scale $\tau_{0}^{2}$: For large $\nu$, this is our best guess for $\bar{c}^{2}$. The mode and mean of Eq. (27) are given by $\nu\tau^{2}/(\nu\pm 2)$, respectively. The mean is only defined for $\nu>2$.

There are no conjugate priors for $\ell$ or $Q$. Reasonable choices could be a strictly positive log-normal prior for $\ell$ and a beta prior for $Q$, since it is restricted to be between zero and one.

Once priors have been chosen for $\bm{\theta}$ and $Q$, the rest of the process is algorithmic. Given order-by-order data $\vec{\mathbf{y}}_{k}\equiv\{\mathbf{y}_{i}\}_{i\leq k}$, one can derive analytic posteriors for all of $\bm{\theta}$ and $Q$. See Appendix A for details on how this is done. From now on we assume that $\ell$ and $Q$ are known a priori or point estimates are obtained from their posteriors. Analytic results for the truncation error posterior follow. Formulae in which $\ell$ and $Q$ are marginalized over are more complicated and not analytic. Point estimates of $\ell$ and $Q$ are often good approximations in our application, making numerical integration a needless complication. MCMC sampling can be used in cases where this approximation is not adequate.

Conjugacy ensures that the posterior for $\mu,\bar{c}^{2}$ has the same functional form as the prior. That is,

[TABLE]

for updated hyperparameters $\eta$, $V$, $\nu$, and $\tau^{2}$ given by Eqs. (74)–(77). The path forward depends on whether point estimates of $\mu$ and $\bar{c}^{2}$ are sufficient. If so, then the mean values $\operatorname{\mathbb{E}}[\cdot]$ of $\mu$ and $\bar{c}^{2}$, given by

[TABLE]

can be used as point estimates in Eq. (21). It is interesting to note that these estimates are Bayesian analogs to standard frequentist estimators for the mean and variance. If one integrates over all possible values of $\mu$ and $\bar{c}^{2}$, then the posterior is a Student-$t$ process (TP) Shah et al. (2014),

[TABLE]

Like a GP, any finite collection of random variables from a TP has a joint Student-$t$ distribution. Equation (31) is centered at our best guess for $b_{\delta k}({x})\mu$ but also includes the uncertainty in $\mu$ via $b_{\delta k}({x})Vb_{\delta k}({x}^{\prime})$ and the uncertainty in $\bar{c}^{2}$ via the heavy tails of the $t$ distribution (which are heavier for smaller $\nu$). The $\kappa({x},{x}^{\prime})$ in $\operatorname{\mathcal{TP}}_{\nu}[m({x}),\kappa({x},{x}^{\prime})]$ is not the covariance function in our notation, but instead a scale function. The covariance between the response at $x$ and $x^{\prime}$ for a TP is given by $\nu\kappa({x},{x}^{\prime})/(\nu-2)$ and is only defined for $\nu>2$. Given this fact, note the relationship between the point estimate in Eq. (30) and the covariance in Eq. (31).

II.4 Application Types

We have shown how to arrive at a physically motivated distribution for the truncation error, which we will now use as a springboard to discuss how our EFT convergence model is applicable to two situations that can arise when fitting and predicting.

For certain physical systems, predictions may be expensive to compute. This can affect the fitting of the LECs, where predictions must be made repeatedly to find a best fit (or posterior), but can also affect predictions with fixed LECs for particularly expensive systems. For example, fitting chiral EFT beyond the two-body sector becomes computationally expensive. 2. 2.

Certain observables may have symmetry constraints, which provides information about contributions from higher orders in the EFT. For example, certain spin observables in the $NN$ sector of chiral EFT must vanish at $\theta=0^{\circ}$ or $180^{\circ}$. Thus, the distribution for the truncation error should reflect this fact.

Here we discuss how our correlated convergence model can solve both of these problems, and how to use this model when fitting an EFT. We assume that point estimates from Eqs. (29) and (30) are used so that everything remains Gaussian, and leave the TP case for Appendix A.

Inexpensive Predictions: If predictions are inexpensive, then $y_{k}({x})$ can be computed at all ${x}$ of interest. In this case, we do not need to interpolate $y_{k}({x})$ and can instead focus on estimating the truncation error, i.e., estimating $\bm{\theta}$ and $Q$. In this data-rich environment, we must be careful when learning the point estimates for these quantities, since too much data within one correlation length $\ell$ can cause ill-conditioned matrices and thus poor point estimates Ababou et al. (1994); Andrianakis and Challenor (2012); Mohammadi et al. (2016); Pepelyshev (2010); Ranjan et al. (2011). We suggest simply choosing a representative subset of data for learning $\bm{\theta}$ and $Q$ and using a small nugget to regularize the matrix inversion. Then $\delta y_{k}$ can be produced at all points of interest without issue using Eq. (21). The distribution for the full observable $y=y_{k}+\delta y_{k}$ has mean and covariance function given by

[TABLE]

An example of the posteriors for $y({x})$ assuming different max orders $k$ is given in Fig. LABEL:sub@fig:gp_applications_inexpensive.

Expensive Predictions: If predictions are expensive, then one may want an estimate of $y_{k}({x})$ at some ${x}$ that was not explicitly computed with the EFT. Thus, we now have both interpolation error and truncation error to worry about. Fortunately, GPs are made for interpolation and are very fast. In this case, we need a distribution for $y_{k}({x})$ as opposed to $\delta y_{k}({x})$. Following Eqs. (3), (18), and (20), one can write

[TABLE]

where

[TABLE]

Now it is a simple matter to estimate $\bm{\theta}$ and $Q$ as before, and then apply Eqs. (11)–(13) to Eq. (34). Because the conditional distribution of Eq. (11) is still Gaussian, one can again use the fact that GPs are closed under addition to compute $y=y_{k}+\delta y_{k}$, which has mean and covariance functions

[TABLE]

[Remember that the tildes above $m$ and $R$ refer to the conditioning shown in Eqs. (11)–(13).] See Fig. LABEL:sub@fig:gp_applications_expensive for an example when only a few points are computed at each order and used to estimate the full prediction with combined interpolation errors and truncation errors.

Symmetry Constraints: Above we discussed how to constrain $y_{k}({x})$ based on known data, but sometimes we could constrain $\delta y_{k}({x})$ based on, e.g., symmetry arguments. That is, we may know that higher-order coefficients are zero at a particular ${x}$. In this case, we need only constrain $\delta y_{k}({x})$ via Eqs. (11)–(13) according to our knowledge of the system. This could apply regardless of whether $y_{k}({x})$ is expensive or inexpensive. Again, since Eq. (11) is still Gaussian, $y=y_{k}+\delta y_{k}$ is a GP. One need only make the replacements $m_{\delta k}\to\tilde{m}_{\delta k}$ and $R_{\delta k}\to\tilde{R}_{\delta k}$ in Eqs. (32)–(33) and (37)–(38). Figure LABEL:sub@fig:gp_applications_constraint shows an example of how information about symmetries can be applied to an inexpensive system.

EFT Fitting: Our model permits—for the first time—a rigorous inclusion of all higher-order terms and their correlations within the fitting procedure. Interpolation of expensive predictions fits seamlessly into this framework, which allows more experimental data to be used than might otherwise be possible. The posterior for the LECs $\vec{a}_{k}$ of the $k$th order EFT, given data for different observables $\mathcal{D}=\{\mathbf{y}_{\textup{exp}}\}_{i}$, is given by Bayes’ theorem

[TABLE]

where $I_{i}$ is the prior information for the $i$th observable, and we have assumed that different observables are conditionally independent of one another given $\vec{a}_{k}$.

Consider one specific observable $y$. Prior information about $y$ could include $\bm{\theta}$ and $Q$ (from previous lower-order fits of the EFT, for example) and also its experimental uncertainty $\Sigma_{\textup{exp}}$. The values of $\vec{a}_{k}$, $\bm{\theta}$, and $Q$ are sufficient to construct the distribution for $y({x};\vec{a}_{k})=y_{k}({x};\vec{a}_{k})+\delta y_{k}({x})$, which is a GP with mean function $m_{\textup{th}}({x};\vec{a}_{k},\bm{\theta},Q)$ and covariance function $\Sigma_{\textup{th}}({x},{x}^{\prime};\vec{a}_{k},\bm{\theta},Q)$. The exact form of these functions is dependent on whether the predictions are interpolated and whether the truncation error is constrained, as discussed above. Once the distribution for $y({x};\vec{a}_{k})$ is determined, then the likelihood for this observable is given by

[TABLE]

This is an extension of the likelihood proposed in Ref. Wesolowski et al. (2019), whose “uncorrelated” and “fully correlated” models correspond to $\ell=0$ and $\ell=\infty$ limits. Applying Eq. (40) to data with finite positive $\ell$ is a subject for future work.

III Model checking

Though the application of our model to any given dataset is simple enough in practice, one must validate that the modeling assumptions made in Sec. II are appropriate for the system under consideration. If this model is inappropriate, then the estimates of $\bm{\theta}$ and $Q$, along with the truncation error distributions, could be biased or even meaningless. Conversely, if one assumes our model describes how an EFT should behave, then it would be of great interest to know when a given EFT is failing. Thus a well chosen set of model checking diagnostics are critical to test our assumptions against the data at hand.

To test our assumptions against reality, we must first enumerate them:

The coefficients $c_{n}$ are i.i.d. realizations of a GP. 2. 2.

The $c_{n}$ GP’s mean and covariance functions—and their parameters $\bm{\theta}$—have been correctly identified. 3. 3.

From the knowledge of a few coefficients, we can create a statistically meaningful distribution for the truncation error.

Model-checking diagnostics can assess whether these assumptions are violated at a statistically significant level. The question of how to diagnose issues with GP fitting and prediction has been addressed thoroughly in Ref. Bastos and O’Hagan (2009). However, our situation is novel compared to that of the everyday GP practitioner: (1) they typically have points along one curve, we have multiple i.i.d. curves; (2) their main focus is often the accuracy of the GP regressor’s predictions, our focus is mainly on the parameters $\bm{\theta}$ themselves; and (3) our model for the $c_{n}$ leads to a truncation error distribution which itself can be tested. Nevertheless, the same techniques discussed in Ref. Bastos and O’Hagan (2009) are directly applicable here.

Throughout the discussion below, we assume that the parameters for a GP have either been fixed a priori, or calibrated/marginalized via training data $\mathbf{f}_{\textup{train}}$. The resulting GP (or TP) is then to be evaluated against test (or validation) data $\mathbf{f}_{\textup{val}}$. One should take care that the training set $\mathbf{f}_{\textup{train}}$ does not overlap with any points in $\mathbf{f}_{\textup{val}}$. We have greater than average freedom with how we partition $\mathbf{f}_{\textup{train}}$ and $\mathbf{f}_{\textup{val}}$ because our data contain multiple i.i.d. curves $c_{n}$ that could be tested against one another and because our truncation error model can be tested against reality. For now, assume that $\mathbf{f}_{\textup{train}}$ contains data on one curve and was used to construct an emulator via Eqs. (11)–(13) to predict $\mathbf{f}_{\textup{val}}$ from that same curve. We will return to the details on how one might choose different emulators, or how to partition the training/test sets differently, in Sec. III.4. Below, we use the generic notation of $\mathbf{m}$ and $K$ for the mean and covariance of the emulator (with $\nu$ degrees of freedom in the TP case).

III.1 Mahalanobis Distance

The common way to measure the loss, or incorrectness, of a prediction is through the sum of squared residuals. This metric as a summary assumes that the errors at different input points ${x}$ are completely uncorrelated; the individual weighted residuals at each point are simply summed up. The Mahalanobis distance is the multivariate analog of this idea that takes into account these correlations. That is,

[TABLE]

Large $\textup{D}_{\textup{MD}}^{2}$ implies that the validation data does not match the emulator. But what defines “large” and “small” in this context? A necessary component of any diagnostic is a reference distribution. A reference distribution sets the scale of variation and defines what normal and abnormal diagnostic values look like.

Suppose that we have validation data at $M$ locations. Then, if the emulator is Gaussian, the reference distribution for $\operatorname{pr}(\textup{D}_{\textup{MD}}^{2})$ is a $\chi^{2}$ distribution with $M$ degrees of freedom. If the emulator is a Student-$t$ distribution with $\nu$ degrees of freedom, the reference distribution is an $F$ distribution $F_{M,\nu}$ that is scaled by $s=(\nu-2)M/\nu$ Bastos and O’Hagan (2009). (That is, $\textup{D}_{\textup{MD}}^{2}/s$ follows a standard $F_{M,\nu}$ distribution.) The reference distributions in either the GP or TP cases do not depend on the specific form of the means and covariance functions: only the number of validation data and, if applicable, the degrees of freedom. Both the $\chi^{2}$ and the $F$-distributions are available in standard statistical libraries, such as SciPy, and permit simple analytic evaluations of credible intervals.

But, in general, the distributions of such diagnostics may not have a simple form. A surefire way to generate a distribution in any case is to simply sample from the emulator and compute the diagnostic for each $\mathbf{f}$ in the sample. Then the validation data $\mathbf{f}_{\textup{val}}$ could then be rejected at a 68% or 95% level, for example, based on the empirical reference distribution.

III.2 Pivoted Cholesky Decomposition

The MD is great as a one-number summary, but more information can be gained by considering decompositions of this quantity. Let $G$ be defined by $K=G{G}^{\intercal}$. Then

[TABLE]

is a variance decomposition of $\textup{D}_{\textup{MD}}^{2}=\mathbf{D}_{G}^{\intercal}\mathbf{D}_{G}$. There is freedom in the exact definition of this diagnostic because $G$ is not unique. The pivoted Cholesky decomposition of $K$ leads to a diagnostic ($\mathbf{D}_{\textup{PC}}$) that pinpoints the data contributing to a failing $\textup{D}_{\textup{MD}}^{2}$ and shows useful patterns when plotted graphically vs. index Bastos and O’Hagan (2009). A misestimated variance leads to unusually sized $\mathbf{D}_{\textup{PC}}$ across all indices, whereas an incorrect estimation of the length scale leads to failing $\mathbf{D}_{\textup{PC}}$ at large index.

The reference for each component of $\mathbf{D}_{G}$ is distributed as a standard Gaussian (for a GP emulator) since the points have been scaled and de-correlated by $G$. For the TP emulator, the reference is a Student-$t$ distribution with $\nu$ degrees of freedom. The desired credible interval (here, $2\sigma$) can be computed accordingly.

III.3 Credible Interval Diagnostic

The credible interval diagnostic tests the accuracy of the emulator. An emulator is accurate if a $100\alpha\%$ credible interval (CI) contains approximately $100\alpha\%$ of the validation data. This diagnostic could be particularly useful when comparing experimental data $\mathbf{y}_{\textup{exp}}$ to the full prediction with truncation uncertainty, since it captures whether the error behaves as advertised. Following Ref. Bastos and O’Hagan (2009), start by building $\operatorname{CI}_{i}(\alpha)$, the $100\alpha\%$ marginal CI at the point ${x}_{i}$, for each of the $M$ corresponding $f_{i}\in\mathbf{f}_{\textup{val}}$. Then compute the accuracy of the credible intervals:

[TABLE]

where $\mathbf{1}$ is an indicator function, i.e., it equates to 1 if its argument is true and equates to 0 otherwise. Note that although the $\operatorname{CI}$ may be informed by training data at many points, the $\textup{D}_{\operatorname{CI}}$ is computed point by point.

This diagnostic is the analog to the “consistency plots” shown in Ref. Melendez et al. (2017), which used the pointwise truncation model. It was shown that in this pointwise case, where the data are assumed uncorrelated with one another, the reference distribution is a simple binomial. Since the data here are curves, the reference is more complicated. Here we compute the reference distribution by simulation, that is, we sample the emulator and compute $\textup{D}_{\operatorname{CI}}$ for each sample. The creation of a reference in this manner takes into account correlations by widening the acceptable region of accuracy, though the diagnostic is still inherently a point-by-point summary. Interestingly, we have found that the simulated reference is approximately a binomial with an effective number of data $M_{\text{eff}}=\Delta{x}/\ell$, where $\Delta{x}$ is the difference between the minimum and maximum values of ${x}$.

III.4 Honorable Mentions

We have tested many more diagnostics than are shown graphically in this work, such as

• •

Variogram: A common tool used to investigate correlation structure by looking at the variability of pairs of points across the domain Cressie (1985); Bowman and Crujeiras (2013). We have found that it is too variable for our application.

• •

KL Divergence: Another one-number summary for emulators. We did not find that it gave any new information beyond $\textup{D}_{\textup{MD}}^{2}$.

• •

Other variance decompositions: In addition to the $\mathbf{D}_{\textup{PC}}$, we considered the standard Cholesky decomposition and the eigendecomposition. The Cholesky decomposition has the benefit that each index corresponds to a single validation datum, but it does not order the indices in a visually intuitive manner for diagnosing problems. The eigendecomposition, on the other hand, displays similar qualities as the pivoted Cholesky decomposition when plotted vs. index, but its indices do not correspond to individual validation points. Thus, we choose to display only the pivoted Cholesky here, which enjoys the best of both worlds.

The above diagnostics are still available in the gsum Melendez package despite not appearing here.

III.5 How to Choose an Emulator

We have described the most common approach for diagnosing GP emulators: $\mathbf{f}_{\textup{train}}$ and $\mathbf{f}_{\textup{val}}$ are function values from one curve evaluated at mutually exclusive sets of input points ${x}$. Now we discuss alternatives to this approach based on the questions our diagnostics are designed to address. This discussion is split into two parts: (1) evaluating coefficients $c_{n}$ and (2) evaluating predictions $y_{n}$.

By testing the coefficients we can address the question of whether the observed convergence pattern matches our assumptions. Because they are assumed to be i.i.d., the parameters of their process can be learned from training data along each $c_{n}$ simultaneously or along a subset of the curves. Then an emulator for some $c_{n}$ could be used to predict test data between the interpolating points, regardless of whether those points along $c_{n}$ were used to estimate $\bm{\theta}$.

A problem occurs when using interpolating processes if the training data are too close Ababou et al. (1994); Andrianakis and Challenor (2012); Mohammadi et al. (2016); Pepelyshev (2010); Ranjan et al. (2011). The bubble-like interpolation uncertainty becomes very small, and can approach the size of the nugget used for numerical stability. In this case the diagnostics become highly sensitive to the nugget regulator, a clear sign of diagnostics gone awry. One way to avoid this issue is to spread the training points out relative to the length scale of the system, but this solution comes with a price. The estimates of $\bm{\theta}$ and the accuracy of the interpolants will suffer.

But the issue of numerical instability can be mitigated by avoiding interpolating processes altogether. Our goal is to test whether the $c_{n}$ are drawn from one underlying process, so why not test the coefficients against this process rather than interpolating? The marginal variance of the underlying process should always be large with respect to the nugget size, which sidesteps the problem of comparable sizes. Because the emulator is not forced to interpolate the training data (only $\bm{\theta}$ is tuned), it may be less crucial to exclude training data from the validation set and would allow for more data to be used in the analysis. Diagnosing against the underlying process is thus the method of choice in this work.

Now we outline the diagnosis procedure for $y_{n}$. The predictions $y_{n}$ allow us to test whether our truncation (and possibly interpolation) model behaves as advertised. The data used to estimate the parameters $\bm{\theta}$ could include all EFT predictions up to the highest order available, with the caveat that data that are too closely spaced could cause numerical instability in their estimates. Then an emulator for the full prediction can be constructed as discussed in Sec. II.4, and the validation data could then consist of experimental data $\mathbf{y}_{\textup{exp}}$. If testing the truncation error, there is no issue computing diagnostics with the same ${x}$ as was used in the training set.

Modulo the possible numerical issues discussed above, poorly behaving diagnostics are a sign that something is wrong with the EFT or the convergence model assumptions. To diagnose such problems requires a working knowledge of what success and failure look like. Hence, it will be more illuminating to see the diagnostics in action rather than to further discuss theory. After the following illustrations it should be clear that statistics is a powerful tool to diagnose physics and EFT failure modes.

IV Toy Application

Figure 1 is a summary of how one takes the theoretical predictions of an EFT and extracts coefficients, which then inform the truncation error model as described in Sec. II. This section takes the predictions from Fig. LABEL:sub@fig:pred_to_coeffs_preds as given, and provides an example workflow for extracting physical insights from the data. By working with toy predictions, we are able to assess the accuracy of both the truncation error bands and the extracted GP parameters $\bm{\theta}$. We test the distribution of the coefficients in Fig. LABEL:sub@fig:pred_to_coeffs_coeffs against our GP model, and then compute the truncation error distributions for each $y_{n}$ in Fig. LABEL:sub@fig:pred_to_coeffs_preds, followed by an assessment against a to-all-orders EFT prediction. Throughout, we show visually how the diagnostics from Sec. III behave and how to know when they are failing.

Let us assume that mapping from predictions $y_{n}$ to coefficients $c_{n}$ is well understood, that is, the values of $y_{\mathrm{ref}}$ and $Q$ are known. (We will return later to relax this assumption.) Then we are left with two projects: estimating the process from whence the coefficients were drawn, and subsequently using that process to compute the truncation error distribution, all the while validating assumptions along the way. The first step is to outline the priors on the GP parameters $\bm{\theta}$. Here we take $\eta_{0}=0$ and $V_{0}=0$ (i.e. the mean of the coefficients is known to be zero), which corresponds to the belief that the coefficients are just as likely to be positive as negative and matches the distribution from which they are drawn. The prior on the variance $\bar{c}^{2}$ is given $\nu_{0}=1$ and $\tau_{0}^{2}=1$, which is fairly uninformative with a mode at $1/3$ and a very heavy tail. Finally, an improper uniform prior on the interval $(0,\infty)$ is taken for the length scale $\ell$.

Once priors have been assigned to $\bm{\theta}$, the next step is to update our beliefs using the coefficients as data. We choose 5 points along each of the 4 $c_{n}$ with $\Delta{x}\approx 0.25$, shown as dots in Fig. LABEL:sub@fig:toy_coeffs_and_diagnostics_coeffs, and use them to compute MAP values for $\bar{c}^{2}$ and $\ell$ [Eqs. (30) and (101)]. This optimization procedure results in estimates that are remarkably close to the true parameter values, see Table 1. Plugging these estimates into Eq. (18), yields an estimate of the underlying process that generated the $c_{n}$, whose mean and marginal variance are plotted on Fig. LABEL:sub@fig:toy_coeffs_and_diagnostics_coeffs as black solid and dashed lines. The interpolating processes can then be computed for each $c_{n}$ [Eqs. (12)–(13)], whose means and marginal variances are given by colored dashed lines and filled areas in Fig. LABEL:sub@fig:toy_coeffs_and_diagnostics_coeffs.

We know that the estimated process is quite close to the generating process, but we have no such guarantee when considering real EFT predictions. This is where the diagnostics proposed in Sec. III come into play. We choose to compare the $c_{n}$ validation data, whose locations are indicated by minor ticks in Fig. LABEL:sub@fig:toy_coeffs_and_diagnostics_coeffs, to the underlying GP, as opposed to each interpolating process. Thus, our diagnostics will indicate whether these data points could have feasibly been drawn from the GP that we have learned from the training points. The results of the diagnostics are shown in Figs. LABEL:sub@fig:toy_coeffs_and_diagnostics_md and LABEL:sub@fig:toy_coeffs_and_diagnostics_pc. As expected, the diagnostics behave quite well because the emulator matches the generating process. The $\textup{D}_{\textup{MD}}^{2}$ does return a small value for $c_{2}$, reflecting our intuition that $c_{2}$ looks abnormally small compared to the others. The pivoted Cholesky decomposition vector, plotted vs. its index, shows an essentially random variation of points, meaning that the variation of the $c_{n}$ about the mean are distributed approximately as expected. Again, $c_{2}$ is rather small, given that the vertical variation should follow a standard normal, but the others show a mixture of large and small values of $\mathbf{D}_{\textup{PC}}$. Such diagnostics show that the coefficients from the computed EFT orders conform to the assumptions of our convergence model.

The computed EFT orders are only half of the story; we desire truncation bands that take into account all not yet computed orders. By inserting point estimates of $\bm{\theta}$ and $Q$ into Eq. (21) [or by using its marginalized analog, Eq. (31)], posteriors for the full EFT summation are just as simple as the posterior for the $c_{n}$. The truncation error bands for each $y_{n}$, using the point estimates in Table 1, are shown in Fig. LABEL:sub@fig:toy_truncation_estimates_bands. As an example, a higher-order prediction $y_{20}$ is taken as the true value of the prediction. Given the truncation error processes and the validation data of $y_{20}$, we can again compute diagnostics to determine whether the truncation errors behave as expected. We show the credible interval diagnostic in Fig. LABEL:sub@fig:toy_truncation_estimates_diagnostics, which assesses the accuracy of the truncation error credible intervals. On average, the error estimates contain the true value of $y$ in the proportion dictated by the credible interval. This is the sign of a healthy EFT convergence pattern.

Thus far we have assumed for simplicity that the pathway from order-by-order predictions to naturally sized coefficients, via Eq. (4), was known in advance. But first there is some work to do to even make this conversion: assumptions must be made for $y_{\mathrm{ref}}$ and $Q$. As we have stated previously, $y_{\mathrm{ref}}$ can be derived based on dimensional analysis for the observable at hand, and $Q$ can be based on a separation-of-scales argument specific to the EFT in use. But this can be easier said than done. We have found that one of the most informative pathways to physics insight is graphically exploring a dataset. Plotting coefficients generated from different assumptions on $y_{\mathrm{ref}}$ and $Q$ can show very different patterns, some of which better conform to the assumption that the $c_{n}$ are identically distributed. Empirical exploration of coefficients’ behavior can often lead to useful insights regarding the parameters that enter the EFT convergence model. For example, the assignment $y_{\mathrm{ref}}=1$ for spin observables in Ref. Melendez et al. (2017) was first hit upon through such exploration. This process should be iterative, with scale arguments from the EFT reinforcing or critiquing such empirical observations. The diagnostics we have presented here lend statistical rigor to this process.

An important tool for exploring the EFT convergence is the model evidence Wesolowski et al. (2016) $\operatorname{pr}(\vec{\mathbf{y}}_{k}\,|\,\ell,Q,y_{\mathrm{ref}})$, given by Eq. (95) (where dependence on $y_{\mathrm{ref}}$ is shown explicitly here). Up to prior factors and an unimportant normalization constant, the model evidence is equivalent to the posterior for $\ell$, $Q$, and $y_{\mathrm{ref}}$. Suppose that $\ell$ and $Q$ are both scalar quantities with uniform priors and $y_{\mathrm{ref}}$ is known. Then the posterior $\operatorname{pr}(\ell,Q\,|\,\vec{\mathbf{y}}_{k})$ can be computed analytically, as is shown in Fig. 9 and summarized in Table 1. Figure 9 is a pathway to learning the convergence pattern, and other EFT quantities, through data. In fact, even if $Q$ is a function of ${x}$ and other parameters, such as the EFT breakdown scale $\Lambda_{b}$, then the posterior for these parameters is still analytic and can be optimized to find a MAP value [see Eq. (102)].

But choosing $Q$ may not be as simple as finding the best $\Lambda_{b}$, rather, one may need to decide between a discrete set of functional forms for $Q({x})$. In this case, differing $Q({x})$ can be thought of as models to choose between. Moreover, one could question the choice, or functional form, of $y_{\mathrm{ref}}$. The model evidence is useful as a diagnostic tool in this case as well. ($\ell$ could be integrated out numerically, if desired.) In this case, one can compute Bayes’ factors, i.e., ratios of the model evidence for different choices $(Q,y_{\mathrm{ref}})$, to determine the choice statistically favored by the EFT convergence pattern Trotta (2008); Wesolowski et al. (2016).

As an example of what can occur when $Q$ is misestimated, we now choose $Q=0.3$ (opposed to the true $Q=0.5$) and repeat the analysis of Fig. 7. The results of this analysis are shown in Fig. 10. Since $Q$ is underestimated, then the $c_{n}$ grow with $n$, as dictated by Eq. (4). Thus $c_{3}$ is large compared to the others, and causes the $\bar{c}$ estimate to increase. In general, the marginal variance is biased upwards by large outliers. This is reflected by the Mahalanobis distance only capturing $c_{3}$, but the other coefficients are abnormally small. The assumption that the curves are identically distributed has broken down. The pivoted Cholesky decomposition in Fig. LABEL:sub@fig:toy_coeffs_and_diagnostics_pc_fail provides further evidence of this breakdown; the lower-order components are too small compared to the reference distribution. Moreover, there is a trend of smaller errors near the right side of the chart: evidence that the length scale has been misestimated. Indeed, the length scale found by optimizing the likelihood is $\ell=0.17$, small enough compared to the true value of $\ell=0.2$ to show up in this diagnostic.

V Application to NN scattering with chiral EFT potentials

Now that we have introduced both the truncation error model and its diagnostics, we dedicate this section to a real-world example in low-energy nuclear physics. In this regime, chiral EFT is a popular tool used to develop a systematically improvable description of nuclear observables. The most popular version of chiral EFT is the original proposed in Weinberg’s seminal works, where a nuclear potential is written down and resummed using the Lippmann-Schwinger (LS) equation Weinberg (1990, 1991). This resummation obscures the power counting that was manifest in the potential. Nevertheless, for certain regulators, the convergence of nucleon-nucleon scattering observables appears consistent with this power counting Furnstahl et al. (2015); Melendez et al. (2017).

But the promising results in Refs. Furnstahl et al. (2015); Melendez et al. (2017) were obtained using the pointwise convergence model. Though it permitted a statistical extraction of the EFT breakdown scale $\Lambda_{b}$, strong conclusions about $\Lambda_{b}$ could not be made due to the lack of correlations inherent in the pointwise approach. Here we readdress a subset of these analyses using our GP approach. A more thorough treatment of this and other systems is a topic for future work.

We consider three neutron-proton ($np$) scattering observables: the differential cross section $\sigma(\theta)$ as a function of center-of-mass angle $\theta$, the total cross section $\sigma_{\text{tot}}$ as a function of lab energy $E_{\rm lab}$, and the spin observable $A$ as a function of $\theta$. Each of the observables are computed using the semilocal potential of Epelbaum, Krebs, and Meißner (EKM) Epelbaum et al. (2015a, b). We find that, for $NN$ data in chiral EFT, the coefficients look very smooth and center around a mean of zero, see Ref. Melendez et al. (2017). Thus, we model the correlation structure using the squared exponential kernel, Eq. (9), and choose a mean function of zero.

Once the coefficients’ mean and covariance functions have been specified, choices must be made for $y_{\mathrm{ref}}$ and $Q$. In Ref. Melendez et al. (2017), $y_{\mathrm{ref}}=y_{0}$ was used for the total and differential cross section, while it was argued that $y_{\mathrm{ref}}=1$ is the appropriate scale for spin observables. Here we argue for a slight modification: we choose $y_{\mathrm{ref}}=y_{5}$ for the differential cross section since the leading-order prediction $y_{0}$ gets dangerously close to 0 in some regions, which causes unnatural peaks in the coefficient curves.

The expansion parameter $Q$ is a ratio of low- to high-energy scales. The high-energy scale is the breakdown scale of the EFT $\Lambda_{b}$, whose value is obscured by the LS equation. For $np$ scattering in chiral EFT, the low-energy scales include quantities such as the beam kinetic energy in the lab frame $E_{\rm lab}$ and the pion mass $m_{\pi}$. In the center of momentum frame, one can rewrite the lab energy as the relative momentum given the proton mass $M_{p}$ and neutron mass $M_{n}$

[TABLE]

With $p_{\rm rel}$ a function of $E_{\rm lab}$, the expansion parameter is defined as

[TABLE]

where $f(p_{\rm rel},m_{\pi})$ is some mapping of $p_{\rm rel},m_{\pi}$.

Again, due to the resummation performed by the LS equation, the exact mapping $f(p_{\rm rel},m_{\pi})$ is not known. One would expect that $\lim_{p_{\rm rel}\gg m_{\pi}}f(p_{\rm rel},m_{\pi})=p_{\rm rel}$ and $\lim_{p_{\rm rel}\ll m_{\pi}}f(p_{\rm rel},m_{\pi})=m_{\pi}$, but many functional forms could satisfy these constraints. In prior work, we have assumed the mapping $f(p_{\rm rel},m_{\pi})\to\max(p_{\rm rel},m_{\pi})$, or a smoothed max function

[TABLE]

with $n=8$ Melendez et al. (2017). Here we choose the smooth max function, Eq. (46).

Let us begin by testing the convergence pattern of the EFT for the total and differential cross sections. To do so, we follow the path set out in Sec. IV, that is, we extract the observable coefficients using the physically motivated assumptions of $y_{\mathrm{ref}}$ and $Q$ described above, and evaluate their features using the diagnostics of Sec. III. The analysis of the differential cross section evaluated at $E_{\rm lab}=150$ MeV is given in Fig. 11. The coefficients look beautifully i.i.d. and stationary until $\theta$ approaches backward angles, at which point $c_{3}$ and $c_{4}$ grow in size. This observation is reflected in the diagnostics computed at the validation points. It is clear that both $c_{3}$ and $c_{4}$ are outliers by the Mahalanobis distance, whereas $c_{2}$ and $c_{5}$ are appropriately sized. The pivoted Cholesky provides some insight into the problem: unusually sized values at high index point to problems in the correlation structure or non-stationarity. And it turns out that these outliers at large index are exactly the validation points at large $\theta$, as can be verified by removing all points greater than $\theta=125^{\circ}$ from the diagnostic. If the backward angles are excluded, both diagnostics show excellent agreement with our model assumptions.

Now we turn to the analysis of the total cross section, which was extensively studied in Refs. Furnstahl et al. (2015); Melendez et al. (2017) for differently regulated potentials. The coefficients and their diagnostics are shown in Fig. 12. The choice of test and validation points merit mention here. The highly non-stationary behavior at low energy is not representative of the high-energy data and would bias the analysis if a stationary process were used. This is a standard artifact that appears in our tested $NN$ observables when plotted against lab energy Melendez et al. (2017). Hence these points were omitted from this preliminary analysis. Additionally, the training and test points must be spaced out—due to the very large length scale in this system—otherwise numerical instabilities appear. Given these caveats, the analysis of the total cross section nevertheless show interesting patterns. The diagnostics show that $c_{5}$ is a clear outlier, even when restricted to high energy data. The pivoted Cholesky points again to the onset of non-stationarity or a misestimation of the length scale.

How does chiral EFT fare when subjected to these diagnostics? Both the total and differential cross sections show behavior that is consistent with our convergence model in some regimes, yet inconsistent in others. In this case, there are plausible physical explanations for the diagnostics of unusual size. Likely explanations for the failures of the differential and total cross sections are a misspecification of $Q$ at backward angles and low energies, respectively. This problem with the choice of $Q$ will then result in incorrect truncation error estimates. It may be that the correct expansion parameter takes into account the momentum transfer for angle-dependent observables like the differential cross section. Meanwhile, the non-stationarity appearing at low energies occurs in the regime where $p_{\rm rel}\approx m_{\pi}$, pointing to the possibility that the crossover region of $Q$ was not parameterized correctly. Barring the wrong choice of $Q$, these diagnostic failure modes could simply point to the possibility that the EFT convergence pattern fails in these regimes according to the definition of convergence proposed in this work.

We have shown how observable coefficients can be analyzed in real EFT predictions; truncation errors follow as in Sec. IV. We illustrate this for the spin observable $A$, which is constrained to be zero at $\theta=0^{\circ}$ if the (small) magnetic-moment interaction is neglected, as it is here. Therefore contributions from all orders of the EFT are zero as well. We can accommodate the constraint using the formalism described in Sec. II.4 and shown in Fig. LABEL:sub@fig:gp_applications_constraint. Two sets of EFT truncation errors are proposed in Fig. 13, without and with the constraint. The bands in Fig. LABEL:sub@fig:spin_obs_coeffs_and_truncation_unc are unrealistically large at $\theta=0^{\circ}$ given the information we know about this system. On the other hand, the bands in Fig. LABEL:sub@fig:spin_obs_coeffs_and_truncation_constrained vanish at $\theta=0^{\circ}$ as expected. This example uses $E_{\rm lab}=96$ MeV, but the point is general. The truncation errors are smoothly connected between the constrained and unconstrained regions in a manner harmonious with the way that lower-order coefficients approach zero. That is, the extent of the domain affected by this constraint is dependent on the length scale $\ell$ extracted from the known coefficients.

VI Summary and Outlook

Order-by-order predictions of a well-behaved EFT should converge regularly toward the all-orders value. We have formalized this idea into a falsifiable EFT convergence model. Rigorous estimates of the truncation error—induced by stopping the EFT expansion at a finite order—are now made possible, as are novel techniques for estimating EFT-related quantities, such as the expansion parameter $Q$, breakdown scale $\Lambda_{b}$, and low-energy constants (LECs). We believe that formalizing and testing convergence is a crucial step when using an EFT, otherwise failure may go unrecognized.

This work answers several open questions faced by EFT practitioners that were left unresolved by our earlier work Furnstahl et al. (2015); Melendez et al. (2017); Wesolowski et al. (2019).

(1) EFTs make predictions $y_{n}({x})$ that may be correlated in ${x}$, which in turn induces correlations in the truncation errors; how should that be taken into account when making predictions and fitting? With the EFT power counting as our guide, we propose a change of variables that allows us to isolate the correlated quantities: observable coefficients $c_{n}({x})$ (see Fig. 1). The truncation error distribution follows by modeling the $c_{n}$ as draws from a Gaussian process (GP) and summing all higher order contributions. Originally, we had treated predictions at each kinematic value in a pointwise manner Furnstahl et al. (2015); Melendez et al. (2017) and later modeled correlations in two extreme limits, which correspond to zero and infinite GP length scale $\ell$ Wesolowski et al. (2019). Now $\ell$ can be estimated using data across multiple curves. This GP model could be adapted to correlations between discrete predictions as well.

(2) Is this Bayesian formalism computationally demanding? Bayesian updating of our convergence model parameters $\bm{\theta}$ is made simple and analytic by the use of conjugate priors. Our self-contained treatment shows how prior information on naturalness and correlations is updated using all available data, and derives a novel analytic posterior for the EFT expansion parameter $Q$ (see Fig. 9 and Appendix A).

(3) What about computationally expensive systems? If it is too expensive to compute an EFT prediction at all points ${x}$ or orders $n$, then $\bm{\theta}$ and $Q$, which control the truncation error distribution, can still be estimated from the available predictions without problem. (Even for inexpensive systems, if some low-order corrections are known to be unrepresentative of those at higher orders, one would not want to use these in estimating the truncation error; for example, $c_{1}$ vanishes by symmetry in chiral EFT.) In these cases, the posteriors of $\bm{\theta}$ and $Q$ are simply more influenced by the priors. After estimating $\bm{\theta}$ and $Q$, our convergence model then seamlessly combines truncation estimates with fast interpolation formulae to make predictions at all desired ${x}$ (see Fig. 6 and Sec. II.4, which also shows how to incorporate symmetry constraints).

(4) How do these results affect EFT fitting? We derived a novel likelihood function [Eq. (40)] that enables the EFT fit to account for all higher-order terms and their correlations in ${x}$. This extends to fitting with interpolation as well, allowing more experimental data to be used than might otherwise be computationally feasible. Because it is Gaussian, this likelihood can easily be incorporated into existing fitting procedures.

(5) How is this falsifiable? We provide a menu of model checking diagnostics that enable anyone to assess whether an EFT is working as advertised. We emphasize that these diagnostics require reference distributions (whose functional forms are given in Sec. III.5) to assign statistical significance to their output. Claims about EFT convergence or truncation errors can then be falsified at any desired level of significance. We demonstrate the usefulness of these diagnostics on a toy example in Sec. IV, which exemplifies the general workflow of testing an EFT. Section V does the same for selected $NN$ scattering observables in chiral EFT. The success and failure modes within chiral EFT point to new research directions to pursue. This ability to probe for failure using statistical diagnostics and visualization means that our framework is not just a way to establish theoretical error bars, but provides powerful tools for identifying limitations in our statistical model and in the physical understanding encoded in an EFT.

The stage is now set for a wide range of applications within chiral EFT; some of the paths in progress are

• •

A reevaluation of the EFTs critiqued with the pointwise convergence model in Ref. Melendez et al. (2017) in light of the GP model. Additionally, new variants of chiral EFT, which may show promising convergence properties in the $NN$ sector and fewer regulator artifacts Reinert et al. (2018); Entem et al. (2017), should be included in the analysis.

• •

A systematic study of EFT convergence beyond $NN$ observables to assess its validity. This includes few-body, nuclear matter, and Compton scattering observables.

• •

An exploration of the effects of our proposed likelihood [Eq. (40)] on the LECs’ posteriors, and in particular the effect on three-body LECs.

The models and tools we have developed apply generally to other EFTs and perturbative theories in general, so we anticipate many other applications. The pointwise and curvewise models and associated model checking diagnostics are implemented in the gsum Melendez Python package hosted at github.com/buqeye/gsum and can be freely downloaded and used. This includes a notebook that creates all figures shown here, among other examples.

Acknowledgements.

We thank Derek Bingham for useful discussions. JM thanks Akira Horiguchi and John Honaker for contributions to an early draft of this project. We thank Natalie Klco for contributions to the early stages of, and discussions on, this work. The work of RJF and JAM was supported in part by the National Science Foundation under Grant No. PHY–1614460 and the NUCLEI SciDAC Collaboration under US Department of Energy MSU subcontract RC107839-OSU. The work of JAM was also supported in part by the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education for the DOE under contract number DE-SC0014664. The work of DRP was supported by the US Department of Energy under contract DE-FG02-93ER-40756 and by the ExtreMe Matter Institute EMMI at the GSI Helmholtzzentrum für Schwerionenphysik, Darmstadt, Germany. The work of SW was supported in part by Salisbury University.

Appendix A Derivations

The mathematical details of the model described in Sec. II are shown here. Due to the choice of conjugate priors on $\mu,\bar{c}^{2}$, almost all posteriors of interest can be provided in closed form.

Before getting into the fine details of the GP model discussed in this work, we believe it instructive to begin with a short aside. We have considered a simpler pointwise model in previous works Furnstahl et al. (2015); Melendez et al. (2017), which is shown graphically in Fig. 3. To our knowledge, no one has yet worked out analytic expressions for the case of conjugate priors, whose derivations mirror the correlated case. Hence we will provide the results below for completeness. Again, conjugacy allows for efficient computation of truncation errors, parameter posteriors, and evidences using standard statistical libraries. We promote the use of these conjugate priors over the options laid out in previous works.

See Table 2 for a summary of the vector notation. Footnote 4 gives a brief overview of the statistics notation.

A.1 Pointwise Model

Consider a set of EFT predictions $\{y_{n}\}$, each of which is a number rather than a function. In this case, the quantity of interest $y_{n}$ could itself be a scalar or we could be considering a functional quantity $y_{n}({x})$ at a specific input point ${x}$ in its domain. Here we wish to predict the truncation error independently of the function values at any surrounding input points. This procedure could be repeated at all points in the domain to produce error bands, as was done in Ref. Melendez et al. (2017). This is what we call the pointwise model.

If the expansion parameter $Q$ and reference scale $y_{\mathrm{ref}}$ are assumed to be known, then the data $\vec{y}_{k}$ can be converted to coefficients $\vec{c}_{k}$,

[TABLE]

Here we assume that the coefficients $c_{n}$ are normal given $\bar{c}^{2}$ and that $\bar{c}^{2}$ has a scaled inverse-$\chi^{2}$ prior Gelman et al. (2013),

[TABLE]

Equation (49) has qualitatively similar properties to priors used in past works Furnstahl et al. (2015); Melendez et al. (2017): it is strictly positive, can penalize $\bar{c}^{2}$ for being or large or too small, and reduces to a scale-invariant prior for $\nu_{0}=0$. Unlike past priors on $\bar{c}^{2}$, the scaled inverse-$\chi^{2}$ permits intuitive analytic posteriors for $\delta y_{k}$ and $Q$ for all choices of $\nu_{0}$ and $\tau_{0}^{2}$.

Analogously to Eq. (21), one can derive from Eqs. (5), (20), and (48) that

[TABLE]

The above equations outline our prior beliefs about the coefficients, their widths, and how they combine to form a discrepancy term for $y_{k}$. If $\bar{c}$ was precisely known from a priori arguments or otherwise, then Eq. (50) is sufficient to estimate the truncation error. Since this is not the case in general, we would like to update these priors based on order-by-order predictions.

Our belief about an unseen coefficient $c_{n}$ is updated via

[TABLE]

By marginalizing in information about $\bar{c}^{2}$, it is clear that the effect of data on $c_{n}$ flows through $\bar{c}^{2}$. To compute the posterior predictive distribution, Eq. (51), we must first understand how the observed coefficients $\vec{c}_{k}$ impact our understanding of $\bar{c}^{2}$. This is where the utility of conjugate priors is realized. The posterior of $\bar{c}^{2}$ has the same functional form as the prior—an inverse $\chi^{2}$ distribution—but with $\nu$ and $\tau^{2}$ that are updated by the data. This property simplifies our analysis greatly: computing the posterior $\bar{c}^{2}\,|\,\vec{c}_{k}$ is reduced to determining $\nu$ and $\tau^{2}$. From Bayes’ theorem and the fact that the $\vec{c}_{k}$ are marginally independent given $\bar{c}^{2}$,

[TABLE]

Above we have defined $n_{c}$ to be the number of coefficients in $\vec{c}_{k}$ that are useful for induction (which is not always $k+1$, see Ref. Furnstahl et al. (2015)). Upon comparison with Eq. (27), we can read off

[TABLE]

Therefore, we have shown that

[TABLE]

with $\nu$ and $\tau^{2}$ defined by Eqs. (53) and (54).

Now that the posterior for $\bar{c}^{2}$ has been established as Eq. (55), we can return our attention to Eq. (51). This integral evaluates to

[TABLE]

Equation (56) is the pdf for the Student-$t$ distribution with degrees of freedom $\nu$, mean [math], and scale $\tau$, i.e.,

[TABLE]

In contrast to the notation of the Gaussian distribution, the variance of the Student-$t$ is given by $\nu\tau^{2}/(\nu-2)$ and is only defined when $\nu>2$.

Let us pause and ponder the implications of what we have derived. Equations (48), (53)–(55), and (57) fit together within one cohesive picture and relate to common estimators employed in frequentist statistics.

• •

The variance of $\bar{c}^{2}$ decreases as the degrees of freedom $\nu$ increases. This uncertainty in $\bar{c}^{2}$ maps directly to the heaviness of the tails in the $t$ distribution, which is also controlled by $\nu$.

• •

A strict prior ($\nu_{0}\gg 1$) or a lot of data ($n_{c}\gg 1$) both lead to $\operatorname{pr}(\bar{c}^{2}\,|\,\vec{c}_{k})$ being sharply peaked at $\bar{c}^{2}\approx\tau^{2}$. This informative distribution for $\bar{c}^{2}$ corresponds to a $\operatorname{pr}(c_{n}\,|\,\vec{c}_{k})$ that approximates a Gaussian with variance $\tau^{2}$.

• •

Given a lot of data, then $\bar{c}^{2}\approx\tau^{2}\approx\vec{c}_{k}^{\,2}/n_{c}$, which is a frequentist estimator known as the sample variance. In this case $\operatorname{pr}(c_{n}\,|\,\vec{c}_{k})$ will approximate a Gaussian with $\vec{c}_{k}^{\,2}/n_{c}$ as its variance.

• •

Rather than integrating over $\bar{c}^{2}$, we could use its mean value as a point estimate for the variance in Eq. (48). The mean of Eq. (55) is given by $\nu\tau^{2}/(\nu-2)$, which is exactly the variance of the $t$ distribution in Eq. (57).

By computing or assuming values for $\nu$ and $\tau^{2}$, we are able to relax the assumption of a known $\bar{c}^{2}$ in Eq. (48) by giving $c_{n}$ a distribution with tails and a scale compatible with our understanding of $\bar{c}^{2}$. A Bayesian perspective allows us to incorporate prior beliefs about naturalness and update them with finite sets of data in a manner that is harmonious with the limiting cases assumed in frequentist statistics.

A better understanding of unseen observable coefficients $c_{n}$ is useful as a means to better understand the truncation error $\delta y_{k}$. The steps in deriving Eq. (57) lead analogously to the posterior

[TABLE]

where $\nu$ and $\tau^{2}$ are computed via Eqs. (53) and (54). The distribution for the full prediction $y=y_{k}+\delta y_{k}$ is thus

[TABLE]

Equation (59) is a novel result from this work and is the most straightforward method to compute posteriors of $y$ in the pointwise model. Statistical packages such as SciPy have built-in $t$ distributions which can easily calculate degree of belief intervals for Eq. (59).

Until now, we have assumed that the expansion parameter $Q$ is known, but in practice this may not be the case. In fact, $Q$ may even be a function of ${x}$ in certain circumstances, which we assume here for generality. Fortunately, conjugacy permits the analytic computation of the (unnormalized) posterior

[TABLE]

The right-hand side is (up to a simple integral over $Q$) the model evidence, which is a useful quantity for model comparison, and the likelihood $\operatorname{pr}(\vec{\mathbf{y}}_{k}\,|\,\mathbf{Q})$ is exactly calculable. The pointwise model assumes that each of the $\vec{y}_{k}({x}_{i})$ are independent and each have their own $\bar{c}_{i}$, so that $\operatorname{pr}(\vec{\mathbf{y}}_{k}\,|\,\mathbf{Q})=\prod_{i}\operatorname{pr}[\vec{y}_{k}({x}_{i})\,|\,Q({x}_{i})]$. We consider each ${x}_{i}$ individually and drop ${x}_{i}$ for convenience. Note that since $y_{0}=y_{\mathrm{ref}}c_{0}$ and $y_{n}=y_{n-1}+y_{\mathrm{ref}}c_{n}Q^{n}$, we can make a change of variables

[TABLE]

where $c_{0}=y_{0}/y_{\mathrm{ref}}$ and $c_{n}=\Delta y_{n}/y_{\mathrm{ref}}Q^{n}$ implicitly depend on $Q$. Thus we need only compute the joint distribution $\operatorname{pr}(\vec{c}_{k})$. Although the distribution of one coefficient $c_{n}$ follows a $t$ distribution Eq. (57), the same is not true of the joint distribution of multiple coefficients $\vec{c}_{k}$. Nevertheless, we can find an analytic form for their distribution using the clever manipulation of normalization constants. Let primes denote unnormalized distributions. Then both

[TABLE]

from Bayes’ theorem, and

[TABLE]

from Eq. (55), are true statements. Since $\operatorname{pr}(\vec{c}_{k})$ does not depend on $\bar{c}^{2}$, equating the normalization coefficients yields the desired quantity

[TABLE]

Through Eqs. (61) and (64) we have completely specified $\vec{y}_{k}\,|\,Q$, and $\vec{\mathbf{y}}_{k}\,|\,\mathbf{Q}$ is simply the product of these quantities.

Finally, we can write an unnormalized expression for the expansion parameter posterior by dropping terms independent of $Q$ and multiplying $\vec{y}_{k}({x}_{i})\,|\,Q({x}_{i})$ at each ${x}_{i}$

[TABLE]

Rather than dealing with the multidimensional object $\mathbf{Q}$, one can instead parameterize $Q$ as a low-energy scale $f({x})$ and a high-energy scale $\Lambda_{b}$, known as the breakdown scale of the EFT. These are related to $Q$ via $Q({x})=f({x})/\Lambda_{b}$. The posterior for $\Lambda_{b}$, assuming $f({x})$ is known, has a simple relationship to the posterior for $Q$:

[TABLE]

This is a more general version of the posterior given in Ref. Melendez et al. (2017), where the restriction $\nu_{0}=0$ was required. Here, thanks to conjugacy, $\nu_{0}$ and $\tau_{0}^{2}$ are free to be chosen by the modeler. The extension to further parameters of $Q$ requires only a modification of the prior $\operatorname{pr}(\Lambda_{b})$ in Eq. (66).

A.2 Gaussian Process Model

Here we will expand on the template provided by Appendix A.1 and derive analogous results for the GP truncation error model introduced in Sec. II. The results in this Appendix follow exactly from the analysis of the Bayesian linear model (see Ref. O’Hagan (1994), chaps. 9 and 10) but are reproduced here for completeness. As a reminder, we have assumed that the $c_{n}$ are GPs with the following form:

[TABLE]

and placed a normal-inverse-$\chi^{2}$ prior on $\mu,\bar{c}^{2}$,

[TABLE]

We have generalized the parametrization as promised in Sec. II.3 by introducing a length $p$ vector of functions $b({x})$ that multiplies a length $p$ vector of regression coefficients $\mu$ to allow for a non-constant mean function as in, e.g., Ref. Bastos and O’Hagan (2009). Both $\eta_{0}$ and $V_{0}$, previously scalars, are promoted to a length $p$ vector and a $p\times p$ matrix, respectively. For the examples of EFT truncations shown in this work, we use $b({x})=1$.

As before, we would like to condition on data to gain a better understanding of the observable coefficients and ultimately the discrepancy $\delta y_{k}({x})$. Again, assume that $Q$ and $y_{\mathrm{ref}}$ are known so that our data are $N$ vectors of coefficients $\vec{c}_{k}$ (each of length $n_{c}$) at the chosen input points. The combined ($n_{c}\times N$)-shaped observations are denoted $\vec{\mathbf{c}}_{k}$ (see Table 2). The posterior predictive distribution for a new curve $c_{n}({x})$ is then

[TABLE]

where $\bm{\theta}=\{\mu,\bar{c}^{2},\ell\}$. As before, the data influence $c_{n}({x})$ by updating our beliefs about $\bm{\theta}$. Because of the conjugate prior on $\mu$ and $\bar{c}^{2}$, they can be updated and integrated out of Eq. (68) analytically, but there is no conjugate prior for $\ell$. For now, assume that $\ell$ is known; we will return to finding its posterior later. With a posterior for $\ell$, one could use its MAP value as a point estimate in the following equations, or marginalize over it numerically.

To make progress on Eq. (68) requires deriving $\operatorname{pr}(\mu,\bar{c}^{2}\,|\,\vec{\mathbf{c}}_{k},\ell)$. Again, since we know that this posterior has the same functional form as the prior, our task is to determine $\eta,V,\nu,\tau^{2}$. As usual, Bayes’ theorem gets the ball rolling:

[TABLE]

Next, by inserting Eqs. (67) and (24),

[TABLE]

where $R_{\ell}=r({\mathbf{x}},{\mathbf{x}};\ell)$ and $B=b^{\intercal}({\mathbf{x}})$ are $N\times N$ and $N\times p$ matrices, respectively. If $b({x})=1$, as is assumed for the examples here, then $B=\bm{1}$, a length $N$ vector of ones.

So far we can read off the updated value of $\nu=\nu_{0}+Nn_{c}$, but must inspect the exponent further for the others. Define $\expectationvalue{\vec{\mathbf{c}}_{k}}$ as an $N\times 1$ vector that is the average over the $n_{c}$ orders,

[TABLE]

The sum over the quadratic form can be rearranged as

[TABLE]

so that we can complete the square for $\mu$ in the exponent of Eq. (70). The result is

[TABLE]

where

[TABLE]

After some algebra, one can rewrite Eq. (77) in an intuitive way O’Hagan (1994):

[TABLE]

which updates the scale $\tau^{2}$ by combining prior information with two other terms: one due to the discrepancy between the prior mean and sample mean, and the other due to the sample variance $s^{2}$ for correlated observations,

[TABLE]

When $n_{c}\gg 1$, then $\tau^{2}\approx s^{2}$. If one only wanted point estimates so that the $c_{n}$ or $\delta y_{k}$ posteriors remain Gaussian, then $\mu\approx\eta$ and $\bar{c}^{2}\approx\tau^{2}$ from Eqs. (74) and (78) provide approximate Bayesian estimators.

Returning to Eq. (68), we must integrate out $\mu$ and $\bar{c}^{2}$ (again, assuming $\ell$ is fixed). Though the integral over $\mu$ is doable, there is neat trick to see the answer right away: rewrite the GP as a mean $\mu$ (whose distribution is Gaussian) plus a zero-mean GP, that is,

[TABLE]

Now, substitute Eq. (81) into Eq. (80) and use Gaussian sum rules to arrive at

[TABLE]

which has marginalized out $\mu$ for us. The integral over $\bar{c}^{2}$, follows from the multivariate analog of Eq. (57), that is, integrating Eq. (82) over $\bar{c}^{2}$ yields a Student-$t$ process $\operatorname{\mathcal{TP}}_{\nu}[m({x}),\kappa({x},{x}^{\prime})]$ Shah et al. (2014)

[TABLE]

Just like GPs, Student-$t$ processes (TPs) are collections of random variables, any finite number of which have a joint $t$ distribution. Specifically, the density of an $N$ dimensional multivariate $t$ distribution is

[TABLE]

With this notation, $\kappa({x},{x}^{\prime})$ is not the covariance function of $\operatorname{\mathcal{TP}}_{\nu}[m({x}),\kappa({x},{x}^{\prime})]$, and, similarly, $K$ is not the covariance matrix in Eq. (A.2). Rather, the covariance is given by $\nu\kappa({x},{x}^{\prime})/(\nu-2)$ and is only defined for $\nu>2$. This differs from the notation in, say, Refs. Bastos and O’Hagan (2009); Shah et al. (2014). TPs, like GPs, can be used in regression, etc., but can better handle outliers in the data. Here, TPs only show up in our convergence model due to our choice of priors, but are convenient objects to work with nevertheless.

For completeness, we provide the conditional distribution for Student-$t$ processes, which is very similar to the GP version in Sec. II.2.2. Suppose $\mathbf{y}_{1}$ and $\mathbf{y}_{2}$ are length $n_{1}$ and $n_{2}$ vectors, respectively, and

[TABLE]

Then

[TABLE]

where $\tilde{\mathbf{m}}_{2}$ and $\tilde{K}_{22}$ are defined as in Eqs. (12) and (13), and $d_{1}={(\mathbf{y}_{1}-\mathbf{m}_{1})}^{\intercal}K_{11}^{-1}(\mathbf{y}_{1}-\mathbf{m}_{1})$. One could derive this by starting with $\mathbf{y}\,|\,\sigma^{2}\sim\mathcal{N}(\mathbf{m},\sigma^{2}K)$ and $\sigma^{2}\sim\operatorname{\chi^{--2}}(\nu,1)$, using the Gaussian conditional rules described in Sec. II.2.2, and then marginalizing over $\sigma^{2}$ at the end.

It is important to note that Eq. (86) is not how we compute conditionals for, e.g., Eq. (83). This is because conditioning and marginalizing over the mean do not commute. Consider the case where $\mathbf{f}={\begin{bmatrix}\mathbf{f}_{1}&\mathbf{f}_{2}\end{bmatrix}}^{\intercal}$ is normally distributed just as $c_{n}$. Then

[TABLE]

from which it is clear that one must first follow the conditional rules of Eqs. (12)–(13) and also update the posterior of $\mu$ and $\bar{c}^{2}$ before integrating them out. Thus, following the same logic as above,

[TABLE]

where, remembering that $B=\bm{1}$ (a vector of ones) is the basis matrix for our coefficients,

[TABLE]

Of course, there could be other i.i.d. data beyond just $\mathbf{f}_{1}$ that could be used to compute $\operatorname{pr}(\mu,\bar{c}^{2}\,|\,\mathbf{f}_{1})$, which would serve only to adjust $\eta$, $V$, $\nu$, and $\tau^{2}$.

The above analysis is sufficient to finish our derivation for full uncertainty quantification of $y_{k}$, though the details will vary depending on the application as described in Sec. II.4. In all cases, the derivation involves simply using $m_{\textup{th}}$ and $\Sigma_{\textup{th}}$ from Sec. II.4 as the mean and covariance for the Gaussian $\operatorname{pr}(\mathbf{f}_{2}\,|\,\mathbf{f}_{1},\mu,\bar{c}^{2})$ in Eq. (87), and then integrating out $\mu$ and $\bar{c}^{2}$. Below we provide the mean and scale function for $\operatorname{\mathcal{TP}}_{\nu}[m_{\textup{th}}({x}),\Sigma_{\textup{th}}({x},{x}^{\prime})]$ in each case.

Inexpensive Predictions. Let primes denote means with $\mu\to\eta$, e.g.,

[TABLE]

Then

[TABLE]

Expensive Predictions. Define $\tilde{b}_{i}({x})$ for $i\in\{k,\delta k\}$ as the non-vectorized version of Eq. (89), using the appropriate $b_{i}({x})$ and $R_{i}$. Then

[TABLE]

Constraints. One only needs to make the replacements $m_{\delta k}^{\prime}\to\tilde{m}_{\delta k}^{\prime}$, $b_{\delta k}\to\tilde{b}_{\delta k}$ and $R_{\delta k}\to\tilde{R}_{\delta k}$ in the above equations.

In the final part of this appendix, we derive the formulas for the posterior $\operatorname{pr}(\ell,\mathbf{Q}\,|\,\vec{\mathbf{y}}_{k})\propto\operatorname{pr}(\vec{\mathbf{y}}_{k}\,|\,\ell,\mathbf{Q})\operatorname{pr}(\ell,\mathbf{Q})$. The quantity $\operatorname{pr}(\vec{\mathbf{y}}_{k}\,|\,\ell,\mathbf{Q})$ is (almost) the marginal likelihood (or, evidence), and can be computed exactly. The first step is to make a change of variables from the data $\vec{\mathbf{y}}_{k}$ to coefficients

[TABLE]

To find $\operatorname{pr}(\vec{\mathbf{c}}_{k}\,|\,\ell)$, we will make use of the same normalizing constant trick as in the uncorrelated case. Again, let primes denote unnormalized distributions and $Z_{\text{dist.}}$ be their normalizing constants. By Bayes’ theorem,

[TABLE]

and by conjugacy

[TABLE]

Given that

[TABLE]

with a similar form for $Z_{\chi}$, then

[TABLE]

The joint posterior $\operatorname{pr}(\ell,\mathbf{Q}\,|\,\vec{\mathbf{y}}_{k})$ follows from Eqs. (95), (99) and a choice of prior $\operatorname{pr}(\ell,\mathbf{Q})$. It follows that the unnormalized marginal posteriors for $\mathbf{Q}$ and $\ell$ are

[TABLE]

The posterior for the breakdown scale $\Lambda_{b}$, where $Q({x};\Lambda_{b})=f({x})/\Lambda_{b}$, follows from the same logic as in Appendix A.1 and is given by

[TABLE]

MAP values can then be found numerically and used as point estimates in the analytic equations derived in this Appendix.