# Learning Perturbations to Explain Time Series Predictions

Joseph Enguehard

###### Abstract

Explaining predictions based on multivariate time series data carries the
additional difficulty of handling not only multiple features, but also time
dependencies. It matters not only what happened, but also when, and the same
feature could have a very different impact on a prediction depending on this
time information. Previous work has used perturbation-based saliency methods
to tackle this issue, perturbing an input using a trainable mask to discover
which features at which times are driving the predictions. However these
methods introduce fixed perturbations, inspired from similar methods on static
data, while there seems to be little motivation to do so on temporal data. In
this work, we aim to explain predictions by learning not only masks, but also
associated perturbations. We empirically show that learning these
perturbations significantly improves the quality of these explanations on time
series data.

Machine Learning, ICML

##  1 Introduction

Figure 1: Illustration of different perturbations on a time-series (first
plot). We aim to evaluate the importance of the third spike here. It is very
likely that this spike is unimportant, as it is a regularity in the data, and
only the last, larger spike, could matter. However, using a Gaussian blur
(second plot) or replacing the original data with an average (third plot)
changes the input significantly, which could lead an explanation method to
wrongly state that this spike is important. Our learned perturbation (last
plot) should on the other hand replace the explained data with another spike,
leading to little difference in the output, and therefore correctly stating
that this data is unimportant.

Explaining neural networks predictions has received increasing attention, as
these models become more embedded in many decision processes. It is indeed
important to understand why such models, which are intrinsically difficult to
explain and are often qualified as “black boxes”, made a specific prediction.
This information is crucial to assess the fairness of an algorithm in
impactful situations, such as when providing a medical diagnosis (Ahuja, 2019)
or computing a credit score (Moscato et al., 2021). Explaining a deep neural
network’s predictions is also important to build trust for the users of such
technologies. In the medical field, this has been proven to be a crucial step
in building this trust (LaRosa & Danks, 2018).

As a result, multiple methods to explain why a model made a specific
prediction have been recently developed. Some of these methods, such as Lime
(Ribeiro et al., 2016) or Shap (Lundberg & Lee, 2017) explain a model’s
predictions by approximating it locally using a transparent method, in this
case a weighted linear regression. Other methods aim to specifically explain a
neural network’s output by leveraging the back-propagation algorithm to
compute the gradient of an output w.r.t. an input (Simonyan et al., 2013). A
feature with an associated high gradient can indeed be interpreted as
important, the sign of this gradient indicating if this feature influences the
prediction positively or negatively. Several variations of this method have
also been developed (Sundararajan et al., 2017; Shrikumar et al., 2017).

Another important class of explanation methods is called perturbation-based.
These methods consist in perturbing a feature or a group of features, and
measuring how the resulting prediction changes. A greater change indicates a
higher importance of the perturbed features. Such methods include Occlusion
(Zeiler & Fergus, 2014), which masks features to estimate their importance.
Extremal Masks (Fong & Vedaldi, 2017; Fong et al., 2019) is another
perturbation-based method, which learns a mask used to perturb the input. We
present this method in more detail in the next section.

However, while many explanation methods have been proposed to explain a neural
network, few have been developed to handle multivariate time series data. Yet,
this type of data is especially important in the medical field, where the data
can be a list of timestamped medical events, or of vitals measurements. There
is therefore a need to adapt explanation methods to handle this temporal
element. These adaptations currently include RETAIN (Choi et al., 2016), an
attention-based model which learns this attention over features and time, or
FIT (Tonekaboni et al., 2020), which estimates the importance of features over
time by quantifying the shift in the predictive distribution. Another method,
DynaMask (Crabbé & Van Der Schaar, 2021), adapts perturbation-based methods to
multivariate time-series. We will present and discuss this method further in
the next section.

In this work, we aim to further adapt perturbation-based methods to
multivariate time-series driven with the following insight. In the works of
Fong & Vedaldi (2017) and Crabbé & Van Der Schaar (2021), while the mask is
learned, the perturbation induced by this mask is fixed. For instance, Fong &
Vedaldi (2017) replaces a feature with a Gaussian blur (a weighted average of
data around the feature) depending on the value of the feature’s mask: the
lower this value, the higher the amount of blur. Crabbé & Van Der Schaar
(2021) adapts this method by blurring the data temporally. This method seems
reasonable with images, where information can be assumed to be local, which
explains why convolutional neural networks (CNNs), which have a limited filter
size, still perform very well on such data. However, multivariate time-series
can have long-term dependencies which makes it less obvious to use a temporal
Gaussian blur as the perturbation. Instead of replacing a masked feature with
a local average, we might want to replace it using data further away in time.
But then, how should we choose the correct perturbation formula?

This calls to replace fixed perturbations with learnable ones. In this work,
we present such a method111An implement of this work can be found at
https://github.com/josephenguehard/time_interpret and empirically show that it
significantly improves the quality of the explanations, evaluated on both
synthetic and real-world data. This study is organised as follows. We first
present in more detail the methods of Fong & Vedaldi (2017) and Crabbé & Van
Der Schaar (2021) in the next section. We then present our method in the
following one. We conduct several experiments in the next section, designed to
compare our method with several baselines, and we provide elements of
discussion in the last section.

##  2 Background Work

In this section, we describe in more detail two methods: one developed by Fong
& Vedaldi (2017) and its adaptation to time series by Crabbé & Van Der Schaar
(2021).

Fong & Vedaldi (2017) propose a perturbation-based method which is defined as
following. A trainable mask, with values restricted between 0 and 1, is used
to generate perturbed data, which is then passed to the neural network to be
explained in order to compute predictions. This mask can then be trained in
two different manners, that the authors call the deletion game and the
preservation game. In the deletion game, we aim to mask as little data as
possible, while trying to reduce as much as possible the predictions, on the
targeted class, of the perturbed data, compared with the original predictions.
This objective can be defined as, for a model
$f:\mathbb{R}^{n}\to\mathbb{R}^{p}$, a mask $\textbf{m}\in[0,1]^{n}$, an input
$\textbf{x}\in\mathbb{R}^{n}$ and a perturbation
$\Phi(\textbf{x},\textbf{m}):\mathbb{R}^{n}\times[0,1]^{n}\to\mathbb{R}^{n}$:

$\operatorname*{arg\,min}_{\textbf{m}\in[0,1]^{n}}\lambda||\textbf{1}-\textbf{m}||_{1}-\mathcal{L}(f(\textbf{x}),f(\Phi(\textbf{x},\textbf{m})))$
(1)

The value n represents the input dimension, and $\lambda$ is a hyperparameter
balancing both goals.

Secondly, in the preservation game, we aim to retain the least amount of data
that will preserve the closest predictions compared with the original ones on
the targeted class. This objective can be defined as:

$\operatorname*{arg\,min}_{\textbf{m}\in[0,1]^{n}}\lambda||\textbf{m}||_{1}+\mathcal{L}(f(\textbf{x}),f(\Phi(\textbf{x},\textbf{m})))$
(2)

Moreover, the perturbation $\Phi(\textbf{x},\textbf{m})$ is fixed given an
input and a mask. Fong & Vedaldi (2017) define several strategies
222Unintuitively, the original data is masked when m = 0. We kept this
notation as it is used in both Fong & Vedaldi (2017) and Crabbé & Van Der
Schaar (2021).:

$\Phi(\textbf{x},\textbf{m})=\begin{cases}\textbf{m}\times\textbf{x}+(\textbf{1}-\textbf{m})\times\mu_{0}\\\
\textbf{m}\times\textbf{x}+(\textbf{1}-\textbf{m})\times\nu\\\ \int
g_{\sigma_{0}\times(1-\textbf{m})}(\textbf{y}-\textbf{x})\,d\textbf{y}\end{cases}$
(3)

The first strategy corresponds to replacing the original masked value x with
an average $\mu_{0}$, the second one consisting in replacing this value with
Gaussian noise: $\nu\sim\mathcal{N}(0,1)$ and the last replaces it with a
Gaussian blur $g_{\sigma}$ around x, given a maximum std $\sigma_{0}$. Fong &
Vedaldi (2017) also add some regularisation to ensure the perturbation to be
more natural in the context of computer vision, but we leave it out, as it is
further away from our topic.

While this method was developed to explain predictions based on images, Crabbé
& Van Der Schaar (2021) adapted it to multivariate time series. They propose
as a result a method they call DynaMask, as the learned mask contains in this
case a time dimension. The input space is now $\mathbb{R}^{\mathrm{T}\times
n}$, and we consider similarly a neural network $f$ and a target class $c$
such as: $f_{c}(\textbf{x}):\mathbb{R}^{\mathrm{T}\times n}\to\mathbb{R}$.
Therefore, the mask $\textbf{m}\in\mathbb{R}^{\mathrm{T}\times n}$ and the
input $\textbf{x}\in\mathbb{R}^{\mathrm{T}\times n}$ are also defined on this
input space.

The main contribution of Crabbé & Van Der Schaar (2021) is to then adapt the
perturbation operator $\Phi$ to account for this temporal information. They
also introduce three strategies:

$\Phi(\textbf{x},\textbf{m})_{t,i}=\begin{cases}m_{t,i}\times
x_{t,i}+(1-m_{t,i})\times\mu_{t,i}\\\ m_{t,i}\times
x_{t,i}+(1-m_{t,i})\times\mu^{p}_{t,i}\\\
\frac{\sum_{t^{\prime}=1}^{\textrm{T}}x_{t^{\prime},i}\times
g_{\sigma(m_{t,i})}(t-t^{\prime})}{\sum_{t^{\prime}=1}^{\textrm{T}}g_{\sigma(m_{t,i})}(t-t^{\prime})}\end{cases}$
(4)

Where $\mu_{t,i}$ is an average of $\textbf{x}_{:,i}$ over a window W around
$t$:

$\mu_{t,i}=\frac{1}{2\textrm{W}+1}\sum_{t^{\prime}=t-\textrm{W}}^{t+\textrm{W}}x_{t^{\prime},i}$
(5)

and $\mu^{p}_{t,i}$ is an average of $\textbf{x}_{:,i}$ over a past element up
to $t$:

$\mu^{p}_{t,i}=\frac{1}{\textrm{W}+1}\sum_{t^{\prime}=t-\textrm{W}}^{t}x_{t^{\prime},i}$
(6)

Finally, the last perturbation is a temporal Gaussian blur:

$g_{\sigma(m_{t,i})}(t)=\exp(-\frac{t^{2}}{2\sigma^{2}});\,\sigma(\textbf{m})=\sigma_{\textrm{max}}(\textbf{1}-\textbf{m})$
(7)

Crabbé & Van Der Schaar (2021) uses these perturbations in a preservation
game, which aims to mask the maximum amount of data while keeping close
predictions compared with the originals. They also leverage further work of
Fong & Vedaldi (2017): Fong et al. (2019), which replaces Equations 1 and 2
with an area constraint. In the preservation mode (the deletion mode can be
adapted similarly), the regulation $\lambda||\textbf{m}||_{1}$ in Equation 2
is replaced with:
$\lambda_{a}(\textbf{m})=||\textrm{vecsort}(\textbf{m})-\textbf{r}_{a}||^{2}$,
where $a$ is a number between 0 and 1, $\textrm{vecsort}(\textbf{m})$ sorts
the values of m from lowest to largest, and $\textbf{r}_{a}$ is a vector
containing $(1-a)\times T\times n$ zeros followed by $a\times T\times n$. As a
result, this constraint allows the user to define how much of the data should
be masked. In practice, Crabbé & Van Der Schaar (2021) use $a$ as a
hyperparameter, which is tuned for each data point to be explained.

##  3 Method

Figure 2: Illustration of our method. The input is passed through a neural
network NN to create a perturbation. A mask m is then used to balance the
amount of perturbed data: NN(x) and unperturbed data: x, resulting in
$\Phi(\textbf{x},\textbf{m})$. Both x and $\Phi(\textbf{x},\textbf{m})$ are
then passed through the model to explain f. Learnable parameters (m and NN(x))
are presented in continuous boxes, while fixed parameters (the model f) are
presented in dashed boxes. The objective of this method is to keep the
predictions of the perturbed data as close as possible to the original ones,
while masking as much data as possible and to keep the perturbations NN(x) as
sparse as possible. The overall goal is therefore to identify which features
are salient enough to be sufficient to recover the original predictions, when
all other features are masked.

While Crabbé & Van Der Schaar (2021) propose temporal perturbations as
adaptations of the ones defined by Fong & Vedaldi (2017) in a computer vision
context, these perturbations are kept fixed and local. They are indeed defined
either as a moving average perturbation, or as a temporal Gaussian blur.
However, temporal data is often characterised by long-term dependencies, and
local information can therefore be insufficient to determine the importance of
a feature at a particular time. For instance, temporal data can include
repetitive patterns, as illustrated on Figure 1, which cannot be taken account
using only temporally local information. Moreover, while the perturbations
proposed by Crabbé & Van Der Schaar (2021) do include the possibility to
include data further away in time, by tuning the size of the window W, or the
parameter $\sigma_{\max}$ for the Gaussian blur, it is not clear how to choose
such parameters nor how this would solve the issue of long term patterns.

This insight calls for a generalised perturbation, which can be tuned to the
data we are aiming to explain. A first idea would be to directly learn this
perturbation $\Phi(\textbf{x})$, without needing a mask, by optimizing a
function similar to Equation 2. However, this method is problematic as it
gives too much liberty to the perturbation model. Indeed, such a model,
incentivised to output sparse explanations, could compress the data
information into a small part of the input space, stating that this part is
important while the rest is uninformative. On the contrary, we need to
constrain the perturbation operator to explain each part of the input data
without changing or moving it.

To overcome this difficulty, we take inspiration from the perturbation
operators of Crabbé & Van Der Schaar (2021) in Equation 4. These perturbations
are generally defined as
$\textbf{m}\times\textbf{x}+(1-\textbf{m})\times\mu(\textbf{x})$, where
$\mu(\textbf{x})$ is a function of the input. In this work, we propose to
replace these fixed functions with a neural network (NN), and to train it in
combination with the mask. Our perturbation is therefore defined as:

$\begin{split}\Phi(\textbf{x},\textbf{m})&=\textbf{m}\times\textbf{x}+(\textbf{1}-\textbf{m})\times\textrm{NN}(\textbf{x})\\\
\textbf{0}&\leq\textbf{m}\leq\textbf{1}\end{split}$ (8)

By keeping m between 0 and 1, we constrain the mask to only learn how
important each feature is. Moreover, we can see that this equation can be
interpreted as a generalisation of the perturbations from Crabbé & Van Der
Schaar (2021) defined in Equation 4. The neural network in the second
component of Equation 8 can indeed, after training, output a Gaussian blur or
an average of x over a window.

In practice, we want to model $\textrm{NN}(\textbf{x})$ as a weighted sum of
$x_{t,i},\,t\in\\{1,...,\textrm{T}\\}$. As a result, we choose this model to
be a bidirectional GRU (Cho et al., 2014). This would correspond to a general
form of a Gaussian blur or a window around each element $x_{t,i}$. We also
compare this choice, in the experiment section, with a unidirectional GRU,
which would be closer to the $\mu^{p}_{t,i}$ average in Crabbé & Van Der
Schaar (2021).

As in Crabbé & Van Der Schaar (2021), we define the objective of the mask and
the GRU combined as a preservation game, aiming to mask as much data as
possible while keeping the closest predictions as possible to the original
ones. Our objective is therefore:

$\operatorname*{arg\,min}_{\textbf{m},\Theta\in\textrm{NN}}\lambda||\textbf{m}||_{1}+\mathcal{L}(f(\textbf{x}),f(\Phi(\textbf{x},\textbf{m})))$
(9)

where $\Theta$ represents the parameters of the neural network, and
$\mathcal{L}$ represents a loss between the original and the perturbed
predictions. This loss can be for instance a mean square error for regression
tasks, or a cross-entropy loss for classification tasks.

One issue that can arise from this objective is that the neural network can be
rewarded to mimic the original x data. Indeed, we can see from Equation 8
that, if $\textbf{m}=\textbf{0}$, then
$\Phi(\textbf{x},\textbf{m})=\textrm{NN}(\textbf{x})$. Moreover, if
$\textrm{NN}(\textbf{x})\approx\textbf{x}$, the objective defined in Equation
9 is approximately zero. To prevent this behavior, we modify Equation 9 with
the following one:

$\operatorname*{arg\,min}_{\textbf{m},\Theta\in\textrm{NN}}\lambda_{1}||\textbf{m}||_{1}+\lambda_{2}||\textrm{NN}(\textbf{x})||_{1}+\mathcal{L}(f(\textbf{x}),f(\Phi(\textbf{x},\textbf{m})))$
(10)

In Equation 10, we therefore force the perturbations to be minimal, being not
null only when there is an incentive to do so. Indeed, in Equation 2, there is
a balance on $\Phi$: $||\textbf{m}||_{1}$ tends to make $\Phi$ uninformative,
while $\mathcal{L}$ does the opposite. Equation 10 differs in that sense from
Equation 2, as $||\textbf{m}||_{1}$ tends to make $\Phi$ close to NN(x), which
is not necessarily uninformative. To entice NN(x) to be uninformative, we add
the loss $||NN(\textbf{x})||_{1}$, using zero as a prior. Therefore, breaking
down the objective of Equation 10, we have:

  * •

$||\textbf{m}||_{1}$ induces $\Phi(\textbf{x})$ to be close to
$\textrm{NN}(\textbf{x})$

  * •

$||\textrm{NN}(\textbf{x})||_{1}$ induces $\Phi(\textbf{x})$ to be close to 0
(uninformative)

  * •

$\mathcal{L}$ induces $f(\Phi(\textbf{x},\textbf{m}))$ to be close to
$f(\textbf{x})$ (informative)

We also set $\lambda_{1}=\lambda_{2}=1$ in our experiments, while an ablation
study on the choice of these hyperparameters can be found in Section 4 and
Appendix A.

Moreover, contrary to Crabbé & Van Der Schaar (2021), we do not use an area
constraint $||\textrm{vecsort}(\textbf{m})-\textbf{r}_{a}||^{2}$, as it is not
clear how to choose the hyperparameter $a$ on usually complex data. In
practice, Crabbé & Van Der Schaar (2021) tune this hyperparameter, which is
computationally expensive, as it requires to train multiple masks. We propose
instead to directly train our model using Equation 10.

##  4 Experiments

Following Tonekaboni et al. (2020) and Crabbé & Van Der Schaar (2021), we
perform experiments on two datasets: a synthetic one, generated using a Hidden
Markov model, and a real-world one, MIMIC-III (Johnson et al., 2016).

###  4.1 Hidden Markov model experiment

We generate data using a 2-state hidden Markov model (HMM), closely following
Crabbé & Van Der Schaar (2021). The state $s_{t}$ can therefore be either 0 or
1, and we generate 200 states: $t\in[1:200]$. Moreover, the input vector has
three features, generated according to the current state:
$\textbf{x}_{t}\sim\mathcal{N}(\boldsymbol{\mu}_{s_{t}},\boldsymbol{\Sigma}_{s_{t}})$.
The label $y_{t}$ is generated only using the last two features, the first one
being irrelevant. The choice of which feature is used to generate the label
depends on the state:

$\begin{split}y_{t}&\sim(1+\exp(x_{2,t})^{-1})\mspace{10.0mu}\textrm{if}\mspace{10.0mu}s_{t}=0\\\
y_{t}&\sim(1+\exp(x_{3,t})^{-1})\mspace{10.0mu}\textrm{if}\mspace{10.0mu}s_{t}=1\end{split}$
(11)

Please refer to Crabbé & Van Der Schaar (2021) for more details on this
dataset, in particular in the choice of $\boldsymbol{\mu}_{s_{t}}$ and
$\boldsymbol{\Sigma}_{s_{t}}$.

We generate 1000 time series using this method, and train a one-layer GRU (Cho
et al., 2014) neural network to predict $y_{t}$ using $\textbf{x}_{t}$, which
we aim to explain.

As we know the true salient features with this dataset, we evaluate our
explanation methods by comparing the similarity between salient features
produced by each method and the ground truth. To do so, we use standard
classification metrics: area under recall (AUR) and area under precision
(AUP). We also use two metrics introduced by Crabbé & Van Der Schaar (2021):
Information:
$I_{\textbf{m}}(\textbf{a})=-\sum_{(t,i)\in\textbf{a}}\ln(1-m_{t,i})$ which is
analogous to the Shannon information content. A higher value indicates a more
informative mask. The second metric is the mask entropy:
$S_{\textbf{m}}(\textbf{a})=-\sum_{(t,y)\in\textbf{a}}m_{t,i}\ln
m_{t,i}+(1-m_{t,i})\ln(1-m_{(t,i)})$ which is analogous to Shannon entropy. In
both metrics, a corresponds to the true salient features.

We compare our method with the following ones: DeepLift (Shrikumar et al.,
2017), DynaMask (Crabbé & Van Der Schaar, 2021), Integrated Gradients (IG)
(Sundararajan et al., 2017), GradientShap (Lundberg & Lee, 2017), Fit (Tjoa &
Guan, 2020), Lime (Ribeiro et al., 2016), Augmented Occlusion (Tonekaboni et
al., 2020), Occlusion (Zeiler & Fergus, 2014) and Retain (Choi et al., 2016).
Furthermore, our method uses a bidirectional GRU for the perturbation model.

Method | AUP $\uparrow$ | AUR $\uparrow$ | I $\uparrow$ | S $\downarrow$
---|---|---|---|---
DeepLift | 0.920 (0.019) | 0.454 (0.011) | 359 (9.55) | 145 (0.949)
DynaMask | 0.711 (0.020) | 0.763 (0.026) | 954 (50.0) | 45.4 (0.781)
IG | 0.918 (0.019) | 0.454 (0.011) | 359 (11.6) | 146 (0.871)
GradientShap | 0.849 (0.030) | 0.414 (0.015) | 335 (14.8) | 138 (2.44)
Fit | 0.421 (0.013) | 0.549 (0.017) | 436 (22.7) | 164 (2.79)
Lime |  0.932 (0.017) | 0.438 (0.008) | 347 (8.46) | 143 (1.47)
Occlusion | 0.866 (0.032) | 0.393 (0.006) | 322 (14.6) | 137 (1.90)
Aug Occlusion | 0.755 (0.043) | 0.388 (0.025) | 364 (9.02) | 165 (1.42)
Retain | 0.645 (0.088) | 0.334 (0.013) | 206 (21.2) | 138 (5.85)
Ours | 0.885 (0.030) |  0.781 (0.013) |  1536 (79.0) |  34.1 (3,70)

Table 1: Results of each explanation method compared with ours. For each
metric, $\uparrow$ indicates that higher is better, and $\downarrow$ that
lower is better. Mean and std are reported over 5 folds.

We present our results in Table 1. These results333In Tables 1 and 4, some
results differ from Crabbé & Van Der Schaar (2021) due to a few issues in
their original implementation. Please refer to issues 4, 8 and 9 in
https://github.com/JonathanCrabbe/Dynamask/issues. show that, although our
method performs slightly lower than some baselines in terms of AUP, it
significantly outperforms all other methods by every other metric. In
particular, while it slightly outperforms DynaMask in terms of AUR, it yields
better results in terms of AUP, Information and Entropy. These results
therefore seem to indicate that using learnable perturbations should be
preferred compared with fixed one when explaining predictions based on
multivariate time series data.

#### Ablation study on the lambdas.

We perform here an ablation study to determine which values of $\lambda_{1}$
and $\lambda_{2}$ should be used in Equation 10. We therefore run our
experiment using various values of $\lambda_{1}$ and $\lambda_{2}$. We report
our results on Table 2.

| $\lambda_{1}$
---|---
| 0.01 | 0.1 | 1 | 10 | 100
$\lambda_{2}$ | 0.01 | 0.51 - 0.81 | 0.76 - 0.44 | 0.78 - 0.09 | 0.35 - 0.17 | 0.39 - 0.18
0.1 | 0.51 - 0.91 | 0.65 - 0.83 | 0.95 - 0.08 | 0.32 - 0.16 | 0.37 - 0.20
1 | 0.51 - 0.89 | 0.63 - 0.83 | 0.89 - 0.75 | 0.30 - 0.16 | 0.35 - 0.18
10 | 0.48 - 0.90 | 0.65 - 0.83 | 0.89 - 0.74 | 0.99 - 0.26 | 0.41 - 0.19
100 | 0.49 - 0.90 | 0.65 - 0.84 | 0.90 - 0.74 | 0.99 - 0.27 | 0.37 - 0.17

Table 2: Influence of $\lambda_{1}$ and $\lambda_{2}$ from Equation 10 on the
results of the HMM experiment. For each pair of parameters, 2 values are
reported: AUP - AUR. The average result over 5 runs is reported.

This table show that first $\lambda_{1}$ needs to be close to 1 to yield good
results. Indeed, a low value means lower regularisation, therefore retaining a
lot of unimportant features. A high value, on the other hand, forces m to be
mostly 0, yielding most features to be considered unimportant. Moreover,
$\lambda_{2}$ needs to be at least 1 to force NN(x) to learn uninformative
perturbations. Otherwise, there is only a weak mechanism to prevent NN from
producing an output similar to x.

#### Learning perturbations as a deletion game.

We also explore here learning perturbation using Equation 1, masking as little
data as possible while changing the model’s predictions as much as possible.
However, we cannot directly use Equation 1 for two reasons. First, the term
$-\mathcal{L}(f(\textbf{x}),f(\Phi(\textbf{x},\textbf{m})))$ is hard to
optimize, as it entices $f(\Phi(\textbf{x},\textbf{m}))$ to be “far” from
$f(\textbf{x})$ while there is no clarity on what “far” should be here. For
this reason, we replace this objective with
$\mathcal{L}(f(\textbf{0}),f(\Phi(\textbf{x},\textbf{m})))$, enticing the
predictions to be close to predictions made using 0, uninformative, as an
input. Second, we need to add the term $||\textrm{NN}(\textbf{x})||_{1}$ in
the loss. This results in the following objective:

$\operatorname*{arg\,min}_{\textbf{m},\Theta\in\textrm{NN}}\lambda_{1}||\textbf{1}-\textbf{m}||_{1}+\lambda_{2}||\textrm{NN}(\textbf{x})||_{1}+\mathcal{L}(f(\textbf{0}),f(\Phi(\textbf{x},\textbf{m})))$
(12)

We present our results on Table 3, comparing the preservation and the deletion
modes. While the second setting outperforms the first one in terms of AUR, it
performs poorly according to every other metrics. This might be due to the use
of $\mathcal{L}(f(\textbf{0}),f(\Phi(\textbf{x},\textbf{m})))$, which amounts
to learning a “change” in the predictions. This is a less straightforward
objective compared with the preservation mode, which aims to retain the
original predictions.

Mode | AUP $\uparrow$ | AUR $\uparrow$ | I $\uparrow$ | S $\downarrow$
---|---|---|---|---
Preservation |  0.885 (0.030) | 0.781 (0.013) |  1536 (79.0) |  34.1 (3,70)
Deletion | 0.346 (0.0034) |  0.863 (0.012) | 1079 (41.5) | 68.0 (5.07)

Table 3: Comparison of using the preservation mode vs deletion mode on the HMM
experiment. The average result over 5 runs is reported.

###  4.2 MIMIC-III experiment

We evaluate our method on the real-world MIMIC-III dataset, following the
works of Tonekaboni et al. (2020) and Crabbé & Van Der Schaar (2021). MIMIC-
III consists of patients in intensive-care units (ICU), for which a number of
vital signs and lab test results have been regularly measured. The task is
here to predict in-hospital mortality of each patient based on 48 hours of
data, discretised over each hour. Missing values are imputed using the
previous available ones. If there is no previous feature, a standard value is
imputed.

We train a one layer GRU with a hidden size of 200 to predict this in-hospital
mortality, and we aim to explain this model. In this dataset, the true salient
features are unknown, and we need to provide different metrics to evaluate our
method. Following Crabbé & Van Der Schaar (2021), we compare the original
predictions to ones where a certain proportion of the features have been
masked. We replace masked features either with an average over time of this
feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$, where T = 48 (hours)
or with zeros: $\overline{x}_{t,i}=0$. We use two metrics proposed by Crabbé &
Van Der Schaar (2021), and we also draw from the work of Shrikumar et al.
(2017) and DeYoung et al. (2019) and propose three additional metrics. These
resulting four metrics are then:

  * •

Accuracy (Acc): We mask the most salient features and compute the resulting
accuracy using this masked data. A lower accuracy means that important
features to make accurate predictions have been removed. Therefore, lower is
better with this metric.

  * •

Cross-Entropy (CE): We mask the most salient features and compute the cross-
entropy between predictions made with this masked data with the original one.
A higher value indicates that the predictions have more significantly changed
and that important features have been removed. Higher is better with this
metric.

  * •

Comprehensiveness (Comp): We mask the most salient features and compute the
average change of the predicted class probability compared with the original
one. Higher is better with this metric.

  * •

Sufficiency (Suff): We only keep the most salient features, and compute the
average change of the predicted class probability compared with the original
one. Lower is better with this metric.

Similar to our previous experiment, we use a bidirectional GRU as our
perturbation model. We compare our method against DeepLift (Shrikumar et al.,
2017), DynaMask (Crabbé & Van Der Schaar, 2021), Integrated Gradients (IG)
(Sundararajan et al., 2017), GradientShap (Lundberg & Lee, 2017), Lime
(Ribeiro et al., 2016), Augmented Occlusion (Tonekaboni et al., 2020),
Occlusion (Zeiler & Fergus, 2014) and Retain (Choi et al., 2016).

Method | Acc $\downarrow$ | Comp $\uparrow$ | CE $\uparrow$ | Suff $\downarrow$
---|---|---|---|---
DeepLift | 0.988 (0.002) | -4.36E-4 (0.001) | 0.097 (0.006) | 2.86E-3 (0.001)
DynaMask | 0.990 (0.001) | 2.21E-4 (0.001) | 0.097 (0.005) | 2.99E-3 (0.001)
IG | 0.988 (0.003) | 2.24E-4 (0.002) | 0.098 (0.006) | 2.21E-3 (0.001)
GradientShap | 0.987 (0.004) | -2.19E-3 (0.001) | 0.095 (0.006) | 3.99E-3 (0.001)
Lime | 0.996 (0.001) | -7.36E-4 (0.001) | 0.094 (0.005) | 3.39E-3 (0.001)
Occlusion | 0.988 (0.001) | -1.93E-3 (0.001) | 0.095 (0.005) | 4.57E-3 (0.001)
Aug Occlusion | 0.989 (0.001) | 4.59E-4 (0.001) | 0.098 (0.005) | 1.90E-3 (0.002)
Retain | 0.989 (0.001) | -3.79E-3 (0.001) | 0.093 (0.005) | 7.70E-3 (0.001)
Ours |  0.981 (0.004) |  1.53E-2 (0.004) |  0.118 (0.008) |  -1.19E-2 (0.004)

Table 4: Results of each explanation method compared with ours, by masking 20%
of the data and replacing masked features with an average over time:
$\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. For each metric, $\uparrow$
indicates that higher is better, and $\downarrow$ that lower is better. Mean
and std are reported over 5 folds.

Method | Acc $\downarrow$ | Comp $\uparrow$ | CE $\uparrow$ | Suff $\downarrow$
---|---|---|---|---
DeepLift | 0.972 (0.003) | -1.19E-3 (0.007) | 0.125 (0.014) | -6.92E-3 (0.006)
DynaMask | 0.975 (0.002) | -1.27E-3 (0.004) | 0.106 (0.009) | 6.57E-3 (0.012)
IG | 0.972 (0.003) | 1.24E-4 (0.007) | 0.127 (0.015) | -7.61E-3 (0.006)
GradientShap | 0.968 (0.006) | -6.28E-3 (0.004) | 0.128 (0.017) | 6.61E-4 (0.005)
Lime | 0.983 (0.003) | -5.22E-3 (0.004) | 0.093 (0.008) | -2.23E-3 (0.019)
Occlusion | 0.971 (0.003) | -4.03E-3 (0.003) | 0.122 (0.008) | -4.97E-3 (0.008)
Aug Occlusion | 0.972 (0.003) | -6.88E-4 (0.004) | 0.121 (0.009) | -4.62E-3 (0.011)
Retain | 0.971 (0.003) | -8.01E-3 (0.006) | 0.0123 (0.009) | 4.90E-4 (0.007)
Ours |  0.943 (0.008) |  1.09E-1 (0.023) |  0.318 (0.057) |  -6.94E-2 (0.006)

Table 5: Results of each explanation method compared with ours, by masking 20%
of the data and replacing masked features with zeros: $\overline{x}_{t,i}=0$.
For each metric, $\uparrow$ indicates that higher is better, and $\downarrow$
that lower is better. Mean and std are reported over 5 folds.

Figure 3: Cross Entropy replacing masked data with an average. We present here
the results in terms of cross-entropy by masking between 10% and 60% of the
data for each patient, and replacing the masked data with the overall average
over time for each feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$.
For clarity, we only plot a subset of the baselines. Higher is better with
this metric.

Figure 4: Cross Entropy replacing masked data with zeros. We present here the
results in terms of cross-entropy by masking between 10% and 60% of the data
for each patient, and replacing the masked data with zeros:
$\overline{x}_{t,i}=0$. For clarity, we only plot a subset of the baselines.
Higher is better with this metric.

We present on Tables 4 and 5 results with our method compared with different
baselines, computing our metrics by masking 20% of the data, and replacing
these features with either an average over time (Table 4) or zeros (Table 5).
We also plot on Figures 3 and 4 the cross-entropy (CE) metrics by masking
different proportion of the data, and replacing masked data with either an
average over time (Figure 3) or zeros (Figure 4). We also perform ablation
studies in Appendix A and provide more results in Appendix B.

Our results show that our method significantly outperforms every other method
on every metric, both using the average over time or zeros as masked data.
This also indicates that using learned perturbations is preferable to using
fixed ones when explaining predictions on multivariate time series data.

#### Choice of the perturbation generator.

While our method seems to perform well compared with existing baselines, we
want here to study the impact of the choice of NN in Equation 8. In this
study, we propose the following models:

  * •

Zero: NN(x) is set to zero everywhere. Equation 8 is then simply:
$\Phi(\textbf{x},\textbf{m})=\textbf{m}\times\textbf{x}$.

  * •

GRU: We use a one layer GRU model:
$\textrm{NN}(\textbf{x})=\textrm{GRU}(\textbf{x})$, which corresponds to a
generalisation of the fixed perturbation $\mu^{p}_{t,i}$ in Crabbé & Van Der
Schaar (2021).

  * •

Bi-GRU: Finally, we use a one layer bidirectional GRU
$\textrm{NN}(\textbf{x})=\textrm{bi-GRU}(\textbf{x})$, which corresponds to a
generalisation of the fixed perturbation $\mu_{t,i}$ in Crabbé & Van Der
Schaar (2021).

Method | Acc $\downarrow$ | Comp $\uparrow$ | CE $\uparrow$ | Suff $\downarrow$
---|---|---|---|---
Zeros | 0.981 (0.003) | 1.36E-2 (0.001) | 0.116 (0.004) | -1.02E-2 (0.002)
GRU |  0.980 (0.004) |  1.76E-2 (0.001) |  0.122 (0.004) |  -1.37E-2 (0.002)
Bi-GRU | 0.981 (0.004) | 1.53E-2 (0.004) | 0.118 (0.008) | -1.19E-2 (0.004)

Table 6: Comparison of different perturbation models, masking 20% of the data
and replacing masked features with an average over time:
$\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. For each metric, $\uparrow$
indicates that higher is better, and $\downarrow$ that lower is better. Mean
and std are reported over 5 folds.

Method | Acc $\downarrow$ | CE $\uparrow$ | Comp $\uparrow$ | Suff $\downarrow$
---|---|---|---|---
Zeros | 0.951 (0.005) | 9.64E-2 (0.013) | 0.305 (0.015) | -6.79E-2 (0.002)
GRU |  0.943 (0.007) |  1.22E-1 (0.008) |  0.344 (0.017) |  -7.40E-2 (0.001)
Bi-GRU |  0.943 (0.008) | 1.09E-1 (0.023) | 0.318 (0.057) | -6.94E-2 (0.006)

Table 7: Comparison of different perturbation models, masking 20% of the data
and replacing masked features with zeros: $\overline{x}_{t,i}=0$. For each
metric, $\uparrow$ indicates that higher is better, and $\downarrow$ that
lower is better. Mean and std are reported over 5 folds.

We present our results on MIMIC-III on Tables 6 and 7, replacing 20% of the
data with either an overall average of each feature over time (Table 6), or
with zeros (Table 7). We use the same metrics as with our main MIMIC-III
experiments. As with the main experiment, we provide more results, masking
different proportions of the data, in Appendix B.

Our results are interesting on several accounts. Firstly, the Zeros method,
which simply perturbs the data by masking non salient features:
$\Phi(\textbf{x},\textbf{m})=\textbf{m}\times\textbf{x}$, performs
significantly better than all other baselines, including DynaMask with fixed
perturbations. As each measure in our dataset is normalised, masking one
measure with the Zeros method amounts to replacing it with its average over
the entire dataset. On the other hand, DynaMask replaced mask data with its
average over time for each individual patient. This good performance of Zeros
could be therefore explained by the fact that many measures do not vary much
over time. As a result, replacing masked data with an overall average would be
much more informative than replacing it with an average over time for each
patient.

Secondly, while using the bidirectional GRU perturbation yields better results
than Zeros, it is itself outperformed by our method with the unidirectional
GRU perturbation. Moreover, using this unidirectional GRU also yields more
stable results with a lower standard deviation. Our intuition was that a
bidirectional GRU would yield better results, as it would be able to produce
outputs based on past and future events. However, it seems that modeling
perturbations ignoring future events seems to yield better and more stable
results. We used a Bi-GRU to produce our results in Tables 4 and 5, as it
corresponds to our original intuition, but we also recommend testing different
types of neural networks for best performance when applying our method.

#### Analysis of salient features.

We present on Figure 5 the most salient features, averaged over every positive
patient, to determine which factors are most important when determining in-
hospital mortality.

Figure 5: Importance of each feature to predict in-hospital mortality. For
each feature, we present its average importance over time and over multiple
patients, using our method with a GRU perturbation network. We infer from
these results that anion gap, bicarbonate level, platelet count, systolic
blood pressure and respiratory rate are most important for our model when
making a prediction. We also plot a 95% confidence interval around these
averages.

This averaged feature importance indicates a few salient features: anion gap,
bicarbonate level, platelet count, systolic blood pressure and respiratory
rate. This seems to be consistent with the literature, which has highlighted
the importance of these features, conducting studies on the saliency of
bicarbonate levels (Lim et al., 2019), platelet count (Zhang et al., 2014) and
systolic blood pressure (Kondo et al., 2011). The influence of anion gap on
in-hospital mortality is less clear, with conflicting studies on this subject
(Glasmacher & Stones, 2015). On the other hand, the respiratory rate is often
neglected despite being an important predictor of serious events (Cretikos et
al., 2008).

However, although the 95% confidence interval associated with these features
importances is small due to a large number of patients, there remains a large
corresponding standard deviation. We can therefore infer that the importance
of each feature greatly depends on each patient. It is indeed possible that a
measure such as the systolic blood pressure, for instance, only matters when
it is outside of a normal range. As a result, its importance will greatly vary
depending on each patient’s condition. This demonstrates the superiority of
perturbation-based methods compared to directly using a simpler interpretable
model such as a decision tree instead of a neural network to predict in-
hospital mortality. Indeed, such methods can only infer feature importance on
average, and not explain each prediction individually.

Figure 6: Average feature importance over time to predict in-hospital
mortality. We average every feature’s importance on all positive patients,
over each of the 48 hours of measurement in hospital. We display the mean with
a 95% confidence interval. Our results show that later measures are most
important predictor of in-hospital mortality.

Figure 7: Positive rate over positive patients by masking first or last data
points. We compare the influence of the first or last data points in time by
masking them successively. We observe that masking latter points in time
yields a much lower positive rate than masking former ones. Masking the last
12 measures indeed results in every positive patient predicted as negative.

In addition to determining which feature is salient, our method can also infer
when it is salient. As a result, we present on Figure 6 the average over
positive patients of the importance of all features at each hour. This figure
shows that later measurements have a larger impact on the outcome compared
with earlier data. To evaluate the accuracy of this finding, we also plot on
Figure 7 the positive rate over (true or false) positive patients, when
masking earlier measures on one hand, and later measures on the other hand. We
can see that masking early features has a minimal impact on the predictions,
while masking late features has on the contrary a dramatic impact on the
outcome. As a result, it seems that, when predicting in-hospital mortality,
the last measurements of each patient is more important to make a prediction,
rather than the overall evolution of the patient.

##  5 Conclusion

In this work, we have presented an extension of Fong & Vedaldi (2017) and
Crabbé & Van Der Schaar (2021) to better explain multivariate time series
predictions using a perturbation-based saliency method. Our main intuition is
that the choice made by Crabbé & Van Der Schaar (2021) of fixed perturbations
is less adapted to temporal data due to the possibility of long-term
dependencies.

Our results show that using learned perturbations yields better explanations
compared with existing methods, including the DynaMask one with fixed
perturbations. We have also studied the choice of the neural network to model
the perturbation and found that, on the in-hospital mortality task of MIMIC-
III, a unidirectional GRU yields better and more stable results than the
bidirectional one.

Using our method, we have also been able to derive some insights on the neural
network predicting in-hospital mortality: which features are on average most
important, as well as which measures in time. Precise temporal attributions
could be derived similarly for each patient, giving further insight on this
model’s behavior.

Moreover, an inherent limitation to perturbation-based methods such as ours is
that it is not able to specify the direction of an explanation. As such, it
can measure if a specific feature is important, but cannot distinguish between
features having a positive or a negative influence on the prediction. Adapting
our method to tackle this issue would prove very beneficial for applications
in healthcare.

##  6 Acknowledgements

The author would like to thank Vitalii Zhelezniak for his insightful remarks
and recommendations during the elaboration of this work. We also thank Anthony
Hu and Thomas Uriot for their detailed initial reviews of this paper.

## References

  * Ahuja (2019) Ahuja, A. S.  The impact of artificial intelligence in medicine on the future role of the physician.  _PeerJ_ , 7:e7702, 2019.
  * Cho et al. (2014) Cho, K., Van Merriënboer, B., Bahdanau, D., and Bengio, Y.  On the properties of neural machine translation: Encoder-decoder approaches.  _arXiv preprint arXiv:1409.1259_ , 2014.
  * Choi et al. (2016) Choi, E., Bahadori, M. T., Sun, J., Kulas, J., Schuetz, A., and Stewart, W.  Retain: An interpretable predictive model for healthcare using reverse time attention mechanism.  _Advances in neural information processing systems_ , 29, 2016.
  * Crabbé & Van Der Schaar (2021) Crabbé, J. and Van Der Schaar, M.  Explaining time series predictions with dynamic masks.  In _International Conference on Machine Learning_ , pp. 2166–2177. PMLR, 2021.
  * Cretikos et al. (2008) Cretikos, M. A., Bellomo, R., Hillman, K., Chen, J., Finfer, S., and Flabouris, A.  Respiratory rate: the neglected vital sign.  _Medical Journal of Australia_ , 188(11):657–659, 2008.
  * DeYoung et al. (2019) DeYoung, J., Jain, S., Rajani, N. F., Lehman, E., Xiong, C., Socher, R., and Wallace, B. C.  Eraser: A benchmark to evaluate rationalized nlp models.  _arXiv preprint arXiv:1911.03429_ , 2019.
  * Fong et al. (2019) Fong, R., Patrick, M., and Vedaldi, A.  Understanding deep networks via extremal perturbations and smooth masks.  In _Proceedings of the IEEE/CVF international conference on computer vision_ , pp. 2950–2958, 2019.
  * Fong & Vedaldi (2017) Fong, R. C. and Vedaldi, A.  Interpretable explanations of black boxes by meaningful perturbation.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 3429–3437, 2017.
  * Glasmacher & Stones (2015) Glasmacher, S. A. and Stones, W.  Anion gap as a prognostic tool for risk stratification in critically ill patients–a systematic review and meta-analysis.  _BMC anesthesiology_ , 16(1):1–13, 2015.
  * Johnson et al. (2016) Johnson, A. E., Pollard, T. J., Shen, L., Lehman, L.-w. H., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., Anthony Celi, L., and Mark, R. G.  Mimic-iii, a freely accessible critical care database.  _Scientific data_ , 3(1):1–9, 2016.
  * Kondo et al. (2011) Kondo, Y., Abe, T., Kohshi, K., Tokuda, Y., Cook, E. F., and Kukita, I.  Revised trauma scoring system to predict in-hospital mortality in the emergency department: Glasgow coma scale, age, and systolic blood pressure score.  _Critical care_ , 15(4):1–8, 2011.
  * LaRosa & Danks (2018) LaRosa, E. and Danks, D.  Impacts on trust of healthcare ai.  In _Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society_ , pp. 210–215, 2018.
  * Lim et al. (2019) Lim, S. Y., Park, Y., Chin, H. J., Na, K. Y., Chae, D.-W., and Kim, S.  Short-term and long-term effects of low serum bicarbonate level at admission in hospitalised patients.  _Scientific reports_ , 9(1):1–7, 2019.
  * Lundberg & Lee (2017) Lundberg, S. M. and Lee, S.-I.  A unified approach to interpreting model predictions.  _Advances in neural information processing systems_ , 30, 2017.
  * Moscato et al. (2021) Moscato, V., Picariello, A., and Sperlí, G.  A benchmark of machine learning approaches for credit score prediction.  _Expert Systems with Applications_ , 165:113986, 2021.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C.  ” why should i trust you?” explaining the predictions of any classifier.  In _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , pp. 1135–1144, 2016.
  * Shrikumar et al. (2017) Shrikumar, A., Greenside, P., and Kundaje, A.  Learning important features through propagating activation differences.  In _International conference on machine learning_ , pp. 3145–3153. PMLR, 2017.
  * Simonyan et al. (2013) Simonyan, K., Vedaldi, A., and Zisserman, A.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  _arXiv preprint arXiv:1312.6034_ , 2013.
  * Sundararajan et al. (2017) Sundararajan, M., Taly, A., and Yan, Q.  Axiomatic attribution for deep networks.  In _International conference on machine learning_ , pp. 3319–3328. PMLR, 2017.
  * Tjoa & Guan (2020) Tjoa, E. and Guan, C.  A survey on explainable artificial intelligence (xai): Toward medical xai.  _IEEE transactions on neural networks and learning systems_ , 32(11):4793–4813, 2020.
  * Tonekaboni et al. (2020) Tonekaboni, S., Joshi, S., Campbell, K., Duvenaud, D. K., and Goldenberg, A.  What went wrong and when? instance-wise feature importance for time-series black-box models.  _Advances in Neural Information Processing Systems_ , 33:799–809, 2020.
  * Zeiler & Fergus (2014) Zeiler, M. D. and Fergus, R.  Visualizing and understanding convolutional networks.  In _European conference on computer vision_ , pp. 818–833. Springer, 2014.
  * Zhang et al. (2014) Zhang, Z., Xu, X., Ni, H., and Deng, H.  Platelet indices are novel predictors of hospital mortality in intensive care unit patients.  _Journal of critical care_ , 29(5):885–e1, 2014\.

##  Appendix A Addition studies using the Mimic3 dataset

First, we perform here an ablation study to determine which values of
$\lambda_{1}$ and $\lambda_{2}$ should be used in Equation 10 on the Mimic3
dataset, similarly to the one done on HMM in Section 4. We run our experiment
using various values of $\lambda_{1}$ and $\lambda_{2}$, running on 5
different seed for each pair of parameters. We report our results on Table 8.

This table shows similarly that a high value of $\lambda_{2}$ should be
chosen, to force NN(x) to learn uninformative perturbations. Interestingly,
using a low value of $\lambda_{1}$ yields stronger results on Mimic3,
indicating that constraining too much m to be close to 0 is not necessarily a
good option. However, this setting does not yield good results on the HMM
dataset, therefore our method should be carefully evaluated when used on low
values of $\lambda_{1}$.

| $\lambda_{1}$
---|---
| 0.01 | 0.1 | 1 | 10 | 100
$\lambda_{2}$ | 0.01 | 0.926 - 0.348 | 0.965 - 0.160 | 0.968 - 0.186 | 0.968 - 0.165 | 0.956 - 0.169
0.1 | 0.893 - 0.503 | 0.935 - 0.328 | 0.961 - 0.166 | 0.961 - 0.261 | 0.936 - 0.271
1 | 0.881- 0.534 | 0.899 - 0.491 | 0.947 - 0.284 | 0.957 - 0.191 | 0.960 - 0.261
10 | 0.857 - 0.540 | 0.905 - 0.480 | 0.940 - 0.372 | 0.950 - 0.266 | 0.964 - 0.163
100 | 0.862 - 0.542 | 0.910 - 0.472 | 0.950 - 0.362 | 0.946 - 0.289 | 0.956 - 0.173
Table 8: Influence of $\lambda_{1}$ and $\lambda_{2}$ from Equation 10 on the
results of the Mimic3 experiment. or each pair of parameters, 2 values are
reported: Accuracy - Cross-entropy. For accuracy, lower is better, while
higher is better for cross-entropy. Each metric is computed by masking 20 % of
the data and replacing masked features with zeros: $\overline{x}_{t,i}=0$. The
average result over 5 runs is reported.

Second, we learn perturbations as a deletion game, similarly to the experiment
conducted on the HMM dataset in Section 4. As such, we use Equation 12 in this
experiment. We report our results on Table 9.

Similarly to the HMM experiment, this table shows that the deletion mode
performs poorly compared with the preservation one. The latter mode should
therefore be preferred to the former.

Mode | Acc $\downarrow$ | Comp $\uparrow$ | CE $\uparrow$ | Suff $\downarrow$
---|---|---|---|---
Preservation |  0.943 (0.008) |  1.09E-1 (0.023) |  0.318 (0.057) |  -6.94E-2 (0.006)
Deletion | 0.977 (0.003) | -0.025 (0.009) | 0.079 (0.012) | 0.053 (0.013)
Table 9: Comparison of using the preservation mode vs deletion mode on the
Mimic3 experiment. The average result over 5 runs is reported.

##  Appendix B Additional results on the in-hospital mortality experiment

We present below more results on the in-hospital mortality experiment, based
on the MIMIC-III dataset. Results in terms of accuracy, comprehensiveness and
sufficiency can be found on Figures 9, 9, 11, 11, 13 and 13.

We also provide here more results of the ablation study, comparing using
Zeros, a GRU or a BiGRU as a perturbation model. Results in terms of accuracy,
cross-entropy, comprehensiveness and sufficiency can be found on Figures 15,
15, 19, 19, 21 and 21.

Figure 8: Accuracy, masking between 10% and 60% of the data for each patient,
and replacing the masked data with the overall average over time for each
feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. For clarity, we only
plot a subset of the baselines. Lower is better with this metric.

Figure 9: Accuracy, masking between 10% and 60% of the data for each patient,
and replacing the masked data with zeros: $\overline{x}_{t,i}=0$. For clarity,
we only plot a subset of the baselines. Lower is better with this metric.

Figure 10: Comprehensiveness, masking between 10% and 60% of the data for each
patient, and replacing the masked data with the overall average over time for
each feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. For clarity, we
only plot a subset of the baselines. Higher is better with this metric.

Figure 11: Comprehensiveness, masking between 10% and 60% of the data for each
patient, and replacing the masked data with zeros: $\overline{x}_{t,i}=0$. For
clarity, we only plot a subset of the baselines. Higher is better with this
metric.

Figure 12: Sufficiency, masking between 10% and 60% of the data for each
patient, and replacing the masked data with the overall average over time for
each feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. For clarity, we
only plot a subset of the baselines. Lower is better with this metric.

Figure 13: Sufficiency, masking between 10% and 60% of the data for each
patient, and replacing the masked data with zeros: $\overline{x}_{t,i}=0$. For
clarity, we only plot a subset of the baselines. Lower is better with this
metric.

Figure 14: Accuracy results, masking between 10% and 60% of the data for each
patient, and replacing the masked data with the overall average over time for
each feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. Lower is better
with this metric.

Figure 15: Accuracy results, masking between 10% and 60% of the data for each
patient, and replacing the masked data with zeros: $\overline{x}_{t,i}=0$.
Lower is better with this metric.

Figure 16: Cross-entropy results, masking between 10% and 60% of the data for
each patient, and replacing the masked data with the overall average over time
for each feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. Higher is
better with this metric.

Figure 17: Cross-entropy results, masking between 10% and 60% of the data for
each patient, and replacing the masked data with zeros:
$\overline{x}_{t,i}=0$. Higher is better with this metric.

Figure 18: Comprehensiveness results, masking between 10% and 60% of the data
for each patient, and replacing the masked data with the overall average over
time for each feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. Higher
is better with this metric.

Figure 19: Comprehensiveness results, masking between 10% and 60% of the data
for each patient, and replacing the masked data with zeros:
$\overline{x}_{t,i}=0$. Higher is better with this metric.

Figure 20: Sufficiency results, masking between 10% and 60% of the data for
each patient, and replacing the masked data with the overall average over time
for each feature: $\overline{x}_{t,i}=\frac{1}{T}\sum_{t}x_{t,i}$. Lower is
better with this metric.

Figure 21: Sufficiency results, masking between 10% and 60% of the data for
each patient, and replacing the masked data with zeros:
$\overline{x}_{t,i}=0$. Lower is better with this metric.


# Gaussian process error modeling for chiral effective-field-theory
calculations of $np\leftrightarrow d\gamma$ at low energies

B. Acharya  Institut für Kernphysik and PRISMA+ Cluster of Excellence,
Johannes Gutenberg-Universität Mainz, 55128 Mainz, Germany  S. Bacca  Institut
für Kernphysik and PRISMA+ Cluster of Excellence, Johannes Gutenberg-
Universität Mainz, 55128 Mainz, Germany  Helmholtz-Institut Mainz, Johannes
Gutenberg-Universität Mainz, D-55099 Mainz, Germany

###### Abstract

We calculate the energy-dependent cross section of the $np\leftrightarrow
d\gamma$ process in chiral effective field theory and apply state-of-the-art
tools for quantification of theory uncertainty. We focus on the low-energy
regime, where the magnetic dipole and the electric dipole transitions cross
over, including the range relevant for big-bang nucleosynthesis. Working with
the leading one- and two-body electromagnetic currents, we study the order-by-
order convergence of this observable in the chiral expansion of the nuclear
potential. We find that the Gaussian process error model describes the
observed convergence very well, allowing us to present Bayesian credible
intervals for the truncation error with correlations between the cross
sections at different energies taken into account. We obtain a 1$\sigma$
estimate of about 0.2% for the uncertainty from the truncation of the nuclear
potential. This is an important step towards calculations with statistically
interpretable uncertainties for astrophysical reactions involving light
nuclei.

Bayesian analysis, uncertainty quantification, big bang nucleosynthesis,
Gaussian process, machine learning

Over the past couple of decades, effective field theories (EFTs) have emerged
as powerful tools to study nuclear structure and reactions within systematic
frameworks where a well-justified quantification of theory uncertainties is
feasible Epelbaum:2008ga . An EFT is a model-independent description of
interacting particles at some low-momentum scale $Q$ in terms of the minimal
degrees of freedom relevant below a breakdown momentum scale, $\Lambda$.
Physical observables are calculated as systematic expansions in $Q/\Lambda$,
with undetermined parameters called low-energy constants (LECs). The EFT can
predict an observable at a given order once all the LECs that appear up that
order have been fixed, _e.g.,_ by fitting to experimental data for other
observables. The truncation of the EFT expansion at a finite order results in
a theory error that can be estimated. Recent adoption of Bayesian methods,
which have found increasing prominence in nuclear theory Schindler:2008fh ;
Furnstahl:2014xsa ; Furnstahl:2015rha ; Wesolowski:2015fqa ; Zhang:2015ajn ;
Melendez:2017phj ; Wesolowski:2018lzj ; Lovell:2018bao ; Neufcourt:2018syo ;
Zhou:2018bon ; King:2019sax ; Melendez:2019izc ; Premarathna:2019tup ;
Zhang:2019odg ; Epelbaum:2019zqc ; Phillips:2020dmw ; Volkotrub:2020lsr ;
Reinert:2020mcu ; Wesolowski:2021cni ; Rijneveen:2020qbc ; Maris:2020qne ;
Wellenhofer:2021eis ; Connell:2021qcd ; Rijneveen:2021bfw , have enabled a
statistically meaningful estimation of truncation errors which can be easily
combined with other uncertainties that arise from experimental data,
calculational methods, and parameter fitting (see, _e.g.,_ Ref.
Wesolowski:2021cni ). Such a complete and consistent accounting of theory
uncertainties is of paramount importance not only when the theoretical results
have to be compared to precise experiments, but also when theoretical
predictions need to be used as proxies for experimental data that are
imprecise or unavailable. The latter situation is common for nuclear reactions
in astrophysics and cosmology, for example in those relevant to big-bang
nucleosynthesis (BBN).

There is, in general, a good agreement between simulations of nuclear reaction
networks and astronomical observations for abundances of light nuclei that
were produced during BBN. This agreement stands not only as a historic
cornerstone of big-bang cosmology, but also as a crucial test for future
extensions to the Standard Model. The ability of BBN to constrain key
quantities such as cosmic baryon density rests heavily on the primordial
Deuterium to Hydrogen abundance ratio, whose uncertainty is currently
dominated by the uncertainties of the Deuterium-burning processes:
2H$(p,\gamma){}^{3}$He, 2H$(d,n){}^{3}$He and 2H$(d,p){}^{3}$H. Improvements
in precision on their rates, exemplified by the recent 3% measurement of the
2H$(p,\gamma){}^{3}$He cross section by the LUNA collaboration Mossa:2020gjc ,
will eventually call for an update of the primordial Deuterium production
reaction, $np\rightarrow d\gamma$. This reaction is the first step in the BBN
network and marks the end of the so-called Deuterium bottleneck. BBN
simulations Serpico_2004 ; Pisanti:2007hk ; Coc:2019rza ; Pitrou:2018cgg  rely
on theoretical predictions for this important cross section because available
experimental data are too sparse and not sufficiently precise in the relevant
energy regime. It is therefore imperative to study it using different
theoretical approaches and to rigorously quantify uncertainties. Furthermore,
ever increasing precision of the constraints on the primordial Deuterium
abundance from cosmic microwave background analysis Planck:2018vyg  and
observational astronomy of low-metallicity stars Cooke:2017cwo  need to be
independently matched by nuclear physics determinations.

This reaction has traditionally been studied using phenomenological models for
nuclear interactions and currents Arenhovel:1990yg ; Carlson:1997qn ;
Marcucci:2004sq . Pionless EFT, the low-energy EFT of quantum chromodynamics
(QCD) that uses nucleons as the only dynamical degrees of freedom and has a
breakdown scale $\Lambda$ at the pion mass $m_{\pi}$, was applied to
$np\rightarrow d\gamma$ at BBN energies in Refs. chen-savage ; rupak ; ando .
By fitting an LEC in the $E1$ amplitude to photodissociation data and an LEC
in the $M1$ amplitude to $p(n,\gamma)d$ measurement cokinos-melkonian  at
energy $E=1.2625\times 10^{-8}~{}\mathrm{MeV}$ in the $np$ center-of-mass
frame, Ref. rupak  quoted a sub-percentage precision for the cross section at
BBN energies. The pionless EFT results of Refs. rupak ; ando , along with the
accompanying uncertainties, are adopted by modern BBN simulations
Pisanti:2007hk ; Pitrou:2018cgg .

This cross section has also been calculated at $E=1.2625\times
10^{-8}~{}\mathrm{MeV}$ in Ref. Beane:2015yha  using lattice QCD, in Ref.
piarulli  using chiral EFT ($\chi$EFT) and in Ref. song-lazauskas-park  using
a hybrid approach that utilized $\chi$EFT currents and phenomenological
potentials. Based on the scale hierarchy $Q\sim m_{\pi}<\Lambda$, $\chi$EFT
provides a description of hadronic and nuclear phenomena in terms of nucleons
and pions as effective degrees of freedom, and is therefore well suited for
investigating this cross section at energies beyond the thermal capture,
including the cosmological [20 keV–200 keV] range PhysRevLett.82.4176  and
higher.

In this work, we present the first $\chi$EFT results for the energy-dependent
$np\leftrightarrow d\gamma$ cross section for a range of energies from
threshold to MeV scale, that encompasses the BBN regime. We extend the
Bayesian procedure for estimating the $\chi$EFT truncation error, which was
earlier applied to various other observables—nucleon-mass expansion
Schindler:2008fh ; Wesolowski:2015fqa , nucleon-nucleon elastic scattering
Furnstahl:2015rha ; Melendez:2017phj ; Wesolowski:2018lzj ; Melendez:2019izc ;
Reinert:2020mcu , nucleon-deuteron scattering Epelbaum:2019zqc ;
Volkotrub:2020lsr , nuclear-matter equation of state PhysRevC.102.054315 ;
PhysRevLett.125.202702 , pion-photoproduction on the nucleon Rijneveen:2020qbc
; Rijneveen:2021bfw , and properties of light nuclei Maris:2020qne ;
Wesolowski:2021cni —to an electromagnetic reaction cross section. We adopt the
error model developed by Melendez _et al._ Melendez:2019izc  that used machine
learning with Gaussian processes (GPs) calibrated by physics-based
hyperparameters to determine the convergence of EFT predictions that may be
correlated across independent variables (_e.g.,_ energy). It allows us to
quantify, for the first time, the error in the $np\rightarrow d\gamma$ cross
section, $\sigma_{np}(E)$, and in that of the reverse process, $\sigma_{\gamma
d}(E)$, from truncation of the $\chi$EFT potential as a Bayesian degree-of-
belief band.

We begin by presenting $\chi$EFT results for the $np\leftrightarrow d\gamma$
cross section and related observables in Sec. I. We then briefly introduce the
Gaussian process error model in Sec. II. The fitted model, diagnostic checks
and predicted uncertainties are presented in Sec. III. We close with a brief
summary and outlook in Sec. IV.

##  I The $np\leftrightarrow d\gamma$ cross section in $\chi$EFT

The detailed-balance principle relates the $np\rightarrow d\gamma$ cross
section $\sigma_{np}$ to the deuteron photodissociation $d\gamma\rightarrow
np$ cross section $\sigma_{\gamma d}$:

$\sigma_{np}=\frac{3}{2}\,\frac{(s-m_{d}^{2})^{2}}{(s-4m^{2})(s-\delta
m^{2})}\,\sigma_{\gamma d}\,,$ (1)

where $m_{d}$ is the deuteron mass, $m$ is the isospin-averaged nucleon mass,
and $\delta m$ is the neutron-proton mass difference. The Mandelstam variable
$s$ can be conveniently expressed in terms of the $np$ relative energy $E$, or
the neutron energy $E_{n}$ in the rest frame of the proton,

$s=(E+2m)^{2}=2m_{p}E_{n}+4m^{2}\,,$ (2)

where $m_{p}$ is the proton mass. The photodissociation cross section can be
expressed as a function of the photon energy $\nu$ (in the rest frame of the
deuteron) through

$\sigma_{\gamma d}(\nu)=\frac{2\pi^{2}}{\nu}\,\alpha\,R_{T}(\nu,\nu)$ (3)

in terms of the transverse response function of the deuteron, defined as

$\displaystyle
R_{T}(\nu,q)=\frac{1}{6}\sum_{M}\sum_{S^{\prime}S_{z}^{\prime}}$
$\displaystyle\sum_{T^{\prime}}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\,\delta(\nu+m_{d}-E_{+}-E_{-})$
$\displaystyle\sum_{\lambda=\pm
1}|\langle\phi_{\mathbf{k},S^{\prime}S_{z}^{\prime},T^{\prime}}|j_{\lambda}|\psi_{M}\rangle|^{2}\,.$
(4)

Here, $j_{\lambda}$ are the spherical components of the electromagnetic
current operator ${\bf j}$ at four-momentum transfer $(\nu,\mathbf{q})$. The
deuteron ground state is denoted by $|\psi_{M}\rangle$, where $M$ is the
projection of the total angular momentum, whereas
$|\phi_{\mathbf{k},S^{\prime}S_{z}^{\prime},T^{\prime}}\rangle$ denotes the
$pn$ scattering state with the relative momentum, total isospin, total spin
and spin projection given by $\mathbf{k}$, $T^{\prime}$, $S^{\prime}$ and
$S_{z}^{\prime}$ respectively. $E_{\pm}$ are the energies of the final-state
nucleons in the rest frame of the deuteron, _i.e.,_
$E_{\pm}=\sqrt{(\mathbf{q}/2\pm\mathbf{k})^{2}+m^{2}}$. The deuteron and the
$np$ scattering-state wave functions are obtained by solving their Lippmann-
Schwinger equations using the momentum-space regularized semilocal potentials
of Ref. Reinert:2017usi  at various orders in the $\chi$EFT expansion. The
electromagnetic currents were first derived in $\chi$EFT in Refs. Park:1995pn
; Pastore:2008ui ; Kolling:2009iq . We use multipole expansions
Bijaya_neutrino ; chi-disp-mu-d  of these operators bibnote1  that contribute
at orders $(Q/\Lambda)^{-3,-2,1,0}$, namely the one-body currents consisting
of the convection and spin-magnetization terms, and the leading one-pion
exchange currents. Although we retain the quadrupole and octupole operators in
the multipole expansions, their contributions are negligible compared to the
dominant $M1$ and $E1$ operators at low energies. To compare the relative
strengths of the $M1$ and $E1$ transitions, we calculate the analyzing power,
defined as
$\Sigma(\theta)\equiv(N_{\parallel}-N_{\perp})/(N_{\parallel}+N_{\perp})$,
where $N_{\parallel}(N_{\perp})$ is the number of outgoing neutrons parallel
(perpendicular) to the photon-polarization plane in a photodissociation
experiment. It is related to the $M1$ and $E1$ contributions to the
photodissociation cross section through ando

$\Sigma(\theta)=\frac{\frac{3}{2}\sigma_{\gamma
d}^{E1}\sin^{2}\theta}{\sigma_{\gamma d}^{M1}+\frac{3}{2}\sigma_{\gamma
d}^{E1}\sin^{2}\theta}\,,$ (5)

where $\theta$ is the angle of the neutrons with respect to the photon beam
axis.

We now present the $\chi$EFT predictions for observables related to
$np\leftrightarrow d\gamma$ starting with the $p(n,\gamma)d$ cross section for
thermal neutrons in Table 1. Here (and throughout this work), we denote cross
sections calculated with an $n$th-order (NnLO) potential of Ref.
Reinert:2017usi  (and electromagnetic currents fixed to one-body plus two-body
one-pion exchange) as $y_{n}$. The contribution $y_{1}$ is zero because there
are no corrections to the $\chi$EFT potential at this order. At this threshold
energy, our results undershoot the experimental value of Ref. cokinos-
melkonian  by a few percent. It was shown in Ref. piarulli , that percentage-
level agreement with the experiment of Ref. cokinos-melkonian  required
current operators at the order $(Q/\Lambda)^{1}$ not considered in this work;
these included leading two-pion exchange, sub-leading corrections to one-pion
exchange and contact operators. Since the goal of this paper is to apply
machine learning tools for the first time to uncertainty quantification in
nuclear electroweak reactions in $\chi$EFT and including these operators
significantly complicates the calculations, we work with fixed currents and
focus on analyzing convergence in the $\chi$EFT potential only, leaving the
inclusion of higher order currents to future work. We discuss the implications
of our findings for such a study further below.

$y_{0}$ |  $y_{2}$ |  $y_{3}$ |  $y_{4}$ |  $y_{5}$ | Experiment
---|---|---|---|---|---
340.6 | 325.0 | 321.8 | 321.0 | 322.6 | 332.6 $\pm$ 0.7
Table 1:  The $p(n,\gamma)d$ cross section $\sigma_{np}$ in mb for thermal
neutrons (corresponding to $E=1.2625\times 10^{-08}$ MeV) calculated using
$\chi$EFT potentials at various orders. The experimental result is from Ref.
cokinos-melkonian .  Figure 1: $\chi$EFT predictions for the
$np\leftrightarrow d\gamma$ cross section and related observables at BBN
energies. (a) The product of $p(n,\gamma)d$ cross section and the neutron
speed versus the energy of the neutron. (b) The deuteron photodissociation
cross section as a function of the photon energy in the rest frame of the
deuteron. (c) The photon analyzing power for photodissociation versus its
energy. Experimental data are from Refs. suzuki  (triangles), nagai  (circle),
hara  (crosses), moreh  (square) and tornow  (diamonds). Beam-energy
resolution errors are not shown. Purple hexagons are calculated using pionless
EFT results of Ref. rupak .

In Fig. 1, we compare our $\chi$EFT predictions at different orders for
various observables related to $np\leftrightarrow d\gamma$ with the pionless
EFT results of Ref. rupak  and available experimental data at energies
relevant for BBN and beyond. The $\chi$EFT predictions agree with the
experimental data for all orders. Since Ref. rupak  uses the experimental
cross section given in Table 1 as input, their $M1$ contribution is larger
than our $\chi$EFT predictions at low energies.

##  II Gaussian process model for correlated EFT errors

Following the formalism introduced by Melendez _et al._ Melendez:2019izc , we
consider an observable $y$ as a function of the kinematic variable $p$, which,
in our case is the $np$ relative momentum. The order-$k$ EFT prediction is
written as

$y_{k}(p)=y_{\mathrm{ref}}(p)\sum_{n=0}^{k}c_{n}(p)\,[Q(p)/\Lambda]^{n}\,,$
(6)

and its EFT truncation error as

$\delta
y_{k}(p)=y_{\mathrm{ref}}(p)\sum_{n=k+1}^{\infty}c_{n}(p)\,[Q(p)/\Lambda]^{n}\,.$
(7)

Here $y_{\mathrm{ref}}(p)$ is a dimensionful quantity that sets the overall
scale such that the _a priori_ unknown dimensionless coefficients $c_{n}(p)$
are smooth naturally-sized (_i.e.,_ of order 1) curves, provided that the EFT
is converging in the expected manner. The Bayesian model for EFT error
quantification seeks to build upon this prior knowledge with the data, the
values of $c_{n\leq k}(p)$ evaluated from order-by-order EFT calculations up
through order $k$, to refine our expectations for $c_{n>k}(p)$, and thereby
obtain an estimate for $\delta y_{k}(p)$. Refs. Furnstahl:2014xsa ;
Furnstahl:2015rha ; Wesolowski:2015fqa ; Melendez:2017phj ; Wesolowski:2018lzj
; Melendez:2019izc  developed a pointwise error model that is applicable when
(i) $y$ is a number and not a function of one or more kinematic variables;
(ii) $y(p)$ is analyzed at only one value of $p$ or at multiple values of $p$
that are sufficiently far apart such that the $y(p)$ values can be safely
assumed to be uncorrelated. Ref. Melendez:2019izc  extended this framework to
study curvewise convergence by employing GPs. Below we first introduce the GP
error model and then present application to $np\leftrightarrow d\gamma$.

### The GP error model—

The basic idea of the error model is to use GPs to build a stochastic
representation (an _emulator_) to serve as a surrogate for the sequence of
deterministic calculations (the _simulator_) that yields order-by-order EFT
predictions up through order $k$. The statistical properties of the emulator
are then exploited to yield estimates for the EFT truncation errors at various
orders, subsuming the contributions of even those terms that have never been
calculated. Specifically, we start with the assumption that $c_{n}(p)$ are
independent draws from an underlying GP, _i.e.,_ they follow a multivariate
normal distribution for every finite set of $p$. The GP is completely
specified by the mean $\mu$ and covariance function
$\bar{c}^{2}r(p,p^{\prime};\ell)$,

$c_{n}(p)\,|\,\bar{c}^{2},\ell\,\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\,\mathcal{GP}[\mu,\bar{c}^{2}r(p,p^{\prime};\ell)]\,.$
(8)

The correlation function $r(p,p^{\prime};\ell)$ is commonly chosen to have a
squared-exponential form,

$r(p,p^{\prime};\ell)=\exp\left[-\frac{(p-p^{\prime})^{2}}{2\ell^{2}}\right]\,.$
(9)

The correlation length $\ell$, the mean $\mu$ and the marginal variance
$\bar{c}^{2}$ are the hyperparameters of the GP, which are learned from the
training data set that comprises order-by-order EFT calculations up through
order $k$. Remarkably, it follows from Eq. (8) that the truncation error
defined by Eq. (7) has the distribution

$\delta y_{k}(p)\,|\,\bar{c}^{2},\ell,Q(p),\Lambda\sim\mathcal{GP}[m_{\delta
k}(p),\bar{c}^{2}R_{\delta k}(p,p^{\prime};\ell)]\,,$ (10)

where

$m_{\delta
k}(p)=\frac{y_{\mathrm{ref}}(p)}{\Lambda^{k}}\frac{Q(p)^{k+1}}{\Lambda-Q(p)}\,\mu\,,$
(11)

and

$\displaystyle R_{\delta
k}(p,p^{\prime};\ell)=\frac{y_{\mathrm{ref}}(p)y_{\mathrm{ref}}(p^{\prime})}{\Lambda^{2k}}\frac{[Q(p)Q(p^{\prime})]^{k+1}}{\Lambda^{2}-Q(p)Q(p^{\prime})}\,r(p,p^{\prime};\ell)\,.$
(12)

With point estimates for $\ell$ and $Q(p)$, the normal-inverse-$\chi^{2}$
prior serves as the conjugate prior for the hyperparameters
$(\mu,\bar{c}^{2})$ of the Gaussian processes above, _i.e.,_ their Bayesian
posteriors can be analytically derived and have the same functional forms as
the priors (see Ref. Melendez:2019izc ). The assumptions made above are known
to impose certain limitations bastos_ohagan ; however their validity can be
assessed from several diagnostic metrics on the validation data set (see
below).

### Application to $\sigma_{np}$ and $\sigma_{\gamma d}$—

We now apply the tools described above to the $np\leftrightarrow d\gamma$
cross section. We begin by partitioning the data, _i.e.,_ the order-by-order
$\chi$EFT results depicted in Fig. 1, into training and validation sets. This
data set consists of 32 values of $p$ in 1 MeV increment starting form the
smallest that corresponds to thermal-neutron $p(n,\gamma)d$. These values span
a range of approximately 0 to 1 MeV (2 MeV) in $E$ ($E_{n}$), and encompass
the $4~{}\mathrm{MeV}<p<14~{}\mathrm{MeV}$ interval most relevant for BBN. We
use every fifth point for training and the rest for validation. We will see
further below that our maximum _a posteriori_ (MAP) value of $\ell$ vindicates
this grid and training-validation splitting. Similar results are obtained for
other reasonable choices. As in Ref. Melendez:2019izc , we take $\Lambda=600$
MeV and $Q(p)=(p^{8}+m_{\pi}^{8})/(p^{7}+m_{\pi}^{7})\,,$ which give an
expansion parameter of approximately 0.23 with a very weak dependence on $p$.
We take $y_{\mathrm{ref}}(p)=y_{0}(p)$. The coefficient $c_{5}(p)$ turns out
to be rather unnaturally sized, particularly at smaller values of $p$. We
therefore exclude it from the statistical analysis. This leaves $c_{2,3,4}(p)$
for GP modeling. Finally, we need a value for the _nugget_ $\sigma_{n}^{2}$,
the variance of the Gaussian white noise added to the data to stabilize matrix
inversions during fitting. We find that $\sigma_{n}^{2}>10^{-8}$ is required
to avoid singularities and that $\sigma_{n}^{2}>10^{-3}$ can introduce noises
similar in size to our precision goal of 0.1-0.2% on the predicted cross
sections. We therefore pick $\sigma_{n}^{2}=10^{-5}$ as this value produces
the best performing model on the validation set under the diagnostic criteria
discussed below.

##  III Model calibration, validation and prediction

Figure 2: GP modeling of the $\chi$EFT expansion coefficients and its
diagnostics. (a) The simulators (solid lines) along with the corresponding GP
emulators (dashed lines) and their $2\sigma$ intervals (bands). The training
data are denoted by filled circles; 4 validation points are located uniformly
between each adjacent pair of training points. (b) The Mahalanobis distances
compared to the mean (interior line), $50\%$ (box) and $95\%$ (whiskers)
credible intervals of the reference distribution. (c) The pivoted Cholesky
diagnostics versus the index along with $95\%$ credible intervals (gray
lines). (d) The credible interval diagnostics with $1\sigma$ (dark gray) and
$2\sigma$ (light gray) bands estimated by sampling 1000 GP emulators.

We now present the results of the GP modeling, obtained using the package gsum
Melendez:2019izc . Fig. 2(a) shows the coefficients $c_{n}(p)$ for the
observable $\sigma_{np}$, along with their GP emulators. The MAP value of
$\ell$ is found to be $\ell_{\mathrm{MAP}}=$10.4 MeV. Our choices of the grid
and training-validation splitting places the training points at 5 MeV
intervals. This results in 3 training data points within one correlation
length, which is optimal for GP modeling. Furthermore, our data set spans
approximately $3\ell_{\mathrm{MAP}}$ which, as a rule of thumb, is the range
beyond which the data points become uncorrelated.

Fig. 2(a) provides a beautiful visual indication that the GP model has done an
excellent job of emulating the actual $\chi$EFT calculations. However,
detailed diagnostic checks are needed to quantitatively assess the adequacy of
the model. To this end, we use the diagnostic metrics proposed by Ref.
bastos_ohagan , and originally implemented in EFT error model by Ref.
Melendez:2019izc .

Fig. 2(b) shows the squared Mahalanobis distances, defined as

$\mathrm{D}^{2}_{\mathrm{MD}}=(\mathbf{f}-\mathbf{m})^{\mathrm{T}}K^{-1}(\mathbf{f}-\mathbf{m})\,,$
(13)

where we have used the notation $\mathbf{f}$ for the vector of validation data
points, $\mathbf{m}$ for the vector of means of the emulator at these points
and $K$ for its covariance matrix. $\mathrm{D}^{2}_{\mathrm{MD}}$ is a
generalization of the sum of squared residuals to the case of data points that
are correlated across the independent variable. Values much larger than its
reference distribution, a $\chi^{2}$, indicate conflict between the emulator
and the simulator. Values much smaller than the reference, as we see in the
case of $c_{2}$ to some extent, indicate that, given statistical fluctuations,
there is an unusually close match between the emulator and the simulator.

We now look at a more informative metric, the pivoted Cholesky diagnostic
$\mathrm{D}_{\mathrm{PC}}$. For a “correct” emulator, it returns draws from a
standard Gaussian at every index, which represents the validation points
arranged such that the first element is the one with the largest predictive
variance, the second element is the one with the largest variance conditioned
on the first element, and so on. A group of unusually large or small
$\mathrm{D}_{\mathrm{PC}}$ values across all indices indicates a misestimated
variance whereas a group of unusually large or small
$\mathrm{D}_{\mathrm{PC}}$ values in the latter part of the sequence indicates
an inappropriate correlation structure. Overall, Fig. 2(c) shows that the
points are distributed as expected, _e.g.,_ 4 out of the 72 points lie outside
the $-2<\mathrm{D}_{\mathrm{PC}}<2$ range. We also notice that there is a
slight indication that the variance on $c_{2}$ ($c_{3}$) might have been
somewhat overestimated (underestimated) by observing the spread of the
corresponding data points.

The credible interval diagnostic involves constructing uncertainty bands at
each order and checking whether it actually encompasses the correction that
enters at the next order. The claimed $(1-\alpha)100\%$ credible intervals are
then plotted against the percentage of validation data points found within the
interval—emulators that output credible intervals containing too few data
points compared to the reference distribution are overconfident and those that
contain too many are underconfident. For uncorrelated data points, the
reference distribution is a binomial. For correlated data points, the
reference distribution is numerically estimated by sampling a large number of
emulators from the underlying process. Fig. 2(d) shows that the model is
performing as expected and that it is important to account for correlations
while assigning truncation errors.

Now that we have demonstrated that the coefficients $c_{2,3,4}$ can be
appropriately described by a GP, we can now use Eq. (10) to compute the
truncation errors. We list the $np\rightarrow d\gamma$ cross section values at
order $k=4$ along with their $1\sigma$ truncation errors at several energies
in Table 2. We note that these errors are different from, and vary much more
smoothly with $E$, than naive estimates PhysRevLett.115.122301  based on
multiplying the largest order-to-order shift with the appropriate power of the
expansion parameter in a pointwise manner at each value of $E$. In Fig. 3(a),
we plot $\sigma_{np}v_{n}$, for which we showed order-by-order results earlier
in Fig. 1(a), versus $E_{n}$. This quantity is proportional to the reaction
rate in BBN. The bands represent 2$\sigma$ truncation errors, _i.e.,_ 95%
Bayesian credible intervals for $\sigma_{np}v_{n}$. In Fig. 3(b), we show
these bands for $\sigma_{\gamma d}$ and compare them with photodissociation
data shown earlier in Fig. 1(b). Reassuringly, the experimental data as well
as the highest-order theoretical calculation, $y_{5}$, are compatible with the
truncation error estimates.

$E$ [MeV] |  $\sigma_{np}$ [mb]
---|---
$1.262500\times 10^{-08}$ |  $321.009\pm 0.71496$
$9.607513\times 10^{-03}$ |  $0.32739\pm 0.00073$
$3.838601\times 10^{-02}$ |  $0.12762\pm 0.00029$
$8.633551\times 10^{-02}$ |  $0.06853\pm 0.00015$
$1.534560\times 10^{-01}$ |  $0.04658\pm 0.00010$
$2.397475\times 10^{-01}$ |  $0.03792\pm 0.00008$
$3.452100\times 10^{-01}$ |  $0.03464\pm 0.00007$
$4.698435\times 10^{-01}$ |  $0.03368\pm 0.00007$
$6.136480\times 10^{-01}$ |  $0.03373\pm 0.00007$
$7.766235\times 10^{-01}$ |  $0.03414\pm 0.00007$
$9.587699\times 10^{-01}$ |  $0.03461\pm 0.00007$
Table 2: $\chi$EFT predictions at order $k=4$ for the $np\rightarrow d\gamma$
cross section at $np$ relative energy $E$, along with their $1\sigma$ errors
from the truncation of the $\chi$EFT potential. Figure 3: The $2\sigma$
truncation error bands on the $\chi$EFT predictions $y_{k}$ at $k=2,3,4$ along
with the prediction $y_{5}$ and data from Fig. 1. $(a)$ The product of
$p(n,\gamma)d$ cross section and the neutron speed versus the energy of the
neutron. $(b)$ The deuteron photodissociation cross section as a function of
the photon energy in the rest frame of the deuteron.

The uncertainties quoted in Table 2 and Fig. 3, which amount to about 0.2$\%$,
only include truncation error in the $\chi$EFT potential. The full theory
uncertainty also comprises statistical error from fitting the LECs to
experimental data. The framework we have adopted allows one to also
incorporate fitting uncertainties on the LECs that appear up to the EFT order
we have calculated Melendez:2019izc . Another missing source of uncertainty,
which is important in the $M1$-dominated regime, is the truncation of the
current. Inclusion of the $M1$ operator at order $(Q/\Lambda)^{1}$ is crucial
for obtaining agreement with the experimental value for threshold capture
given in Table 1 piarulli ; Phillips:2016mov , although it introduces several
new LECs and thus poses significant challenges for rigorous uncertainty
analysis. The fitting strategy that uses minimal assumptions about the short-
distance behavior of the current operator, among several explored by Ref.
piarulli , is the one that constrains the LECs $d_{1,2}^{V}$ simultaneously to
$\sigma_{np}(E=1.2625\times 10^{-08}~{}\mathrm{MeV})$ and the isovector
combination of the $A=3$ magnetic moments. However, it was found that this
yields unnatural values for $d_{1,2}^{V}$. This fine-tuning can be mitigated
by including the theory uncertainty we calculated in this paper, as well as
the experimental error, in the fit. Such a strategy for performing parameter
estimation with $\chi$EFT truncation error included as a guard against
overfitting was recently successfully pursued by Ref. Wesolowski:2021cni  in
the context of constraining 3N interactions from properties of $A=3,4$ nuclei.
A calculation of $np\leftrightarrow d\gamma$ along these lines is a subject
for future work.

##  IV Summary and outlook

We performed the first $\chi$EFT calculations of the energy-dependent
$np\leftrightarrow d\gamma$ cross section at low energies, including the range
relevant to BBN, and the first Bayesian analysis of $\chi$EFT truncation error
to a nuclear reaction cross section. Working with fixed one- and two-body
electromagnetic current operators, we studied the convergence of this
observable in the EFT expansion of the potential. By harnessing recent
progress in Bayesian analysis of EFT uncertainty, we were able to provide
statistical estimates, amounting to 0.2$\%$, for the theory uncertainty that
stems from truncation of the $\chi$EFT potential.

At the $\chi$EFT order up to which we work, our calculations are pure
predictions as no new LECs enter that need to be fixed. Inclusion of the next
order in the current operator adds new parameters that are not well
constrained at present and will most likely require fitting to electromagnetic
observables in $A=2,3$ systems. To make predictions for this cross section
with subpercentage-level precision even at threshold, we will therefore need
further investigation into the nature of the current operator at short
distances and calculation of the electromagnetic observables for $A=3$ nuclei
with NN and 3N interactions truncated consistently. Finally, uncertainty
analysis in theoretical calculations of Deuterium burning processes will be an
important development because these strongly affect the primordial Deuterium
abundance and there is currently some discrepancy between theory and
experiments, most notably for 2H$(p,\gamma){}^{3}$He Pisanti:2020efz .

## Acknowledgements

We are grateful to Jordan Melendez, Richard J Furnstahl and Daniel R Phillips
for fruitful discussions. We would also like to thank Daniel R Phillips for a
critical reading of the manuscript, Evgeny Epelbaum for providing us with
computer programs for chiral potentials used in this work and the BUQEYE
collaboration for making the library gsum publicly accessible. This work was
supported by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation), through the Cluster of Excellence [Precision Physics, Fundamental
Interactions, and Structure of Matter] (PRISMA+ EXC 2118/1) within the German
Excellence Strategy (Project ID 39083149).

## References

  * (1) E. Epelbaum, H.-W. Hammer, U.-G. Meißner, Modern Theory of Nuclear Forces, Rev. Mod. Phys. 81 (2009) 1773.  doi:10.1103/RevModPhys.81.1773.
  * (2) M. R. Schindler, D. R. Phillips, Bayesian Methods for Parameter Estimation in Effective Field Theories, Annals Phys. 324 (2009) 682, [Erratum: Annals Phys. 324 (2009) 2051].  doi:10.1016/j.aop.2008.09.003.
  * (3) R. J. Furnstahl, D. R. Phillips, S. Wesolowski, A recipe for EFT uncertainty quantification in nuclear physics, J. Phys. G 42 (2015) 034028.  doi:10.1088/0954-3899/42/3/034028.
  * (4) R. J. Furnstahl, N. Klco, D. R. Phillips, S. Wesolowski, Quantifying truncation errors in effective field theory, Phys. Rev. C 92 (2015) 024005.  doi:10.1103/PhysRevC.92.024005.
  * (5) S. Wesolowski, N. Klco, R. J. Furnstahl, D. R. Phillips, A. Thapaliya, Bayesian parameter estimation for effective field theories, J. Phys. G 43 (2016) 074001.  doi:10.1088/0954-3899/43/7/074001.
  * (6) X. Zhang, K. M. Nollett, D. R. Phillips, Halo effective field theory constrains the solar 7Be + p $rightarrow$ 8B + $\gamma$ rate, Phys. Lett. B 751 (2015) 535.  doi:10.1016/j.physletb.2015.11.005.
  * (7) J. A. Melendez, S. Wesolowski, R. J. Furnstahl, Bayesian truncation errors in chiral effective field theory: nucleon-nucleon observables, Phys. Rev. C 96 (2017) 024003.  doi:10.1103/PhysRevC.96.024003.
  * (8) S. Wesolowski, R. J. Furnstahl, J. A. Melendez, D. R. Phillips, Exploring Bayesian parameter estimation for chiral effective field theory using nucleon–nucleon phase shifts, J. Phys. G 46 (2019) 045102.  doi:10.1088/1361-6471/aaf5fc.
  * (9) A. E. Lovell, F. M. Nunes, Constraining Transfer Cross Sections Using Bayes’ Theorem, Phys. Rev. C 97 (2018) 064612.  doi:10.1103/PhysRevC.97.064612.
  * (10) L. Neufcourt, Y. Cao, W. Nazarewicz, F. Viens, Bayesian approach to model-based extrapolation of nuclear observables, Phys. Rev. C 98 (2018) 034318\.  doi:10.1103/PhysRevC.98.034318.
  * (11) S. Zhou, P. Giulani, J. Piekarewicz, A. Bhattacharya, D. Pati, Reexamining the proton-radius problem using constrained Gaussian processes, Phys. Rev. C 99 (2019) 055202.  doi:10.1103/PhysRevC.99.055202.
  * (12) G. B. King, A. E. Lovell, L. Neufcourt, F. M. Nunes, Direct comparison between Bayesian and frequentist uncertainty quantification for nuclear reactions, Phys. Rev. Lett. 122 (2019) 232502.  doi:10.1103/PhysRevLett.122.232502.
  * (13) J. A. Melendez, R. J. Furnstahl, D. R. Phillips, M. T. Pratola, S. Wesolowski, Quantifying Correlated Truncation Errors in Effective Field Theory, Phys. Rev. C 100 (2019) 044001.  doi:10.1103/PhysRevC.100.044001.
  * (14) P. Premarathna, G. Rupak, Bayesian analysis of capture reactions ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$, Eur. Phys. J. A 56 (2020) 166.  doi:10.1140/epja/s10050-020-00113-z.
  * (15) X. Zhang, K. M. Nollett, D. R. Phillips, $S$-factor and scattering-parameter extractions from ${}^{3}\mathrm{He}+{}^{4}\mathrm{He}\rightarrow{}^{7}\mathrm{Be}+\gamma$, J. Phys. G 47 (2020) 054002.  doi:10.1088/1361-6471/ab6a71.
  * (16) E. Epelbaum, et al., Towards high-order calculations of three-nucleon scattering in chiral effective field theory, Eur. Phys. J. A 56 (2020) 92.  doi:10.1140/epja/s10050-020-00102-2.
  * (17) D. R. Phillips, et al., Get on the BAND Wagon: A Bayesian Framework for Quantifying Model Uncertainties in Nuclear Dynamics, J. Phys. G 48 (2021) 072001\.  doi:10.1088/1361-6471/abf1df.
  * (18) Y. Volkotrub, J. Golak, R. Skibiński, K. Topolnicki, H. Witała, E. Epelbaum, H. Krebs, P. Reinert, Uncertainty of three-nucleon continuum observables arising from uncertainties of two-nucleon potential parameters, J. Phys. G 47 (2020) 104001.  doi:10.1088/1361-6471/ab9b07.
  * (19) P. Reinert, H. Krebs, E. Epelbaum, Precision determination of pion-nucleon coupling constants using effective field theory, Phys. Rev. Lett. 126 (2021) 092501\.  doi:10.1103/PhysRevLett.126.092501.
  * (20) S. Wesolowski, I. Svensson, A. Ekström, C. Forssén, R. J. Furnstahl, J. A. Melendez, D. R. Phillips, Rigorous constraints on three-nucleon forces in chiral effective field theory from fast and accurate calculations of few-body observables, Phys. Rev. C 104 (2021) 064001\.  doi:10.1103/PhysRevC.104.064001.
  * (21) J. Rijneveen, N. Rijneveen, H. Krebs, A. M. Gasparyan, E. Epelbaum, Radiative pion photoproduction in covariant chiral perturbation theory, Phys. Rev. C 103 (2021) 045203.  doi:10.1103/PhysRevC.103.045203.
  * (22) P. Maris, et al., Light nuclei with semilocal momentum-space regularized chiral interactions up to third order, Phys. Rev. C 103 (2021) 054001.  doi:10.1103/PhysRevC.103.054001.
  * (23) C. Wellenhofer, C. Drischler, A. Schwenk, Effective field theory for dilute Fermi systems at fourth order, Phys. Rev. C 104 (2021) 014003.  doi:10.1103/PhysRevC.104.014003.
  * (24) M. A. Connell, I. Billig, D. R. Phillips, Does Bayesian model averaging improve polynomial extrapolations? Two toy problems as tests, J. Phys. G 48 (2021) 104001.  doi:10.1088/1361-6471/ac215a.
  * (25) N. Rijneveen, A. M. Gasparyan, H. Krebs, E. Epelbaum, Pion photoproduction in chiral perturbation theory with explicit treatment of the $\Delta(1232)$ resonance arXiv:2108.01619.
  * (26) V. Mossa, et al., The baryon density of the Universe from an improved rate of deuterium burning, Nature 587 (2020) 210.  doi:10.1038/s41586-020-2878-4.
  * (27) P. D. Serpico, S. Esposito, F. Iocco, G. Mangano, G. Miele, O. Pisanti, Nuclear reaction network for primordial nucleosynthesis: a detailed analysis of rates, uncertainties and light nuclei yields, Journal of Cosmology and Astroparticle Physics 12 (2004) 010.  doi:10.1088/1475-7516/2004/12/010.
  * (28) O. Pisanti, A. Cirillo, S. Esposito, F. Iocco, G. Mangano, G. Miele, P. D. Serpico, PArthENoPE: Public Algorithm Evaluating the Nucleosynthesis of Primordial Elements, Comput. Phys. Commun. 178 (2008) 956.  doi:10.1016/j.cpc.2008.02.015.
  * (29) A. Coc, C. Pitrou, J.-P. Uzan, E. Vangioni, A Public Code for Precision Big Bang Nucleosynthesis with Improved Helium-4 Predictions, Springer Proc. Phys. 219 (2019) 327.  doi:10.1007/978-3-030-13876-9_57.
  * (30) C. Pitrou, A. Coc, J.-P. Uzan, E. Vangioni, Precision big bang nucleosynthesis with improved Helium-4 predictions, Phys. Rept. 754 (2018) 1.  doi:10.1016/j.physrep.2018.04.005.
  * (31) N. Aghanim, et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [Erratum: Astron.Astrophys. 652, C4 (2021)].  doi:10.1051/0004-6361/201833910.
  * (32) R. J. Cooke, M. Pettini, C. C. Steidel, One Percent Determination of the Primordial Deuterium Abundance, Astrophys. J. 855 (2018) 102.  doi:10.3847/1538-4357/aaab53.
  * (33) H. Arenhovel, M. Sanzone, Photodisintegration of the deuteron: A Review of theory and experiment, Few Body Syst. Suppl. 3 (1991) 1.
  * (34) J. Carlson, R. Schiavilla, Structure and Dynamics of Few Nucleon Systems, Rev. Mod. Phys. 70 (1998) 743.  doi:10.1103/RevModPhys.70.743.
  * (35) L. E. Marcucci, K. M. Nollett, R. Schiavilla, R. B. Wiringa, Modern theories of low-energy astrophysical reactions, Nucl. Phys. A 777 (2006) 111.  doi:10.1016/j.nuclphysa.2004.09.008.
  * (36) J.-W. Chen, M. J. Savage, $np\rightarrow d\gamma$ for big bang nucleosynthesis, Phys. Rev. C 60 (1999) 065205.  doi:10.1103/PhysRevC.60.065205.
  * (37) G. Rupak, Precision calculation of $np\rightarrow d\gamma$ cross section for big-bang nucleosynthesis, Nuclear Physics A 678 (2000) 405.  doi:https://doi.org/10.1016/S0375-9474(00)00323-7.
  * (38) S. Ando, R. H. Cyburt, S. W. Hong, C. H. Hyun, Radiative neutron capture on a proton at big-bang nucleosynthesis energies, Phys. Rev. C 74 (2006) 025809.  doi:10.1103/PhysRevC.74.025809.
  * (39) D. Cokinos, E. Melkonian, Measurement of the 2200 m/sec neutron-proton capture cross section, Phys. Rev. C 15 (1977) 1636.  doi:10.1103/PhysRevC.15.1636.
  * (40) S. R. Beane, E. Chang, W. Detmold, K. Orginos, A. Parreño, M. J. Savage, B. C. Tiburzi, Ab initio Calculation of the $np\rightarrow{}d\gamma$ Radiative Capture Process, Phys. Rev. Lett. 115 (2015) 132001.  doi:10.1103/PhysRevLett.115.132001.
  * (41) M. Piarulli, L. Girlanda, L. E. Marcucci, S. Pastore, R. Schiavilla, M. Viviani, Electromagnetic structure of $A=2$ and 3 nuclei in chiral effective field theory, Phys. Rev. C 87 (2013) 014006.  doi:10.1103/PhysRevC.87.014006.
  * (42) Y.-H. Song, R. Lazauskas, T.-S. Park, Up to $\mathrm{N}^{3}\mathrm{LO}$ heavy-baryon chiral perturbation theory calculation for the $M1$ properties of three-nucleon systems, Phys. Rev. C 79 (2009) 064002.  doi:10.1103/PhysRevC.79.064002.
  * (43) S. Burles, K. M. Nollett, J. W. Truran, M. S. Turner, Sharpening the predictions of big-bang nucleosynthesis, Phys. Rev. Lett. 82 (1999) 4176.  doi:10.1103/PhysRevLett.82.4176.
  * (44) C. Drischler, J. A. Melendez, R. J. Furnstahl, D. R. Phillips, Quantifying uncertainties and correlations in the nuclear-matter equation of state, Phys. Rev. C 102 (2020) 054315.  doi:10.1103/PhysRevC.102.054315.
  * (45) C. Drischler, R. J. Furnstahl, J. A. Melendez, D. R. Phillips, How well do we know the neutron-matter equation of state at the densities inside neutron stars? a bayesian approach with correlated uncertainties, Phys. Rev. Lett. 125 (2020) 202702.  doi:10.1103/PhysRevLett.125.202702.
  * (46) P. Reinert, H. Krebs, E. Epelbaum, Semilocal momentum-space regularized chiral two-nucleon potentials up to fifth order, Eur. Phys. J. A 54 (2018) 86.  doi:10.1140/epja/i2018-12516-4.
  * (47) T.-S. Park, D.-P. Min, M. Rho, Chiral Lagrangian approach to exchange vector currents in nuclei, Nucl. Phys. A596 (1996) 515–552.  arXiv:nucl-th/9505017, doi:10.1016/0375-9474(95)00406-8.
  * (48) S. Pastore, R. Schiavilla, J. L. Goity, Electromagnetic two-body currents of one- and two-pion range, Phys. Rev. C78 (2008) 064002.  arXiv:0810.1941, doi:10.1103/PhysRevC.78.064002.
  * (49) S. Kolling, E. Epelbaum, H. Krebs, U.-G. Meißner, Two-pion exchange electromagnetic current in chiral effective field theory using the method of unitary transformation, Phys. Rev. C80 (2009) 045502.  doi:10.1103/PhysRevC.80.045502.
  * (50) B. Acharya, S. Bacca, Neutrino-deuteron scattering: Uncertainty quantification and new ${L}_{1,A}$ constraints, Phys. Rev. C 101 (2020) 015505.  doi:10.1103/PhysRevC.101.015505.
  * (51) B. Acharya, V. Lensky, S. Bacca, M. Gorchtein, M. Vanderhaeghen, Dispersive evaluation of the lamb shift in muonic deuterium from chiral effective field theory, Phys. Rev. C 103 (2021) 024001.  doi:10.1103/PhysRevC.103.024001.
  * (52) T. Suzuki, Y. Nagai, T. Shima, T. Kikuchi, H. Sato, T. Kii, M. Igashira, First measurement of a $p(n,\gamma)d$ reaction cross section between 10 and 80 kev, The Astrophysical Journal 439 (1995) L59.
  * (53) Y. Nagai, T. S. Suzuki, T. Kikuchi, T. Shima, T. Kii, H. Sato, M. Igashira, Measurement of ${}^{1}\mathrm{H}(n,\gamma{)}^{2}\mathrm{H}$ reaction cross section at a comparable $M1/E1$ strength, Phys. Rev. C 56 (1997) 3173.  doi:10.1103/PhysRevC.56.3173.
  * (54) K. Y. Hara, H. Utsunomiya, S. Goko, H. Akimune, T. Yamagata, M. Ohta, H. Toyokawa, K. Kudo, A. Uritani, Y. Shibata, Y.-W. Lui, H. Ohgaki, Photodisintegration of deuterium and big bang nucleosynthesis, Phys. Rev. D 68 (2003) 072001.  doi:10.1103/PhysRevD.68.072001.
  * (55) R. Moreh, T. J. Kennett, W. V. Prestwich, ${}^{2}\mathrm{H}$($\gamma$,n) absolute cross section at 2754 kev, Phys. Rev. C 39 (1989) 1247.  doi:10.1103/PhysRevC.39.1247.
  * (56) W. Tornow, N. Czakon, C. Howell, A. Hutcheson, J. Kelley, V. Litvinenko, S. Mikhailov, I. Pinayev, G. Weisel, H. Witała, Low-energy photodisintegration of the deuteron and big-bang nucleosynthesis, Physics Letters B 574 (2003) 8.  doi:https://doi.org/10.1016/j.physletb.2003.08.078.
  * (57) L. S. Bastos, A. O’Hagan, Diagnostics for gaussian process emulators, Technometrics 51 (2009) 425.  doi:10.1198/TECH.2009.08019.
  * (58) E. Epelbaum, H. Krebs, U.-G. Meißner, Precision Nucleon-Nucleon Potential at Fifth Order in the Chiral Expansion, Phys. Rev. Lett. 115 (2015) 122301.  doi:10.1103/PhysRevLett.115.122301.
  * (59) D. R. Phillips, Electromagnetic Structure of Two- and Three-Nucleon Systems: An Effective Field Theory Description, Ann. Rev. Nucl. Part. Sci. 66 (2016) 421–447.  doi:10.1146/annurev-nucl-102014-022321.
  * (60) O. Pisanti, G. Mangano, G. Miele, P. Mazzella, Primordial Deuterium after LUNA: concordances and error budget, JCAP 04 (2021) 020.  doi:10.1088/1475-7516/2021/04/020.
  * (61) K. Borah, R. J. Hill, G. Lee, O. Tomalak, Parametrization and applications of the low-$Q^{2}$ nucleon vector form factors, Phys. Rev. D 102 (2020) 074012.  doi:10.1103/PhysRevD.102.074012.
  * (62) To include nucleon-structure effects, we adopt a recent model-independent parametrization of the nucleon form factors Borah:2020gte  in lieu of summing the corresponding diagrams in $\chi$EFT.