Non-parametric Bayesian approach to extrapolation problems in configuration interaction methods

TL;DR

This paper introduces a non-parametric Bayesian extrapolation technique using constrained Gaussian processes for configuration interaction methods, enabling accurate predictions with uncertainty quantification, especially effective with small data sets.

Contribution

The paper presents a novel non-parametric Bayesian approach that incorporates physical constraints into Gaussian process models for extrapolation in configuration interaction calculations.

Findings

Effective extrapolation with small data sets

Incorporates physics-guided constraints

Provides quantified uncertainty in predictions

Abstract

We propose a non-parametric extrapolation method based on constrained Gaussian processes for configuration interaction methods. Our method has many advantages: (i) applicability to small data sets such as results of {\it ab initio} methods, (ii) flexibility to incorporate constraints, which are guided by physics, into the extrapolation model, (iii) providing predictions with quantified extrapolation uncertainty, etc. In the present study, we show an application to the extrapolation needed in full configuration interaction method as an example.

Figures

Equations

$$\displaystyle\Sigma=\left[\begin{array}[]{cc}K_{XX}&K_{XX^{*}}\\ K^{T}_{XX^{*}}&K_{X^{*}X^{*}}\end{array}\right].$$

$$k_{logMat52}(x_i,x_j)$$

$$P(\mathbf{y},{\mathbf{y}}^{\ast }|\mathbf{\theta })$$

$$P(\mathbf{y},{\mathbf{y}}^{\ast }|\mathbf{\theta })=P({\mathbf{y}}^{\ast }|\mathbf{y},\mathbf{\theta })P(\mathbf{y}|\mathbf{\theta }).$$

$$P({\mathbf{y}}^{\ast }|\mathbf{y},\mathbf{\theta })$$

$${\mathbf{\mu }}_{{\mathbf{y}}^{\ast }|\mathbf{y}}(\mathbf{\theta })$$

$${\Sigma }_{{\mathbf{y}}^{\ast }|\mathbf{y}}(\mathbf{\theta })$$

$$P(\mathbf{y}|\mathbf{\theta })$$

$$P({\mathbf{y}}^{\ast }|\mathbf{y})$$

$$P({\mathbf{y}}^{\ast }|\mathbf{y},\alpha ,\beta ,\dots )\propto \int$$

$$\times P(\alpha ,\beta ,\dots |{\mathbf{y}}^{\ast },\mathbf{y})d\mathbf{\theta },$$

$$P({\mathbf{y}}^{\ast }|\mathbf{y},\alpha ,\beta ,\dots )\simeq \sum _{i=1}^{N_p}{w}^{(i)}P({\mathbf{y}}^{\ast (i)}|\mathbf{y},{\mathbf{\theta }}^{(i)}),$$

$$w^{(i)}\equiv \frac{P(\mathbf{y},{\mathbf{\theta }}^{(i)})P(\alpha ,\beta ,\dots |{\mathbf{y}}^{\ast (i)},\mathbf{y})}{{\sum }_{j=1}^{N_p}P(\mathbf{y},{\mathbf{\theta }}^{(j)})P({\mathbf{y}}^{\ast (j)}|\mathbf{y},{\mathbf{\theta }}^{(j)})P(\alpha ,\beta ,\dots |{\mathbf{y}}^{\ast (j)},\mathbf{y})}.$$

$$P(\alpha |{\mathbf{y}}^{\ast },\mathbf{y})=$$

$$\times \Phi (y_{P-1}^{\ast }-y_P^{\ast }).$$

$$\Phi (z)$$

$$r(x_j)\equiv$$

$$r(x_1^{\ast })$$

$$r(x_{P-1}^{\ast })$$

$$P(\beta |{\mathbf{y}}^{\ast },\mathbf{y})=$$

$$\cdots \times \Phi \left(R_E-r(x_{P-1}^{\ast })\right),$$

$$R_E$$

$$r_{\mathrm{mean}}$$

$$r_{\mathrm{std}}$$

$$k_{\mathrm{RBF}}(x_i,x_j)$$

$$k_{\mathrm{Mat}\acute{\mathrm{e}}\mathrm{rn}}(x_i,x_j;\nu )$$

$$\xi$$

$$k_{\mathrm{logRBF}}$$

$$k_{logMat52}$$

$$k_{logMat32}$$

$$\mathrm{BIC}=\ln P(\mathbf{y}|{\mathbf{\theta }}_{\mathrm{ML}})-\frac{1}{2}M\ln n,$$

$$\Delta \equiv {\log }_{10}\left(max(|{\mathbf{\mu }}_{{\mathbf{y}}^{\ast }|\mathbf{y}}^{\mathrm{own}}-{\mathbf{\mu }}_{{\mathbf{y}}^{\ast }|\mathbf{y}}^{lib.}|)\right),$$

$$P({\mathbf{y}}^{\ast },\mathbf{y},\mathbf{c})$$

$$P({\mathbf{y}}^{\ast }|\mathbf{y},\mathbf{c})$$

$$=\int d\mathbf{\theta }\frac{P(\mathbf{c}|{\mathbf{y}}^{\ast },\mathbf{y},\mathbf{\theta })P({\mathbf{y}}^{\ast },\mathbf{y},\mathbf{\theta })}{P(\mathbf{y},\mathbf{c})},$$

$$\propto \int d\mathbf{\theta }P(\mathbf{c}|{\mathbf{y}}^{\ast },\mathbf{y})P({\mathbf{y}}^{\ast }|\mathbf{y},\mathbf{\theta })P(\mathbf{y}|\mathbf{\theta })P(\mathbf{\theta }).$$

Full text

††preprint: …

I Introduction

Recent developments in nuclear potentials based on chiral effective field theory (chiral EFT) [1, 2] and ab initio methods [3, 4, 5, 6, 7] have provided deep insights into nuclear many-body problems starting from the fundamental interaction between nucleons, see very recent works [8, 9, 10, 11] and references therein.

The Full configuration interaction (FCI) method, which is also known as no-core full configuration (NCFC)/no-core shell model (NCSM) [12, 13], is one of the successful ab initio methods. In FCI, wave functions are represented in a truncated sub-space, and the truncation is typically specified by the parameter $N_{\mathrm{max}}$ which defines the maximum number of harmonic oscillator quanta allowed in the many-body states above the lowest configuration for a target nucleus.

Despite enormous efforts for developing efficient codes [14, 15, *MFDn2, *MFDn3, 18, *KSHELL2, 20, 21] and advances in computing power, the currently available $N_{\mathrm{max}}$ for upper $p$-shell nuclei is around $10$ (see e.g. [22]) and this is still far from $N_{\mathrm{max}}=\infty$ corresponding to exact calculations. One usually extrapolates the sequence of results with different $N_{\mathrm{max}}$ to $N_{\mathrm{max}}=\infty$ to estimate the exact value. In previous studies, several extrapolation methods were proposed and the dependence on them was analyzed [23, 22, 24, 25, 26, 27, 21]. The most intuitive example is one based on an exponential function [22].

In addition to FCI, such extrapolation techniques are also required in CI calculations for a valence space, which is the so-called shell model, using additional truncations. Representative examples of the truncations are importance truncation scheme [28, *Roth_IT1b, *Roth_IT1c, 31] and Monte Carlo shell model [32, 33]. In those calculations, the rapid growth of many-body basis for a valence space is alleviated by selecting a small subset of the many-body basis states which is physically more relevant. These truncation schemes have been successfully applied to valence CI and also FCI calculations in previous works such as [34, 35, 36, 37, 38, 39].

In these studies, extrapolation is performed with some specific functions such as exponential or polynomials, and its coefficients are determined so as to minimize the $\chi^{2}$ deviation from the given calculated data points. While any of these offers intuitively reasonable extrapolated results, there is a risk of overfitting, i.e., lack of predictive power for true exact values. This overfitting is because $\chi^{2}$ minimization of a parametric model and point estimation of the parameter leads to too deterministic predictions due to the limited expression power of the model. This is a problematic situation if one intends to discuss quantitative issues like a level ordering of states with small energy differences, the positions of proton and neutron drip lines, and so on.

In the present study, we propose a novel non-parametric extrapolation method for CI calculations using constrained Gaussian processes (GPs), which can give extrapolated results with quantified uncertainty in a systematic manner.

Evaluating extrapolation uncertainties are helpful for breaking down possible origins of discrepancy between FCI results and experimental data, though the major source of uncertainty is, at the moment, from the input potentials: the low-energy constants and the truncated expansion in chiral EFT. Although GPs are also not free from overfitting, GPs allow us to consider a wider class of functions, and then it is expected to alleviate underestimating extrapolation uncertainties due to the specific choice of the extrapolation function.

We also note that extrapolation techniques using an artificial neural network (ANN) are proposed recently [40, 41]. To train networks, one usually requires large data sets. However, it is still tough to achieve an enormous number of ab initio calculations while varying their inputs such as $N_{\mathrm{max}}$ and harmonic oscillator parameter $\hbar\Omega$. In future applications of full CI and also valence CI methods with importance truncation to heavier systems, it is strongly desired to develop an extrapolation technique applicable even to sparse data. The proposed method is applicable to small data sets too.

We demonstrate the validity of our model by taking the ground state energies obtained by FCI calculations as an example. The code is available on the author’s GitHub page [42].

II Formulation of constrained Gaussian Process

Gaussian Process (GP) is a popular statistical method as a non-parametric regression model [43]. It is also becoming popular in physics due to its flexibility (see e.g. recent applications published in APS journals [44, 45, 46, 47, 48, 49, 50]). The GP regression can be interpreted as a method to describe distribution over a function space and to perform inference of the probability for each function. This is just what we need because this enables us to consider an ensemble of many possible functions for the extrapolation, infer probability of each sample function, and then quantify uncertainties in the extrapolated value in a statistical manner.

Interestingly, GPs are mathematically equivalent or related to many other models such as ANN, support vector machines, spline models, and so on. We refer the interested reader to e.g. [51, 43, 52].

Here we introduce some notations. As in statistics literature, $P(a|b)$ denotes the probability distribution of $a$ under the condition $b$, and we use $\mathcal{N}(\bm{\mu},\bm{\Sigma})$ to express the multivariate Gaussian distribution with mean vector $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$. In what follows, we consider two variable sets, data and prediction. The terminology data, which is distinguished from experimental data, is used to express a set of $X=\{x_{i}|i=1,...,D\}$ and $Y=\{y_{i}|i=1,...,D\}$. Here we assumed that we have $D$ input points. The prediction represents positions $X^{*}=\{x^{*}_{i}|i=1,...,P\}$ and values $Y^{*}=\{y^{*}_{i}|i=1,...,P\}$ for $P$ points where the target values are not known. Especially in our applications, $X$ denotes currently computable $N_{\mathrm{max}}$, and $X^{*}$ is a set of $N_{\mathrm{max}}$ at which FCI calculations have not done (e.g. larger $N_{\mathrm{max}}$).

In Gaussian processes, it is assumed that the two target values at the two arbitrary points in the vicinity must be similar, and the so-called kernel functions express the similarities. Then the data values $Y$ and prediction values $Y^{*}$ are assumed to be generated from the multivariate Gaussian distribution $\mathcal{N}(\bm{\mu},\bm{\Sigma})$ whose covariance matrix $\bm{\Sigma}$ is given as

[TABLE]

Here $K_{XX}$, $K_{XX^{*}}$, and $K_{X^{*}X^{*}}$ are respectively $D\times D$, $D\times P$, and $P\times P$ matrices, and these elements are evaluated with a kernel function. It is a common practice to choose this kernel function as the radial basis function (RBF) or the Matérn kernel with $\nu=3/2$ (Mat32) or $\nu=5/2$ (Mat52) (see Appendix A).

In this work, we use the logMat52 kernel function for the reasons described in Appendix A. The logMat52 kernel is defined for, e.g. two data points $x_{i}$ and $x_{j}$, as follows:

[TABLE]

where $r\equiv|x_{i}-x_{j}|$, and the global strength $\tau$ and correlation length $\ell$ are the hyperparameters. Let $\bm{\theta}$ denote the vector of hyperparameters. We will revisit the issue of hyperparameters later.

Once the kernel function and its hyperparameters are fixed, one can define the joint covariance matrix $\bm{\Sigma}$ in Eq. (3) for data/prediction as a function of $\bm{\theta}$. Then, the joint distribution of data $\bm{y}$ and predictions $\bm{y}^{*}$ under the hyperparameters is given as

[TABLE]

It is a common practice for mean vectors to be normalized, i.e., $\bm{\mu}$ and $\bm{\mu}^{*}$ are fixed as zero vectors and the data is scaled to have zero mean and unit variance. The dependence on the choice of mean vectors is also discussed later.

By definition of conditional probabilities, the left-hand side of Eq. (7) can be rewritten as

[TABLE]

Under given $\bm{\theta}$, one can write down $P(\bm{y}^{*}|\bm{y},\bm{\theta})$ and $P(\bm{y}|\bm{\theta})$ in a closed form:

[TABLE]

Regarding the hyperparameters, the so-called maximum a posteriori (MAP), i.e. one to maximize the hyperparameter posterior $P(\bm{\theta}|\bm{y})$, is often used. However, we do not use a single value for the hyperparameters to receive benefit of GPs; various hyperparameters gives us a much wider class of sample functions. We do inference of their probability distributions by a sampling scheme to integrate out the hyperparameter dependence. In this case, the posterior distribution of $\bm{y}^{*}$ for unobserved input $\bm{x}^{*}$ is written as

[TABLE]

In addition to this, we extend this formulation to more general one to incorporate physics information into GPs. In many practical situations, the target function is known to have shape constraints (e.g. monotonicity or convexity) or inequality constraints. That is also the case with problems of interest, i.e. energy eigenvalues in FCI are monotonic and (almost) convex with respect to $N_{\mathrm{max}}$. In general, the accuracy of a statistical model like GP is improved by including such physics information. To this end, we extend Eq. (13) to

[TABLE]

where $P(\alpha,\beta,\ldots|\bm{y}^{*},\bm{y})$ is the probability that the constraints $\alpha,\beta,...$ are satisfied under the given $\bm{y}^{*}$ and $\bm{y}$. The contribution to the integral is determined by the balance among the posterior for the prediction, the likelihood for the hyperparameters, the hyperparameter prior, and the fidelity to the constraints. This expression is justified when the $\bm{\theta}$ is independent of the constraints (see Appendix B). This is rather general expression, i.e. constraints can be introduced independently for each problem of interest.

In general, the integration in Eq. (14) cannot be evaluated analytically. Therefore, some approximation or sampling scheme is required. We evaluate the integration in Eq. (14) by weighted $N_{p}$ samples as follows:

[TABLE]

We employ the particle filtering method [53] (also known as Sequential Monte Carlo) as a sampling scheme to evaluate the summation in Eq. (15). In our particle filtering algorithm, states $\{\bm{\theta}^{(i)}$, $\bm{y}^{*(i)}\}$ are assigned to particles labeled by $i=1,2,...,N_{p}$, and those particles are evolved independently according to the Metropolis-Hastings method with the so-called resampling scheme; at a certain step of the algorithm, the particles which do not respect the physics constraints are discarded.

III Problems of interest

III.1 FCI results

In what follows, we apply the constrained GP model to extrapolation problems in FCI calculations. We analyze published FCI results of the ground state energy of 6Li using JISP16/$\mathrm{NNLO}_{\mathrm{opt}}$ interaction with $\hbar\Omega=17.5$ MeV [54] and N3LO interaction with $\hbar\Omega=16.0$ MeV which is softened by similarity renormalization group (SRG) method with a flow parameter $\lambda=2.02$ fm*-1* [55]. The results are summarized in FIG. 2 as a function of $N_{\mathrm{max}}$.

Let $\{(x_{1},y_{1}),(x_{2},y_{2}),...,(x_{D},y_{D})|x_{1}<x_{2}<\cdots<x_{D}\}$ denote the data, i.e. $(x_{1},y_{1})=(6,-28.602)$, …, $(x_{D},y_{D})=(14,-31.977)$ in the case of N3LO results. Unlike least squares fitting of parametric models in which one should remove outliers from data, there is no reason to reduce data in the GP model and we use all $N_{\mathrm{max}}$ results as data unless otherwise mentioned.

The extrapolation problem addressed below is to estimate the ground state energies at $N_{\mathrm{max}}$ larger than $x_{D}$ and we express them as $\{(x^{*}_{1},y^{*}_{1}),(x^{*}_{2},y^{*}_{2}),...,(x^{*}_{P},y^{*}_{P})|x^{*}_{1}<x^{*}_{2}<\cdots<x^{*}_{P}\}$; Here $x^{*}_{1}=x_{D}+2$ and $P$ is large even integer. For the sake of simplicity, we restrict ourselves to consider the ground state of 6Li with natural parity, i.e. we only consider even $N_{\mathrm{max}}$. In practice, we truncate at certain finite $P$ value where predictions are converged with respect to $N_{\mathrm{max}}$. A detailed discussion about this $P$ will be given later.

III.2 Constraints on extrapolation function

As minimal constraints on the extrapolation function to capture the asymptotic behavior of FCI results, we impose the following two constraints, $\alpha$ and $\beta$.

The constraint $\alpha$ is variational property, i.e. the monotonicity of energy eigenvalues with respect to $N_{\mathrm{max}}$:

[TABLE]

Here we introduced the Probit function

[TABLE]

where $\kappa$ is the controlling parameter of the strictness of constraints. This probit function approaches the step function at $z=0$ when $\kappa\to\infty$. This $\kappa$ is set as about $10^{6}$ in our code, which is large enough to impose the constraints with satisfactory accuracy less than $0.01$ keV. More precisely, $\kappa$ is gradually increased in our code to $\sim 10^{6}$ so as to avoid possible localization at the early steps of the Monte Carlo sampling. We confirmed that the form of this $\kappa$ as a function of the Monte Carlo step does not affect the extrapolation results other than the sampling efficiency.

The second constraint $\beta$ is about the convergence pattern. We use a ratio of the absolute value of energy gains as a measure of convergence of g.s. energy in FCI. This ratio $r$ at a certain point $x_{j}$ is defined as follows:

[TABLE]

We plot $\{r(x_{j})\}$ for the given data in FIG. 2. For the energy eigenvalues by the FCI method, the denominator and numerator in Eq. (19) are both positive. If the calculated results of the g.s. energy exactly obey an exponential function, it means $r$ is a constant, which is not the case as shown in FIG. 2.

This $r$ can be extended to include predictions $\{y^{*}_{j}\}$:

[TABLE]

We impose the constraint on $\{r\}$ as follows:

[TABLE]

where $R_{E}$ is an upper threshold determined as follows:

[TABLE]

We use $\sigma_{r}=0.1$ throughout this work for simplicity. As can be expected from FIG. 2, this is a rather moderate constraint on the convergence pattern. When $\sigma_{r}$ is large enough, results agree with ones with only the constraint $\alpha$. We refer to the GP extrapolation model using constraints as the constrained Gaussian process (cGP) model.

III.3 Choice of the mean function

Here we introduce two different mean functions $\bm{\mu}^{(*)}$ needed in Eqs.(7)–(12):

• •

(case a) zero mean: $\bm{\mu}=\bm{0}_{D},\bm{\mu}^{*}=\bm{0}_{P}$

• •

(case b) B3 fit: mean of data and prediction are both determined by B3 fit [22], i.e. minimizing $\chi^{2}$ deviation between the largest three $N_{\mathrm{max}}$ data and the exponential function in the form of $E_{\infty}+c_{0}\exp{(-c_{1}N_{\mathrm{max}})}$ with three free parameters $(E_{\infty},c_{0},c_{1})$. In this choice, it can be said that preliminary knowledge on the behavior of the quantity is included in terms of the mean function of GPs.

We refer to these as cGP-a and cGP-b, respectively and analyze both cases below.

IV Results and Discussions

IV.1 Extrapolation of g.s. energies

We now present the results of the cGP predictions for the ground-state energies of 6Li in FIG. 3. The extrapolated values at a certain $N_{\mathrm{max}}$ are shown by histograms. The cGP-a results are shown by the transparent histograms colored in green, and the hatched histograms colored in pink are for the cGP-b results.

We note that the prediction is truncated up to a certain $N_{\mathrm{max}}$ where the mean value is converged within 0.2 keV. It means that possible differences between predictions at the finite $N_{\mathrm{max}}$ and one at $N_{\mathrm{max}}=\infty$ are suppressed below 1 keV because of the constraint $\beta$. The $N_{\mathrm{max}}$ giving converged results are 44,42 and 34 for NNLOopt, JISP16, and N3LO, respectively. These numbers are consistent with the intuition that harder interaction requires larger $N_{\mathrm{max}}$ to obtain converged results, while it should be noted that these results are for different $\hbar\Omega$. In addition to the mean values, the 68$\%$ and 99$\%$ credible intervals are shown and plotted below the histograms. In this manuscript, the 68$\%$ interval is determined from the 16th and 84th percentile of the distribution, and the $99\%$ interval is defined in a similar manner. As a whole, the cGP-b gives smaller credible intervals than the cGP-a.

Other extrapolated values are also shown. The B3 denotes the exponential fit using the largest three $N_{\mathrm{max}}$ data. The conventional B3 extrapolation is always included as a special case of two cGP results. The values for Literature are from Ref. [54, 55]. Here we note that the literature value for N3LO [55] might be obtained by an exponential fit using all five data, which is not explicitly stated. It must also be noted that it is a highly non-trivial task to fairly compare the results with different extrapolation techniques, because data is truncated in parametric models and some use data with multiple $\hbar\Omega$ as in the A5 extrapolation [54].

We note that our sampling scheme with the particle filtering gives converged results within a few keV in case of $20,000$ particles after 2,000 times Metropolis-Hastings updates for each particle, and that independent runs reproduce the same results within the Monte Carlo error.

IV.2 Impact of the constraints

Here we discuss how extrapolation results are influenced by the two constraints imposed. In FIG. 4, the impact of the constraints is shown. All symbols are the same as in FIG. 2, and the cGP-a and cGP-b results are summarized in the left and right regions, respectively. The bands correspond to the $68\%$ credible intervals of GP predictions with $\alpha$ and $\beta$ (red), with the only $\alpha$ (orange), and without constraints (green), respectively. Since the $68\%$ errors for the cGP-b prediction is an order of ten keV, the red lines in the right panel are very narrow.

As shown in the textbook [43], predictions of unconstrained zero-mean GP at points far from the data domain converge to zero with a fixed standard error. For this reason, we omit the case of cGP-a without $\alpha$ and $\beta$, which is obviously not appropriate for the current purpose. When one assumes that the wave function is dominated by relatively lower $N_{\mathrm{max}}$ configurations, predictions with both constraints $\alpha$ and $\beta$ are expected to be more reliable than the others.

IV.3 Data dependence

We have used all $N_{\mathrm{max}}$ results as data so far. Here we explore the dependence of extrapolated values on the used data to test the potential predictive power of our extrapolation method. In FIG. 5, extrapolated values for both cGP-a and cGP-b are shown as a function of the maximum $N_{\mathrm{max}}$ used as data, i.e. $x_{D}=\mathrm{max}(N_{\mathrm{max}})$.

The mean predictions by the cGP models are shown by horizontal line in white, and the shaded areas show the posterior distributions obtained by 20,000 particles. The 68$\%$ and 99$\%$ credible intervals, respectively, are shown by the thick and thin vertical lines, respectively. We note that the width of the shaded area is scaled to be the same for each area, and $\mathrm{NNLO}_{\mathrm{opt}}$ (N3LO) results are slightly shifted to left (right) for visibility.

As a whole, the size of credible intervals for the cGP-a is larger than that for the cGP-b, and the credible intervals become smaller as higher $N_{\mathrm{max}}$ data is added with only one exception, which is the cGP-b result for N3LO with max($N_{\mathrm{max}}=10$). This exception can be understood from Eqs. (8)–(12). In this case, the exponential function exactly fits the given three data and then $\bm{\mu}_{\bm{y}^{*}|\bm{y}}$ in Eq. (10) is identical with $\bm{\mu}^{*}$. Any fluctuation of the joint mean value $\bm{\mu}_{\bm{y}^{*}|\bm{y}}$ is not allowed, and this significantly reduces probability weights for functions other than the B3 fit.

These plots with quantified uncertainties tell us one criterion of where to stop the massive FCI calculations under the given extrapolation model; It is inadvisable to carry out FCI calculations while increasing $N_{\mathrm{max}}$ forever, whereas the point to stop must depend on the desired accuracy.

In the rest of this subsection, let us regard the mean values at rightmost max($N_{\mathrm{max}}$) in FIG. 5, i.e. max($N_{\mathrm{max}}$)=18 for $\mathrm{NNLO}_{\mathrm{opt}}$ and JISP16, and max($N_{\mathrm{max}}$)=14 for N3LO, as the (temporary) exact values.

For NNLOopt and JISP16, an important remark is that the exact values are covered by the cGP-a predictions with a relatively lower max($N_{\mathrm{max}}$). The sign of convergence can be seen around max($N_{\mathrm{max}}$)=10-12, i.e. one can choose these as points to stop the calculation. From the behavior of the credible intervals, the cGP-b seems to underestimate the uncertainties than the cGP-a. In other words, the cGP-b takes account of fluctuation of the functional form only around the exponential function, while the cGP-a would include a wider class of functions.

For N3LO, any predictions by cGP-a, cGP-b, and B3 fit with lower max($N_{\mathrm{max}}$) have almost no overlap with the tentative exact value. It means that all the models fail to estimate the extrapolation uncertainty, while the cGP-a could be better than the others. The extrapolated values are much more sensitive to the max($N_{\mathrm{max}}$) than the results with other interactions; the extrapolation for the N3LO results seems to be a more non-trivial problem than the others. This can be understood from the behavior of the ratio of energy gains $r$. Especially in the N3LO case, the $r$ is unstable with respect to $N_{\mathrm{max}}$, as seen in FIG. 2. There are at least two possible explanations for this non-flat behavior of $r$. One possibility is that the calculation have not yet converged, i.e. the additional bindings by increasing $N_{\mathrm{max}}$ cannot be regarded as a simple asymptotic behavior. The other one is the SRG evolution of the input nuclear potential. There is not yet enough open data to conclude the origin of the non-flat behavior of $r$. If one could figure out an additional constraint on the behavior of extrapolated values as a function of max($N_{\mathrm{max}}$), such a difficulty in extrapolation, which could be observed in particular nuclei and/or particular interactions, would be alleviated.

V Summary and Outlook

We introduced an extrapolation method for CI-type calculations using constrained Gaussian processes. This method has the following advantages that are required for future generations of ab initio studies to make more quantitative discussions on observables of interest and on the quality of adopted nuclear interactions.

Firstly, our extrapolation method does not need to remove outliers and has applicability to sparse data sets, which are strongly needed for future FCI calculations. Secondly, one can naturally incorporate domain knowledge into the model. It is often the case especially in physics that one knows in advance behavior of the target quantity at a certain level, which is ranging from empirical laws to physical principles. One can expect that imposing such information into the extrapolation model improves the accuracy of the predictions. This flexibility might be useful to alleviate difficulties in the extrapolation for some particular cases. Thirdly, uncertainty in extrapolation can be quantified in a systematic manner. Although the main source of uncertainties in FCI calculations is the input nuclear potential, evaluating extrapolation uncertainties are helpful for further understandings about the nuclear observables.

Regarding uncertainties from input parameters in nuclear many-body methods, the tremendous efforts to propagate input uncertainty to the observables have been made in the last decade, see e.g. [56, 57, 58, 59, 60, 61, 62, 63, *UQpot_1e, 65, 66, 67, 68, 69, 70, 71]. In addition to these, it has been shown that eigenvector continuation (EC), which is introduced in Ref. [72], can be used as an efficient emulator of ab initio methods, then used for uncertainty quantification and sensitivity analyses on input parameters such as low-energy constants in the chiral EFT potentials [73, 74]. It is expected that EC (and some other method) facilitates comprehensive studies of uncertainty propagation in ab initio methods combined with an extrapolation method with quantified uncertainties.

The benefits of the uncertainty quantification are not limited to putting error bars in predictions. If one properly propagate uncertainties from the input interaction and also quantifies uncertainties such as that due to the extrapolation, it enables us to visualize non-linear relation between input and output of many-body calculations and capability of the many-body method. Then it would provide us with footholds to understand some missing contributions, if there were any.

In this work, we discussed only the ground state energies obtained by FCI calculations. When it comes to the extrapolation problem of other quantities or in other systems, the main problem is to find minimal constraints to capture the asymptotic behavior of the quantities well. It is a possible future direction along this line to extend the cGP model to a higher dimension. In case of FCI calculations, for example, one can impose the following additional constraint on GP by extending the formulation to $(N_{\mathrm{max}},\hbar\Omega)$ space: extrapolated values with different $\hbar\Omega$ should converge to the same value to some extent. The extension of the formulation to a multi-dimensional space is rather straightforward, while it is expected that one needs more technical analyses in numerical studies such as positive semi-definiteness of covariance matrices. Our model can also be applied to valence CI techniques using an importance-truncation in which extrapolation function is much more non-trivial than the FCI case. It would also be interesting to apply this kind of cGP to finite-size scaling analyses in other systems.

Acknowledgements.

The author thanks Noritaka Shimizu for fruitful discussions on future applications to the extrapolation problem in Monte Carlo shell model calculations. This work was supported by JSPS KAKENHI (Grants No. 17J06775).

Appendix A Kernel selection

In this section, we discuss the technical details to choose the kernel function. As mentioned in the main text, popular choices are RBF kernel:

[TABLE]

and Matérn kernel:

[TABLE]

where $\Gamma$ is the gamma function and $K_{\nu}$ is the modified Bessel function of the second kind. Especially, Matérn kernel with $\nu=3/2$ and $\nu=5/2$ are commonly used. The RBF kernel corresponds to the special case of the Matérn kernel with $\nu=\infty$. For the Matérn kernel, a sample function is $k$-times mean square differentiable if and only if $\nu>k$ [43]. For that reason, the Matérn kernels with $\nu>3/2$ are thought to be appropriate for our purpose, i.e. extrapolation of FCI results, but we include $\nu=3/2$ case too.

In what follows, we consider these Mat32 (Matérn kernel with $\nu=3/2$), Mat52 (Matérn kernel with $\nu=5/2$), RBF (Radial basis function), and their counterparts with logarithm distances:

[TABLE]

where $\eta\equiv|\ln x_{i}-\ln x_{j}|$. Taking the logarithm distance, i.e. replacing $|x_{i}-x_{j}|$ in e.g. Eq. (25) by $|\ln{x_{i}}-\ln{x_{j}}|$, makes results independent on the scale of the x-axis, and allows to capture the non-stationary nature of FCI results (results rapidly converge to certain values as functions of $N_{\mathrm{max}}$) [75].

To compare the six kernels, we here use the Bayesian information criterion (BIC) defined as follows:

[TABLE]

where $\bm{\theta}_{\mathrm{ML}}$ is the maximum likelihood estimation under the given data, $M$ is the dimension of the hyperparameter $\bm{\theta}$, and $n$ is the number of data. In this case, the values of the latter term are common among the results with the given interaction. In FIG. 6, we summarize the BIC for the six kernels. Colors and symbols for three interactions are the same as FIG. 2. The $\bm{\theta}_{\mathrm{ML}}$ is obtained through updates of 2,000 independent particles by the Metropolis-Hastings method, and the $\bm{\theta}_{\mathrm{ML}}$ is converged less than $0.1\%$ accuracy. As a whole, taking the logarithm of the distance increases the maximum loglikelihood, which is now equivalent to the maximum BIC.

As stated in the main text, we concluded to use the logMat52 kernel. This is partly because the logMat52 gives the highest BIC in total. The other reason is its numerical stability over the (log)RBF kernel. The RBF is one of the most popular choices for the kernel. However, its smoothness of the sample functions is often regarded as too high [76], and, in practice, this too smooth nature sometimes breaks down the positive semi definiteness of covariant matrices in numerical calculations due to rounding errors. We briefly demonstrate that in the following.

In applications of GPs, one needs to obtain conditional mean vectors Eq. (10) and covariance matrices Eq. (11). The covariance matrices such as $K_{XX}$ and $\bm{\Sigma}$ in Eq. (3) must be positive semi-definite (PSD) to achieve e.g. the Cholesky decomposition for $K^{-1}_{XX}$ and to gemerate samples from the posterior distribution.

In some cases (e.g. points are located too close to each other), covariance matrices become non-PSD due to rounding errors, although it must be PSD mathematically. This is also true for our extrapolation method. The typical prescription to this non-PSD is to add an infinitesimal diagonal matrix to the $K_{XX}$ and/or $\bm{\Sigma}$. Let us call this the $\epsilon$ prescription in the following. This $\epsilon$ prescription is mathematically equivalent to assume the observation and/or prediction to have noise.

In the current case the target quantity is the calculated energies of light nuclei by ab initio full configuration interaction method, i.e. in the order of a few tens of MeV. On the other hand, the typical convergence tolerance of the Lanczos method in shell-model codes on the market is 1.e-5 MeV or better [18, *KSHELL2, 20], i.e. the problems of interest are almost noiseless.

We calculated the conditional mean vectors and covariance matrices while varying the hyperparameters $\tau$ and $\ell$ in Eq. (28)–Eq. (30). In Fig. 7, the positive semi-definiteness of the covariance matrices are summarized. The diamond symbols (red) correspond to the case that both $K_{XX}$ and $\Sigma_{\bm{y}^{*}|\bm{y}}$ are non-PSD, and the cross symbols (blue) are assigned if only the $\Sigma_{\bm{y}^{*}|\bm{y}}$ is non-PSD. In terms of the length scale $\ell$, the logRBF gives non-PSD matrices easier than logMat52 by one or two orders of magnitude.

Now we show the effects of the $\epsilon$ prescription on the predictions. In FIG. 8, the following quantity is shown:

[TABLE]

where the superscripts own and lib. mean the conditional mean vector calculated by our own code and one using the GaussianProcesses.jl library [77], respectively. For the latter, we fixed the size of observation noise as zero and allowed to use the default $\epsilon$ prescription implemented in the library. We note that both $\bm{\mu}^{\mathrm{own}}$ and $\bm{\mu}^{\mathrm{lib.}}$ are now given in a unit of MeV. The hatched regions with dots for smaller $\ell$ means that the results of the two codes are exactly the same. At the points shown by the diamond symbol in FIG. 7, our code without the $\epsilon$ prescription cannot give the mean vectors. Those correspond to the white region appeared in FIG. 8. When we increase the length scale $\ell$ at which the posterior covariance becomes non-PSD, i.e. points shown by the cross symbol in FIG. 7, the $\Delta$ becomes larger because of the immune system in the library. In some cases, the deviations in mean vectors reach a few MeV, which are obviously non-negligible.

One should pay much attention to how the PSD and the observation noise are treated in the codes and whether or not one really can neglect the impact of the $\epsilon$ prescriptions on the predictions, especially when one would like to integrate out the hyperparameter dependence. For more detailed analyses and codes to reproduce the FIG. 7 and FIG. 8, we refer the reader to Ref. [78].

From the analyses above, we concluded that the logMat52 is the most appropriate choice for our extrapolation method.

Appendix B Extension of the posterior distribution under the constraints

Eq. (14) is derived in the followings. Let the $\bm{c}$ denote the physics constraints to be imposed.

[TABLE]

Appendix C Contributions to the posterior distribution

In Figures 10–12, the hyperparameter distributions are shown. Colors of each dot point show the relative size of contributions to the likelihood $P(\bm{\theta}|\bm{y})$ and to the integral in Eq. (14). In our sampling method, particles are evolved according to the random walk Metropolis-Hastings method. We can see from FIG. 10 and FIG. 12 that particles distribute like an ellipse in hyperparameter space. The particles do not distribute over the non-PSD region discussed in Sec. B.