© G. Schnabel, published by EDP Sciences, 2018

1 Introduction

Despite theoretical advances, nuclear models are in general not able to reproduce all features of trustworthy experimental data. Because experiments alone do not provide sufficient information to solve problems of nuclear engineering, models are still needed to fill the gaps.

Being in need of reliable nuclear data, a pragmatic solution to deal with imperfect models is the introduction of a bias term on top of the model. The true prediction is then given as the sum of model prediction and bias term. The form of the bias term is in principle arbitrary and a low order polynomial could be a possible choice. However, assumptions about how models deviate from reality are usually very vague, which makes it difficult to justify one parametrization over another one.

Methods of non-parametric statistics help to a certain extent to overcome this specification problem. In particular, Gaussian process (GP) regression (also known as Kriging, e.g., [1]) enjoys popularity in various fields, such as geostatistics, remote sensing and robotics, due to its conceptual simplicity and sound embedding in Bayesian statistics. Instead of providing a parameterized function to fit, one specifies a mean function and a covariance function. The definition of these two quantities induces a prior probability distribution on a function space. Several standard specifications of parametrized covariance functions exist, e.g., [[1], Chap. 4], whose parameters regulate the smoothness and the magnitude of variation of the function to be determined by GP regression.

The optimal choice of the values for these so-called hyperparameters is problem dependent and can also be automatically performed by methods such as marginal likelihood optimization and cross validation, e.g., [[1], Chap. 5]. In addition, features of various covariance functions can be combined by summing and multiplying them. Both the automatic determination of hyperparameters and the combination of covariance functions will be demonstrated.

From an abstract viewpoint, GP regression is a method to learn a functional relationship based on samples of input-output associations without the need to specify a functional shape. GP regression naturally yields besides estimates also the associated uncertainties. This feature is essential for evaluating nuclear data, because uncertainties of estimates are as important for the design of nuclear facilities as are the estimates themselves. Prior knowledge about the smoothness and the magnitude of the model bias can be taken into account. Furthermore, many existing nuclear data evaluation methods, e.g., [2–4], can be regarded as special cases. This suggests that the approach discussed in this paper can be combined with existing evaluation methods in a principled way.

The main hurdle for applying GP regression is the bad scalability in the number of data points N. The required inversion of an N × N covariance matrix leads to computational complexity of N³. This limits the application of standard GP regression to several thousand data points on contemporary desktop computers. Scaling up GP regression to large datasets is therefore a field of active research. Approaches usually rely on a combination of parallel computing and the introduction of sparsity in the covariance matrix, which means to replace the original covariance matrix by a low rank approximation, see e.g., [5].

In this paper, I investigate the sparse approximation introduced in [6] to estimate the model bias of inclusive double-differential neutron spectra over the complete nuclide chart for incident protons above 100 MeV. The predictions are computed by the C++ version of the Liège intranuclear cascade model (INCL) [7,8] coupled to the Fortran version of the evaporation code ABLA07 [9]. The experimental data are taken from the EXFOR database [10].

The idea of using Gaussian processes to capture deficiencies of a model exists for a long time in the literature, see e.g., [11], and has also already been studied in the context of nuclear data evaluation, e.g., [12–14]. The novelty of this contribution is the application of GP regression to a large dataset with isotopes across the nuclide chart, which is possible thanks to the sparse approximation. Furthermore, the way GP regression is applied enables predictions for isotopes without any data.

The exemplary application of sparse GP regression in this paper indicates that the inclusion of hundred thousands of data points may be feasible and isotope extrapolations yield reasonable results if some conditions are met. Therefore, sparse GP regression is a promising candidate to perform global assessments of models and to quantify their imperfections in a principled way.

The structure of this paper is as follows. Section 2 outlines the theory underlying sparse GP regression. In particular, Section 2.1 provides a succinct exposition of standard GP regression, Section 2 sketches how to construct a sparse approximation to a GP and how this approximation is exploited in the computation, and Section 2.3 explains the principle to adjust the hyperparameters of the covariance function based on the data.

The application to INCL/ABLA and inclusive double-differential neutron spectra is then discussed in Section 3. After a brief introduction of INCL and ABLA in Section 3.1, the specific choice of covariance function is detailed in Section 3.2. Some details about the hyperparameter adjustment are given in Section 3.3 and the results of GP regression are shown and discussed in Section 3.4.

2 Method

2.1 GP regression

GP regression, e.g., [1], can be derived in the framework of Bayesian statistics under the assumption that all probability distributions are multivariate normal. Let the vector y₁ contain the values at the locations $\mathrm{}\{x_i^p\}\mathrm{}$ of interest and y₂ the observed values at the locations $\mathrm{}\{x_j^o\}\mathrm{}$.

For instance, in the application in Section 3, the elements in the vector y₁ represent the relative deviations of the “truth” from the model predictions for neutron spectra at angles and energies of interest. The vector y₂ contains the relative deviations of the available experimental data from the model predictions for neutron spectra at the angles and energies of the experiments. The underlying assumption is that the experimental measurements differ from the truth by an amount compatible with their associated uncertainties. If this assumption does not hold, model bias should be rather understood as a combination of model and experimental bias.

Given a probabilistic relationship between the vectors y₁ and y₂, i.e. we know the conditional probability density function (pdf) of y₂ given y₁, ρ(y₂|y₁) (e.g., because of continuity assumptions), the application of the Bayesian update formula yields $$\rho \mathrm{}(y_1|\mathrm{}{y}_2\mathrm{})\propto \rho \mathrm{}(y_2\mathrm{}|y_1)\mathrm{}\rho \mathrm{}(y_1)\mathrm{}\mathrm{}\text{.}$$(1)

The occurring pdfs are referred to as posterior $\rho \mathrm{}(y_1|\mathrm{}{y}_2\mathrm{})$, likelihood $\rho \mathrm{}(y_2\mathrm{}|y_1)\mathrm{}$, and prior $\rho \mathrm{}(y_1)\mathrm{}$. The posterior pdf represents an improved state of knowledge.

The form of the pdfs in the Bayesian update formula can be derived from the joint distribution $\rho \mathrm{}(y_1,\mathrm{}{y}_2)\mathrm{}$. In the following, we need the multivariate normal distribution $$\mathcal{N}(y|\mu ,K)=\frac{1}{\sqrt{{(2\pi )}^N\det K}}\times \exp \left(-\frac{1}{2}{(y-\mu )}^T{K}^{-1}(y-\mu )\right)\text{,}$$(2) which is characterized by the center vector μ and the covariance matrix K. The dimension of the occurring vectors is denoted by N. Under the assumption that all pdfs are multivariate normal and centered at zero, the joint distribution of y₁ and y₂ can be written as $$\rho (y_1,y_2)=\mathcal{N}\left(\left(\begin{array}{c}\hfill y_1\hfill \\ \hfill y_2\hfill \end{array}\right)|\left(\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \end{array}\right),\left(\begin{array}{cc}\hfill K_{11}\hfill & \hfill K_{21}^T\hfill \\ \hfill K_{21}\hfill & \hfill K_{22}\hfill \end{array}\right)\right)\text{.}$$(3)

The compound covariance matrix contains the blocks K₁₁ and K₂₂ associated with $y_1$ and y₂, respectively. The block K₂₁ contains the covariances between the elements of y₁ and y₂. Centering the multivariate normal pdf at zero is a reasonable choice for the estimation of model bias. It means that an unbiased model is a priori regarded as the most likely option.

The posterior pdf is related to the joint distribution by $$\rho \mathrm{}(y_1|\mathrm{}{y}_2)=\mathrm{}\frac{\rho \mathrm{}(y_1\mathrm{},y_2)\mathrm{}}{\int \rho \mathrm{}(y_1\mathrm{},y_2)\mathrm{}d{y}_1}\mathrm{}\text{.}$$(4)

The solution for a multivariate normal pdf is another multivariate normal pdf. For equation (3), the result is given by $$\rho (y_1|y_2)=\mathcal{N}(y_1|y{\text{'}}_1,K{\text{'}}_{11})\text{,}$$(5) with the posterior mean vector $y{\text{'}}_1$ and covariance matrix K'₁₁ (e.g., [[1], A.3]), $$y{\text{'}}_1\mathrm{}=K_{21}^T{K}_{22}^{-1}{y}_2\mathrm{}$$(6) $$K{\text{'}}_{11}=K_{11}-K_{21}^T{K}_{22}^{-1}{K}_{21}\text{.}$$(7)

The important property of these equations is the fact that the posterior moments depend only on the observed vector y₂. The introduction of new elements into y₁ and associated columns and rows into K₁₁ in equation (3) has no impact on the already existing values in $\mathit{y}{\text{'}}_1$ and K'₁₁. In other words, one is not obliged to calculate all posterior expectations and covariances at once. They can be calculated sequentially.

GP regression is a method to learn a functional relationship $f(x)$ based on samples of input-output associations {(x₁, f₁), (x₂, f₂), … }. The vector y₂ introduced above then contains the observed functions values of $f\mathrm{}(x)\mathrm{}$, i.e. $y_2\mathrm{}=(f_1,\mathrm{}{f}_2\mathrm{}\dots ){\mathrm{}}^T$. Assuming all prior expectations to be zero, the missing information to evaluate equations (6) and (7) are the covariance matrices. Because predictions for $f(x)\mathrm{}$ should be computable at all possible locations x and the same applies to observations, covariances between functions values must be available for all possible pairs $\mathrm{}(x_i\mathrm{},x_j\mathrm{})$ of locations. This requirement can be met by the introduction of a so-called covariance function $\kappa \mathrm{}(x_i,\mathrm{}{x}_j\mathrm{})$. A popular choice is the squared exponential covariance function $$\kappa \mathrm{}(x_i\mathrm{}{x}_j)\mathrm{}={\delta }^2\exp \left(-\frac{{\mathrm{}(x_i-x_j)\mathrm{}}^2}{2{\lambda }^2}\right)\mathrm{}\text{.}$$(8)

The parameter δ enables the incorporation of prior knowledge about the range functions values are expected to span. The parameter λ regulates the smoothness of the solution. The larger λ the slower the covariance function decays for increasing distance between x₁ and x₂ and consequently the more similar function values are at nearby locations.

If the set $\mathrm{}\{x_i^p\}\mathrm{}$ contains the locations of interest and $\mathrm{}\{x_j^o\}\mathrm{}$ the observed locations, the required covariance matrices in equations (6) and (7) are given by $$K_{22}=\left(\begin{array}{ccc}\hfill \kappa (x_1^o,x_1^o)\hfill & \hfill \kappa (x_1^o,x_2^o)\hfill & \hfill \cdots \hfill \\ \hfill \kappa (x_2^o,x_1^o)\hfill & \hfill \kappa (x_2^o,x_2^o)\hfill & \hfill \cdots \hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill \end{array}\right)$$(9) and $$K_{21}=\mathrm{}\left(\begin{array}{ccc}\hfill \kappa \mathrm{}(x_1^o,\mathrm{}{x}_1^p)\mathrm{}\hfill & \hfill \kappa \mathrm{}(x_1^o,\mathrm{}{x}_2^p)\mathrm{}\hfill & \hfill \cdots \hfill \\ \hfill \kappa \mathrm{}(x_2^o,\mathrm{}{x}_1^p)\mathrm{}\hfill & \hfill \kappa \mathrm{}(x_2^o\mathrm{},x_2^p)\mathrm{}\hfill & \hfill \cdots \hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill \end{array}\right)\mathrm{}\text{.}$$(10)

2.2 Sparse Gaussian processes

The main hurdle for applying GP regression using many observed pairs (x_i, f_i) is the inversion of the covariance matrix k₂₂ in equations (6) and (7). The time to invert this N × N matrix with N being the number of observations increases proportional to N³. This limits the application of standard GP regression to several thousand observations on contemporary desktop computers. Parallelization helps to a certain extent to push this limit. Another measure is the approximation of the covariance matrix by a low rank approximation. I adopted the sparse approximation described in [6], which will be briefly outlined here. For a survey of different approaches and their connections consult e.g., [5].

Suppose that we have not measured the values in y₂ associated with the locations $\mathrm{}\{x_j^o\}\mathrm{}$, but instead a vector y₃ associated with some other locations $\mathrm{}\{x_k^{psi}\}\mathrm{}$. We refer to y₃ as vector of pseudo-inputs. Now we use equation (6) to determine the hypothetical posterior expectation of y₂, $$y{\text{'}}_2\mathrm{}={\displaystyle \stackrel{S}{{\displaystyle \overbrace{K_{23}{K}_{33}^{-1}}}}}{y}_3\mathrm{}\text{.}$$(11)

The matrices K₂₃ and K₃₃ are constructed analogous to equations (10) and (9), respectively, i.e. ${(K_{23})}_{ij}=\kappa (x_i^o,x_j^{psi})$ and ${(K_{33})}_{jk}=\kappa (x_j^{psi},x_k^{psi})$ with i = 1..N, j = 1..M and N being the number of observations and M the number of pseudo-inputs. Noteworthy, y'₂ is a linear function of y₃. Under the assumption of a deterministic relationship, we can replace the posterior expectation y'₂ by y₂. Using the sandwich formula and $\rho (y_3)=\mathcal{N}(y_3|0,K_{33})$, we get for the covariance matrix of y₂ $${{\displaystyle \tilde{K}}}_{22}=S{K}_{33}{S}^T=K_{23}{K}_{33}^{-1}{K}_{23}^T\text{.}$$(12)

Given that both K₃₃ and K₂₃ have full rank and the number of observations N is bigger than the number of pseudo-inputs M, the rank of ${{\displaystyle \tilde{K}}}_{22}$ equals M. The approximation is more rigid than the original covariance matrix due to the lower rank.

In order to restore the flexibility of the original covariance matrix, the diagonal matrix $D_2=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[K_{22}-{{\displaystyle \tilde{K}}}_{22}]$ is added to ${{\displaystyle \tilde{K}}}_{22}$. This correction is essential for the determination of the pseudo-input locations via marginal likelihood optimization as explained in [[6], Sect. 3]. Furthermore, to make the approximation exhibit all properties of a GP, it is also necessary to add $D_1=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[K_{11}-{{\displaystyle \tilde{K}}}_{11}]$ to ${{\displaystyle \tilde{K}}}_{11}=K_{31}^T{K}_{33}^{-1}{K}_{31}$ as explained in [[5], Sect. 6].

Making the replacements $K_{22}\to {{\displaystyle \tilde{K}}}_{22}+D_2$, $K_{21}\to {{\displaystyle \tilde{K}}}_{21}=K_{23}{K}_{33}^{-1}{K}_{31}$, and $K_{11}\to {{\displaystyle \tilde{K}}}_{11}+D_1$ in equations (6) and (7), we obtain $$y_{1^{\prime }}={{\displaystyle \tilde{K}}}_{21}^T{\left({{\displaystyle \tilde{K}}}_{22}+D_2\right)}^{-1}{y}_2\text{,}$$(13) $$K{\text{'}}_{11}={{\displaystyle \tilde{K}}}_{11}+D_1-{{\displaystyle \tilde{K}}}_{21}^T{\left({{\displaystyle \tilde{K}}}_{22}+D_2\right)}^{-1}{{\displaystyle \tilde{K}}}_{21}\text{.}$$(14)

Using the Woodbury matrix identity (e.g., [1, A. 3]), these formulas can be rewritten as (e.g.,[5, Sect. 6]), $$y_{1^{\prime }}=K_{31}^T{\left(K_{33}+K_{23}^T{D}_2^{-1}{K}_{23}\right)}^{-1}{K}_{23}^T{D}^{-1}{y}_2\text{,}$$(15) $$K{\text{'}}_{11}=D_1+K_{31}^T{\left(K_{33}+K_{23}^T{D}_2^{-1}{K}_{23}\right)}^{-1}{K}_{31}\text{.}$$(16)

Noteworthy, the inversion in these expressions needs only to be performed for an M × M matrix where M is the number of pseudo-inputs. Typically, M is chosen to be in the order of magnitude of hundred. The computational cost for inverting the diagonal matrix D is negligible.

Using equations (15) and (16), the computational complexity scales linearly with the number of observations N due to the multiplication by K₂₃. This feature enables to increase the number of observations by one or two orders of magnitude compared to equations (6) and (7) in typical scenarios.

2.3 Marginal likelihood maximization

Covariance functions depend on parameters (called hyperparamters) whose values must be specified. In a full Bayesian treatment, the choice of values should not be informed by the same observations that enter into the GP regression afterwards. In practice, this ideal is often difficult to achieve due to the scarcity of available data and therefore frequently abandoned. A full Bayesian treatment may also become computationally intractable if there are too much data.

Two popular approaches to determine the hyperparameters based on the data are marginal likelihood maximization and cross validation, see e.g., [1, Chap. 5]. A general statement which approach performs better cannot be made. I decided to use marginal likelihood optimization because it can probably be easier interfaced with existing nuclear data evaluation methods.

Given a covariance function depending on some hyperparameters, e.g., κ(x_i, x_j2|δ, λ) with hyperparameters δ and λ as in equation (8), the idea of marginal likelihood maximization is to select values for the hyperparameters that maximize the probability density for the observation vector ρ(y₂). In the case of the multivariate normal pdf in equation (3), it is given by (e.g., [1, Sect. 5.4.1) $$\ln \rho (y_2)=-\frac{N}{2}\log (2\pi )-\frac{1}{2}\log \det K_{22}-\frac{1}{2}{y}_2^T{K}_{22}^{-1}{y}_2\text{.}$$(17)

The first term is a constant, the second term is up to a constant the information entropy of the multivariate normal distribution, and the third term is the generalized χ²-value. The maximization of this expression amounts to balancing two objectives: minimizing the information entropy and maximizing the χ²-value.

The partial derivative of equation (17) with respect to a hyperparameter is given by (e.g., [[1], Sect. 5.4.1) $$\frac{\partial \ln \rho (y_2)}{\partial \lambda }=-\frac{1}{2}Tr\left(K_{22}^{-1}\frac{\partial K_{22}}{\partial \lambda }\right)+\frac{1}{2}{y}_2^T{K}_{22}^{-1}\frac{\partial K_{22}}{\partial \lambda }{K}_{22}^{-1}{y}_2\text{,}$$(18) which enables the usage of gradient-based optimization algorithms. In this paper, I use the L-BFGS-B algorithm [15] because it can deal with a large number of parameters and allows to impose restrictions on their ranges.

Due to the appearance of the determinant and the inverse of K₂₂, the optimization is limited to several thousand observations on contemporary desktop computers. However, replacing K₂₂ by the approximation ${{\displaystyle \tilde{K}}}_{22}+D_2$ in equations (17) and (18) enables to scale up the number of observations by one or two orders of magnitude. The structure of the approximation is exploited by making use of the matrix determinant lemma, the Woodbury identity, and the trace being invariant under cyclic permutations.

The approximation to K₂₂ is not only determined by the hyperparameters but also by the location of the pseudo-inputs $\left\{x_k^{psi}\right\}$. Hyperparameters and pseudo-input locations can be jointly adjusted by marginal likelihood maximization. The number of pseudo-inputs is usually significantly larger (e.g., hundreds) than the number of hyperparameters (e.g., dozens). In addition, the pseudo-inputs are points in a potentially multi-dimensional space and their specification requires a coordinate value for each axis. For instance, in Section 3 sparse GP regression is performed in a five dimensional space with three hundred pseudo-inputs, which gives 1500 associated parameters. Because equation (18) has to be evaluated for each parameter in each iteration of the optimization algorithm, its efficient computation is important.

The mathematical details are technical and tedious and hence only the key ingredient for efficient computation will be discussed. Let x_kl be the lth coordinate of the kth pseudo-input. The crucial observation is that ∂K₃₃/∂ x_kl yields a matrix in which only the lth and kth column and row contain non-zero elements. A similar statement holds for ∂K₂₃/∂ x_kl. This feature can be exploited in the multiplications and the trace computation in equation (18) (where K₂₂ is substituted by ${{\displaystyle \tilde{K}}}_{22}+D_2$) to achieve $\mathcal{O}(NM)$ per coordinate of a pseudo-input with M being the number of pseudo-inputs, and N being the number of observations. This is much more efficient than $\mathcal{O}(dN{M}^2)$ for the partial derivative with respect to a hyperparameter.

3 Application

3.1 Scenario

The model bias was determined for the C++ version of the Liège INCL [7], a Monte Carlo code, coupled to the evaporation code ABLA07 [9] because this model combination performs very well according to an IAEA benchmark of spallation models [16,17] and is used in transport codes such as MCNPX and GEANT4. This suggests that the dissemination of a more quantitative performance assessment of INCL coupled to ABLA potentially helps many people to make better informed decisions. Because some model ingredients in INCL are based on views of classical physics (as opposed to quantum physics), the model is mainly used for high-energy reactions above 100 MeV.

The ability of a model to accurately predict the production of neutrons and their kinematic properties may be regarded as one of the most essential features for nuclear engineering applications. Especially for the development of the innovative research reactor MYRRHA [18] driven by a proton accelerator, these quantities need to be well predicted for incident protons.

For these reasons, I applied the approach to determine the model bias in the prediction of inclusive double-differential neutron spectra for incident protons and included almost all nuclei for which I found data in the EXFOR database [10]. Roughly ten thousand data points were taken into account. Table 1 gives an overview of the data.

3.2 Design of the covariance function

The covariance function presented in equation (8) is probably too restrictive to be directly used on double-differential spectra. It incorporates the assumption that the model bias spans about the same range for low and high emission energies. Because the neutron spectrum quickly declines by orders of magnitude with increasing emission energy, it is reasonable to use a covariance function with more flexibility to disentangle the systematics of the model bias associated with these two energy domains.

Assuming the two covariance functions κ₁(x_i, x_j) and κ₂(x_i, x_j), a more flexible covariance function can be constructed in the following ways (e.g., [1, Sect. 4.2.4]) $${\kappa }_{1+2}(x_i,x_j)={\kappa }_1(x_i,x_j)+{\kappa }_2(x_i,x_j)\text{,}$$(19) $${\kappa }_{1\times 2}(x_i,x_j)={\kappa }_1(x_i,x_j){\kappa }_2(x_i,x_j)\text{.}$$(20)

Taking into account that the purpose of a covariance function is to compute elements of a covariance matrix, the construction in equation (19) is analogous to the possible construction of an experimental covariance matrix: An experimental covariance matrix can be assembled by adding a diagonal covariance matrix reflecting statistical uncertainties to another one with systematic error contributions. Please note that sparse GP regression as presented in this paper works only with a diagonal experimental covariance matrix. Equation (20) will be used to achieve a transition between the low and high energy domains.

To state the full covariance function used in this paper, the following covariance function on a one-dimensional input space is introduced, $$\kappa (x,x^{\prime }|\lambda )=\exp \left(-\frac{{(x-x^{\prime })}^2}{2{\lambda }^2}\right)\text{.}$$(21)

The meaning of the hyperparameter λ was explained below equation (8).

The inclusive double-differential neutron spectra for incident protons over the complete nuclide chart can be thought of as a function of spectrum values associated with points in a five dimensional input space. The coordinate axes are incident energy (En), mass number (A), charge number (Z), emission angle (θ) and emission energy (E). Using equation (21) we define two covariance functions κ_L and κ_H associated with the low and high energy domain of the emitted neutrons. Given two input vectors $$x_i={(E{n}_i,A_i,Z_i,{\theta }_i,E_i)}^T\text{,}$$(22) $$x_j={(E{n}_j,A_j,Z_j,{\theta }_j,E_j)}^T\text{,}$$(23) the form of the covariance function for both κ_L and κ_H is assumed to be $${\kappa }_{\mathit{C}}(x_i,x_j|{\delta }_{\mathit{C}},{\lambda }_{\mathit{C}})={\delta }_{\mathit{C}}^2\kappa (E{n}_i,E{n}_j|{\lambda }_{En}^C)\kappa (A_i,A_j|{\lambda }_A^{\mathit{C}})\kappa (Z_i,Z_j|{\lambda }_Z^{\mathit{C}})\kappa ({\theta }_i,{\theta }_j|{\lambda }_{\theta }^{\mathit{C}})\kappa (E_i,E_j|{\lambda }_E^{\mathit{C}})\text{.}$$(24)

The hyperparameters ${\mathrm{\lambda }}_{\mathit{x}}^{\mathit{C}}$ are for the coordinate axis indicated by x ∈ {En, A, …} and can take different values for κ_L and κ_H (C ∈ {L, H}).

The transition between the two energy domains, which means to switch from κ_L to κ_H, is established by the logistic function $$\sigma \mathrm{}(x\mathrm{}|k,\mathrm{}{x}_0)\mathrm{}=\frac{1}{1+\exp (\mathrm{}-k\mathrm{}(x-x_0)\mathrm{})}\mathrm{}\text{.}$$(25)

Noteworthy, all the variables are vectors. The equation $k(\mathrm{}x-x_0)=0$ defines a (hyper)plane with normal vector k and distance $\mathrm{}|k{x}_0|\mathrm{}/|k\mathrm{}|$ to the origin of the coordinate system. The function in equation (25) attains values close to zero for x far away from the plane on one side and close to one if far away on the other side. Within which distance to the plane the transition from zero to one occurs depends on the length of k. The larger |k|, the faster the transition and the narrower the window of transition around the plane.

With the abbreviations $${\tau }_L(x_i,x_j|k,x_0)=\sigma (x_i|k,x_0)\sigma (x_j|k,x_0)\text{,}$$(26) $${\tau }_H(\mathrm{}{x}_i\mathrm{}{x}_j\mathrm{}|k\mathrm{},x_0)=(1-\sigma \mathrm{}(x_i\mathrm{}|k\mathrm{}{x}_0))\mathrm{}$$(27) $$\times (1-\sigma (x_j|k,x_0))\text{,}$$(28) the full covariance function κ_full is given by $${\kappa }_{\text{full}}\mathrm{}(x_i\mathrm{},x_j)\mathrm{}={\tau }_L(\mathrm{}{x}_i,\mathrm{}{x}_j|\mathrm{}k,\mathrm{}{x}_0)\mathrm{}{\kappa }_L(\mathrm{}{x}_i\mathrm{},x_j|\mathrm{}{\delta }_L\mathrm{},{\lambda }_H)\mathrm{}+{\tau }_H(x_i,x_j|k,x_0){\kappa }_H(x_i,x_j|{\delta }_H,{\lambda }_H)\text{.}$$(29)

Finally, the GP regression is not performed on the absolute difference Δ between a model prediction σ_mod and experimental data point σ_exp, but on the transformed quantity $${\displaystyle \tilde{\mathrm{\Delta }}}=\frac{({\sigma }_{\exp }-{\sigma }_{\mathrm{mod}})}{\max ({\sigma }_{\mathrm{mod}},0.1)}\text{.}$$(30)

In words, relative differences are taken for model predictions larger than 0.1 and absolute differences scaled up by a factor of ten for model predictions below 0.1. Relative differences fluctuate usually wildly for spectrum values close to zero–especially for a Monte Carlo code such as INCL–and the switch to absolute values helps GP regression to find more meaningful solutions with a better ability to extrapolate.

Due to the number of roughly ten thousand data points, the covariance matrices computed with equation (29) were replaced by the sparse approximation outlined in Section 2.2. I introduced three hundred pseudo-inputs and placed them randomly at the locations of the experimental data. Their locations were then jointly optimized with the hyperparameters, which will be discussed in Section 3.3.

The diagonal matrix D₂ occurring in the approximation was changed to $${{\displaystyle \tilde{D}}}_2=D_2+B+P\text{,}$$(31) to accommodate statistical uncertainties of the model prediction (due to INCL being a Monte Carlo code) and the experimental data. Both B and P are diagonal matrices. The matrix B contains variances corresponding to 10% statistical uncertainty for all experimental data points. The matrix P contains the estimated variances of the model predictions.

A diagonal matrix for the experimental covariance matrix B can certainly be challenged because the important systematic errors of experiments reflected in off-diagonal elements are neglected. This is at the moment a limitation of the approach.

3.3 Marginal likelihood maximization

The hyperparameters appearing in equation (29) and the locations of the three hundred pseudo-inputs $\left\{x_k^{psi}\right\}$ determining the approximation in equation (11) were adjusted via marginal likelihood maximization described in Section 2.3. To be explicit, the hyperparameters considered were δ_L, δ_H, $${\mathrm{\lambda }}_L=({\lambda }_{\text{En}}^L,{\lambda }_A^L,{\lambda }_Z^L,{\lambda }_{\theta }^L,{\delta }_E^L)\text{,}$$(32) $${\mathrm{\lambda }}_H\mathrm{}=({\lambda }_{\text{En}}^H\mathrm{},{\lambda }_A^H\mathrm{},{\lambda }_Z^H\mathrm{},{\lambda }_{\theta }^H,\mathrm{}{\delta }_E^H),\mathrm{}$$(33) and also x₀ and k of the logistic function. The vector k was forced to remain parallel to the axes associated with En, A, and Z, i.e. $\mathit{k}={(\mathrm{0,0,0},k_{\theta },k_E)}^T$. Further, polar coordinates k_θ = k_csinγ, k_E = k_ccosγ were introduced. The direction of the vector x₀ was taken equal to that of k, which removes ambiguity in the plane specification without shrinking the set of possible solutions. Because of these measures, it sufficed to consider the length x_c = |x₀| as hyperparameter. Counting both hyperparameters and pseudo-inputs, 1515 parameters were taken into account in the optimization.

I employed the L-BFGS-B algorithm [15] as implemented in the optim function of R [19], which makes use of an analytic gradient, can deal with a large number of variables and permits the specification of range restrictions. The optimization was performed on a cluster using 25 cores and was stopped after 3500 iterations, which took about 10 h. The obtained solution corresponds to χ²/N = 1.03 and is with a two-sided p-value of 0.04 reasonably consistent in a statistical sense. Restrictions of parameter ranges were established to introduce prior knowledge and to guide the optimization procedure. Noteworthy, lower limits on length-scales, such as ${\lambda }_{\theta }^L$ and ${\lambda }_E^L$; were introduced to counteract their dramatic reduction due to inconsistent experimental data in the same energy/angle range. Table 2 summarizes the optimization procedure. The evolution of the pseudo-inputs projected onto the (A,Z)-plane is visualized in Figure 1.

Fig. 1 Dispersion of the pseudo-inputs during the optimization procedure shown in the (A,Z)-projection. The black points are the initial positions and the green points the final positions of the pseudo-inputs.

A thorough study of the optimization process exceeds the scope of this paper and hence I content myself with a few remarks. The length scales associated with the emission angle, i.e. ${\lambda }_{\theta }^{\mathit{C}}$, experienced significant changes. Their increase means that the model bias is similar for emission angles far away from each other and the GP process is able to capture them. Concerning the length scales associated with the emission energy, the small value ${\lambda }_E^L=4.8$ compared to the larger value ${\lambda }_E^H=43$ indicates that the features of the model bias of the low and high energy domain are indeed different. The most striking feature, however, is that the large length scales ${\lambda }_A^C$ and ${\lambda }_Z^C$ “survived” the optimization, which means that the model bias behaves similar over large regions of the nuclide charts. Examples of isotope extrapolations will be given in Section 3.4.

As a final remark, a more rigorous study of the optimization procedure is certainly necessary and there is room for improvement. This is left as future work.

3.4 Results and discussion

The values of the hyperparameters and pseudo-inputs obtained by marginal likelihood maximization were used in the covariance function in equation (29). This covariance function was employed to compute the required covariance matrices in equations (15) and (16) based on the experimental data summarized in Table 1. Because the hyperparameters are determined during hyperparameter optimization before being used in the GP regression, they will be referred to as prior knowledge.

Equations (15) and (16) enable the prediction of a plethora of spectrum values and their uncertainties for combinations of incident energy, mass number, charge number, emission angle, and emission energy. The few selected examples of predictions in Figure 2 serve as the basis to discuss general features of GP regression, the underlying assumptions, and its accuracy and validity.

How well we can interpolate between data and how far we can extrapolate beyond the data depends on the suitability of the covariance function for the problem at hand. The building block for the full covariance function in equation (29) is the one-dimensional covariance function in equation (21). Using the latter imposes the assumption that possible solutions have derivatives of any order and hence are very smooth [1, Chap. 4]. Interpolations between data points included in the regression are determined by this smoothness property and values of the length scales ${\mathit{\lambda }}_{\mathit{x}}^{\mathit{C}}$.

The length scales reflect the prior assumption about similarity between spectrum values of points a certain distance away from each other. This prior assumption directly impacts the uncertainty of the predictions. The farther away a prediction is from an observation point, the higher the associated uncertainty. If a prediction is already multiples of any length scale away from all observations, the uncertainty reverts to its prior value given by either δ_L or δ_H depending on the energy domain.

In the case of the sparse approximation, the uncertainty is related to the distance to the pseudo-inputs. Because only few pseudo-inputs are located at very high and very low emission energies, the 2σ uncertainty bands in Figure 2 in those energy domains are rather large despite the presence of experimental data.

The important finding in this specific application is that the length scales related to emission angle, ${\mathit{\lambda }}_{\mathit{E}}^{\mathit{C}}$, mass number, ${\mathit{\lambda }}_{\mathit{A}}^{\mathit{C}}$, and charge number, ${\mathit{\lambda }}_{\mathit{E}}^{\mathit{C}}$, are very large. GP regression is therefore able to interpolate and extrapolate over large ranges of these axes.

The interpolation of spectrum values between angles and emission energies of an isotope with available data may be considered rather standard. For instance, one can do a χ²-fit of a low order Legendre polynomial to the angular distributions for emission energies with data. The coefficients of the Legendre polynomial for intermediate emission energies without data can then be obtained by linear interpolation.

The important difference between such a conventional fit and GP regression is the number of basis functions. Whereas their number is fixed and limited in a conventional fit, GP regression amounts to a fit with an infinite number of basis functions [1, Sect. 2.2]. The length scales regulate the number of basis functions that effectively contribute to the solution. Hyperparameter optimization decreases the length scales for unpredictable data which leads to a greater number of contributing basis functions and consequently to greater flexibility and larger uncertainties in the predictions. This feature sets GP regression apart from a standard χ²-fit. In the latter, the uncertainty of the solution depends to a much lesser extent on the (un)predictability of the data and much more on the number of data points and the fixed number of basis functions.

One truly novel element in the approach is the inclusion of the mass and charge number, which enables predictions for isotopes without data. We can easily imagine that different isotopes differ significantly by their physical properties. From this perspective, the idea to extrapolate the model bias to other isotopes should be met with healthy skepticism.

To get a first grasp on the validity of isotope extrapolations, let us consider again the hyperparameter optimization discussed in Section 3.3. The hyperparameters were adjusted on the basis of the isotopes in Table 1. These data are spread out over the periodic table and cover a range from carbon to thorium. In these data, similar trends of the model bias persist across significant ranges of the mass and charge number, which was the reason that the associated length scales retained high values during optimization. For instance, the experimental data of carbon and indium in Figure 2 show comparable structures of the model bias despite their mass differences.

However, the isotopes considered in the optimization are not very exotic and gaps between them are at times large. Further, these isotopes are certainly not a random sample taken from the periodic table and therefore most theoretical guarantees coming from estimation theory do not hold. So how confident can we be about isotope extrapolation?

To provide a basis for the discussion of this question, Figure 2 also contains predictions for oxygen and cadmium. Importantly, the associated experimental data have not been used for the hyperparameter optimization and in the GP regression. The predictions follow well the trend of the experimental data. The supposedly 2σ uncertainty bands, however, include less than 95% of the data points. This observation points out that there are systematic differences between isotopes. Therefore, due to the sample of isotopes not being a random sample, uncertainty bands should be interpreted with caution.

One way to alleviate this issue could be to add a covariance function to equation (29) which only correlates spectrum values of the same isotope but not between isotopes. This measure would lead to an additional uncertainty component for each isotope, which only decreases if associated data are included in the GP regression.

A very extreme mismatch up to 500% between the prediction and experimental data occurs for oxygen at 60° and emission energies below 1 MeV. The creation of low energy neutrons is governed by de-excitation processes of the nucleus suggesting that the angular distribution of emitted particles is isotropic. This property holds for the model predictions but not for the experimental data. The experimental data exhibits a peak at 60° which is about a factor five higher than at 30°, 120° and 150°. The origin of this peculiarity may deserve investigation but is outside the scope of this work. For the sake of argument I assume that it is indeed a reflection of some property of the nucleus.

Because the data in Table 1 do not include any measurements below 1 MeV emission energy in this mass range, both hyperparameter optimization and GP regression had no chance to be informed about such a feature. In this case, the predictions and uncertainties are determined by nearby values or – if there are not any – by the prior expectation.

If we would have been aware of such effects, we could have used it as a component in the covariance function to reflect our large uncertainty about the spectrum for low emission energies. Otherwise, the only sensible data-driven way to provide predictions and uncertainties for unobserved domains is to assume that effects are similar to those in some observed domains. Of course, it is our decision which domains are considered similar. The employed covariance function in equation (29) incorporates the assumption that nearby domains in terms of mass charge, angle, etc. are similar. However, we are free to use any other assumption about similarity to construct the covariance function. As a side note, also results of other models relying on different assumptions could serve as a basis for uncertainties in unobserved regions. Such information can be included in GP regression in a principled way.

Fig. 2 Model bias of INCL in terms of equation (30) in the (p,X)n double differential spectra for 800 MeV incident protons and different isotopes as predicted by GP regression. A missing mass number behind the isotope symbol indicates natural composition. The uncertainty band of the prediction and the error bars of the experimental data denote the 2σ confidence interval. Carbon and indium were taken into account in the GP regression but not cadmium and oxygen. The experimental data are colored according to the associated access number (ACCNUM) in the EXFOR database. This shows that all displayed data come from just three experiments [20–22].

4 Summary and outlook

Sparse GP regression, a non-parametric estimation technique of statistics, was employed to estimate the model bias of the C++ version of the Liège INCL coupled to the evaporation code ABLA07. Specifically, the model bias in the prediction of double-differential inclusive neutron spectra over the complete nuclide chart was investigated for incident protons above 100 MeV. Experimental data from the EXFOR database served as the basis of this assessment. Roughly ten thousand data points were taken into account. The hyperparameter optimization was done on a computing cluster whereas the GP regression itself on a desktop computer. The obtained timings indicate that increasing the number of data points by a factor of ten could be feasible.

For this specific application, it was shown that GP regression produces reasonable results for isotopes that have been included in the procedure. It was argued that the validity of predictions and uncertainties for isotopes not used in the procedure depends on the validity of the assumptions made about similarity between isotopes. As a simple benchmark, the isotopes oxygen and cadmium, which have not been taken into account in the procedure, were compared to the respective predictions. The agreement between prediction and experimental data was reasonable but the 95% confidence band sometimes misleading and should therefore be interpreted with caution. Accepting the low energy peak of oxgyen at 60° in the data as physical reality, the low energy spectrum of oxygen served as an example where the similarity assumption between isotopes did not hold.

As for any other uncertainty quantification method, it is a hard if not impossible task to properly take into account unobserved phenomena without any systematic relationship to observed ones. The existence of such phenomena are unknown unknowns, their observation is tagged shock, surprise or discovery, and they potentially have a huge impact where they appear.

Even though GP regression cannot solve the philosophical problem associated with unknown unknowns, knowledge coming from the observation of new effects can be taken into account by modifying the covariance function. For instance, the peculiar peak in the oxygen spectrum suggests the introduction of a term in the covariance function which increases the uncertainty of the spectrum values in the low emission energy domain. The ability to counter new observations with clearly interpretable modifications of the covariance function represents a principled and transparent way of knowledge acquisition.

The scalability and the possibility to incorporate prior assumptions by modeling the covariance function makes sparse GP regression a promising candidate for global assessments of nuclear models and to quantify their uncertainties. Because the formulas of GP regression are structurally identical to those of many existing evaluation methods, GP regression should be regarded more as a complement than a substitute, which can be interfaced with existing evaluation methods.

Acknowledgments

This work was performed within the work package WP11 of the CHANDA project (605203) financed by the European Commission. Thanks to Sylvie Leray, Jean-Christophe David and Davide Mancusi for useful discussion.

References

All Tables

All Figures

	Fig. 1 Dispersion of the pseudo-inputs during the optimization procedure shown in the (A,Z)-projection. The black points are the initial positions and the green points the final positions of the pseudo-inputs.
In the text

Fig. 2 Model bias of INCL in terms of equation (30) in the (p,X)n double differential spectra for 800 MeV incident protons and different isotopes as predicted by GP regression. A missing mass number behind the isotope symbol indicates natural composition. The uncertainty band of the prediction and the error bars of the experimental data denote the 2σ confidence interval. Carbon and indium were taken into account in the GP regression but not cadmium and oxygen. The experimental data are colored according to the associated access number (ACCNUM) in the EXFOR database. This shows that all displayed data come from just three experiments [20–22].

In the text