aainstitutetext: Fermi National Accelerator Laboratory, Batavia,
USAbbinstitutetext: Department of Physics & Astronomy, Rutgers, The State
University of New Jersey,
Piscataway, USAccinstitutetext: Department of Physics, Carnegie Mellon
University, Pittsburgh, USA

# Model selection and signal extraction using Gaussian Process regression

Abhijith Gandrakota, b Amitabh Lath b Alexandre V. Morozov c Sindhu Murthy
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a novel computational approach for extracting weak signals, whose
exact location and width may be unknown, from complex background distributions
with an arbitrary functional form. We focus on datasets that can be naturally
presented as binned integer counts, demonstrating our approach on the CERN
open dataset from the ATLAS collaboration at the Large Hadron Collider, which
contains the Higgs boson signature. Our approach is based on Gaussian Process
(GP) regression - a powerful and flexible machine learning technique that
allowed us to model the background without specifying its functional form
explicitly, and to separate the background and signal contributions in a
robust and reproducible manner. Unlike functional fits, our GP-regression-
based approach does not need to be constantly updated as more data becomes
available. We discuss how to select the GP kernel type, considering trade-offs
between kernel complexity and its ability to capture the features of the
background distribution. We show that our GP framework can be used to detect
the Higgs boson resonance in the data with more statistical significance than
a polynomial fit specifically tailored to the dataset. Finally, we use Markov
Chain Monte Carlo (MCMC) sampling to confirm the statistical significance of
the extracted Higgs signature.

##  1 Introduction

Analyzing data from physical experiments or observations often involves
fitting computational models in order to extract a signal in the presence of
both background effects and random noise. For example, such a setting appears
naturally in the analysis of X-ray diffraction patterns from crystalline
samples that contain contributions both from distinctive Bragg peaks and
diffuse background scattering Sivia1994 ; David2001 , inference of transiting
exoplanet parameters in astronomy exoplanet_astro ; gp_astro_general , and the
discovery of the Higgs boson and search for the new physics at the Large
Hadron Collider (LHC) at CERN Aad_2012 . The data from LHC and similar
experiments usually comes in the form of binned integer counts Zyla:2020zbs .
Traditionally, modeling such data is performed under the assumption of the
Poisson distribution, by employing a parametric fit Zyla:2020zbs . The fitted
models are subsequently used to estimate the background contributions and
extract the signal of interest. However, the choice of parametric functions is
often ad-hoc and the degree of model complexity requires a delicate balancing
act between overfitting and underfitting the data bishop2006pattern ;
Mehta2019 ; Rocks2020 .

For data analyses of the type performed on the LHC bin counts, the optimal
complexity of the model is usually evaluated by performing Wilk’s tests
Wilks1938 . Most analyses that employ this technique test it on a fraction of
the full dataset, usually 10% or less. This process is called “blinding” and
is meant to reduce biases in the analysis. However, the model selection
process must be repeated periodically as more data becomes available, and
often the functional form employed in the model has to be updated as well
Aad_2016 . Using non-parametric methods such as Gaussian Process regression is
an effective way of alleviating these concerns pmlr-v5-titsias09a ;
bishop2006pattern ; gpml ; Kersting2007 ; ohagan .

Gaussian Process (GP) regression is a well-established machine learning
technique gpml ; bishop2006pattern  commonly used in various fields such as
astrophysics, gravitational wave detection, and high energy physics Iyer_2019
; Moore_2016 ; golchi2015 ; frate2017modeling . In particular, in high energy
physics GP regression was used to model the smooth continuum background from
quantum chromodynamics (QCD) in searches for dijet resonances in LHC data
frate2017modeling . The authors argued that using GP for background estimation
was more robust with respect to increasing luminosity compared to parametric
fitting methods. GP regression’s advantages over more conventional methods,
which employ a linear expansion over a fixed set of basis functions such as
polynomials or Gaussians, are due to its non-parametric flexibility and a
principled Bayesian framework. Instead of explicit basis functions, GP
regression is defined in terms of kernel functions which specify the degree of
correlations between two points in the dataset. The GP approach allows us to
perform inference using a much broader class of functions, including those
which would otherwise require an infinite basis set bishop2006pattern . GP
regression is robust with respect to the size of the dataset frate2017modeling
.

Nevertheless, the flexibility of GP regression can be a double-edged sword. In
GP regression, kernel functions typically depend on several hyperparameters
that are varied to fit the data, typically through non-Bayesian techniques
such as maximizing the marginal likelihood of the observed data gpml ;
bishop2006pattern . The hyperparameters describing the kernel function control
the flexibility of the resulting model, while the type of the kernel function
determines the success in capturing certain features in the data, such as
periodic oscillations and long-term trends gpml . Thus, the universality and
the power of the GP approach may come at the cost of overfitting with respect
to both kernel type and kernel hyperparameter choices. Therefore, a method is
required that can constrain the flexibility of the GP regression in a
controlled manner.

Previous work in this area has focused on “kernel learning” to address the
issues of flexibility and robustness, with several techniques proposed that
aim at constructing composite kernels for Support Vector Machines
Diosan2007EvolvingKF , Relevance Vector Machines Bing2010 , and GP regression
duvenaud2013structure ; Wilson2013  using a library of base kernels.
Semiparametric regression attempts to combine interpretability of parametric
models with flexibility of non-parametric models by combining the two
approaches in a single framework Ruppert2003 . However, none of the above
approaches focus specifically on integer count data or on the processes that
are naturally viewed as localized signals superimposed on the smooth
background. A previous application of GP regression to LHC data
frate2017modeling  employed both standard and custom-built kernels motivated
by physical considerations. In contrast, in this work we have developed both a
model selection procedure suitable for GP regression and an approach for
estimating the statistical significance of the extracted signal.

New signals in physical observations of particle resonances in LHC data often
appear as localized features (“bumps”) superimposed on a smooth background.
Accurate modeling of the background spectrum is therefore essential to both
extracting the signal and assessing its statistical significance. In this
paper, we present a rigorous approach to model selection in GP regression
applied to binned integer data, which we expect to be a superposition of a
localized signal and a smooth background of unknown functional form. We
exploit the flexibility of the GP regression by determining the kernel
hyperparameters through the fit to background-only data, with the signal
window masked out. These parameters are subsequently used to extrapolate the
background contribution across the signal window, enabling us to separate the
background from the signal contribution. We describe procedures for kernel
type selection based on both Bayesian and Akaike information criteria. We also
propose a method for estimating the statistical significance of the signal by
performing a hypothesis test with data devoid of signal as the null hypothesis
and data containing signal as the alternate hypothesis. While similar in
spirit to standard hypothesis-testing approaches, our significance test takes
into account both the uncertainties inherent in the Bayesian nature of GP
regression and the sampling noise related to generating integer bin counts
from GP-predicted, real-valued Poisson rates.

In this work, we illustrate our procedure by detecting for the Higgs boson
resonance in the open data collected by the ATLAS experiment at the LHC ATL-
OREACH-PUB-2020-001 . We show that using GP regression leads to extracting the
Higgs boson signature at a higher level of statistical significance compared
to parametric fits. Our computational pipeline can be applied for background
estimation and signal detection in any dataset where a localized signal is
obscured by background processes.

##  2 Model selection and signal extraction procedure

###  2.1 Gaussian process regression

In regression, the observed data $y=(y_{1}\dots y_{N})$ (in our case, integer
counts in $N$ bins) is modeled by $z=(z(x_{1})\dots z(x_{N}))$:

$y_{k}=z(x_{k})+\epsilon_{k},$ (1)

where $k=1\dots N$ ($N$ is the total number of datapoints), $X=(x_{1}\dots
x_{N})$ is a vector of input variables (in our case, centers of the bins with
integer counts), and $\epsilon_{k}$ is a random noise variable independently
sampled from a Gaussian distribution ${\cal N}(\epsilon|0,\sigma^{2}_{i})$ for
each data point, where $\sigma^{2}_{i}$ is the noise variance in bin $i$.

In the GP framework, the model $z(x)$ is not represented as an explicit linear
expansion over a set of pre-determined basis functions. Instead, we directly
consider the marginal likelihood $p(y|X)$ integrated over all possible models
Goldberg1998 ; bishop2006pattern :

$p(y|X)=\int dy~{}p(y|z)p(z|X)={\cal N}(y|m,C),$ (2)

where $m=(m(x_{1})\dots m(x_{N}))$ is a vector of values of the mean function
$m(x)$ for all datapoints, and the covariance matrix $C=K+\Sigma$, where $K$
is the Gram matrix and $\Sigma$ is a diagonal $N\times N$ matrix with
$\Sigma_{ii}=\sigma^{2}_{i}$. The elements of the Gram matrix are values of
the kernel function $k(x,x^{\prime})$ evaluated for all pairs of input
variables: $K_{ij}=k(x_{i},x_{j})$. Note that, consistent with Eq. (1),
$p(y|z)={\cal N}(y|z,\Sigma)$, while $p(z|X)={\cal N}(z|m,K)$ from the
definition of the Gaussian Process.

Thus, GP regression is defined by the mean function $m(x)$ and the kernel
function $k(x,x^{\prime})$ gpml ; bishop2006pattern , which determines the
degree of correlation between any two datapoints. In general, the kernel
function depends on a set of $n$ hyperparameters
$\theta=(\theta_{1},\theta_{2},\dots,\theta_{n})$, whose number and meaning
depend on the kernel type. The hyperparameters of a given kernel are usually
optimized by maximizing the marginal likelihood in Eq. (2), a non-Bayesian
procedure gpml ; bishop2006pattern . Ordinarily, kernel hyperparameters
include $\sigma^{2}_{i}$, which represent the amount of experimental noise in
each bin. However, since this would introduce too many hyperparameters and
make their optimization difficult or impossible, we estimate $\sigma^{2}_{i}$
directly from the data using Garwood intervals which allow us to extract two-
sided confidence intervals from the number of events in each bin under the
assumption that the events are Poisson-distributed Garwood1936 ; Demortier2013
. Thus, $\sigma^{2}_{i}$ are estimated independently and are not treated as
hyperparameters in our approach.

The strength of the GP approach stems from the fact that the joint marginal
probability of observing a set of datapoints is Gaussian (Eq. (2)). Moreover,
the predictive probability $p(\tilde{y}_{i}|y)$, the conditional probability
distribution of observing a real-valued “count” $\tilde{y}_{i}$ in bin $i$
given a dataset with $N$ previous observations, is also Gaussian, with the
mean $f(x_{i})$ and the variance $V(x_{i})$ given by:

$\displaystyle f(x_{i})$ $\displaystyle=$ $\displaystyle
m(x_{i})+\tilde{k}^{T}C^{-1}\tilde{y}_{i},$ (3) $\displaystyle V(x_{i})$
$\displaystyle=$ $\displaystyle\alpha-\tilde{k}^{T}C^{-1}\tilde{k},$

where $\tilde{k}=(k(x_{1},x_{i}),\dots,k(x_{N},x_{i}))$ and
$\alpha=k(x_{i},x_{i})+\sigma^{2}_{i}$.

In this work, we consider two types of GP regression: with
$m(x_{i})=0,~{}\forall i$ for modeling the background-only distribution, and
with the Gaussian mean function for modeling signal+background datasets, where
the signal component is represented by:

$m(x_{i})=\frac{A}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\left(x_{i}-\mu\right)^{2}}{2\sigma^{2}}\right).$
(4)

Here, $A$ defines the signal strength, while $\mu$ and $\sigma$ represent
signal mean and width, respectively. The rounded value of $A$ can be
interpreted as the total number of signal events. Note that when the Gaussian
mean function is introduced, the set of model hyperparameters $\theta$ needs
to be augmented with $\\{A,\mu,\sigma\\}$.

###  2.2 Model selection

A key issue in GP regression is the choice of a kernel and, given a kernel,
derivation of the optimal set of hyperparameters $\hat{\theta}$ for it.
Typically, the optimal set of hyperparameters is obtained by maximizing the
marginal log-likelihood $\log p(y|X,\theta,K_{i})$ gpml ; bishop2006pattern ,
where $p(y|X,\theta,K_{i})$ is given by Eq. (2) and its dependence on the set
of hyperparameters $\theta$ and the kernel type $K_{i}$ are made explicit for
clarity. Since this step is non-Bayesian, a question of kernel selection
arises which would take into account both kernel complexity (i.e., the amount
of signal smoothing provided by a given kernel) and the number of kernel
hyperparameters. A standard way for carrying out model comparison is based on
the Bayesian Information Criterion (BIC) for the marginal log-likelihood lhr ;
bishop2006pattern :

$-\log p(y|X,K_{i})\simeq-\log p(y|X,\hat{\theta},K_{i})+\frac{n}{2}\log
N\equiv\text{BIC}^{\text{naive}},$ (5)

where $p(y|X,K_{i})$ is the model evidence (likelihood marginalized over
hyperparameters), $n$ is the number of model parameters, and $N$ is the number
of datapoints. Note that the second term on the right-hand side penalizes
model complexity, such that lower BIC scores are more preferable. The
derivation of BIC relies on a number of approximations whose validity depends
on the details of the system under consideration. Specifically, the derivation
employs the Laplace approximation to estimate the integral over the
hyperparameters and assumes that $N$ is so large (or the Gaussian prior
distribution over the hyperparameters is so broad) that the effects of the
hyperparameter priors are negligible, resulting in:

$-\log p(y|X,K_{i})\simeq-\log
p(y|X,\hat{\theta},K_{i})+\frac{1}{2}\log|H|\equiv\text{BIC},$ (6)

where $H=-\vec{\nabla}_{\theta}\vec{\nabla}_{\theta}\log
p(y|X,\theta,K_{i})|_{\hat{\theta}}$ is the Hessian in the model
hyperparameter space evaluated at the hyperparameter values that maximize the
marginal log-likelihood. If $N$ is large and the Hessian has full rank, the
second term on the right-hand side can be roughly approximated as
$\frac{1}{2}\log|H|\simeq\frac{n}{2}\log N$, yielding Eq. (5).

An elegant alternative approach to model selection is based on the Akaike
Information Criterion (AIC), which accounts for the fact that the log-
likelihood computed on a training dataset provides an estimate of the
prediction error that is too optimistic, because the same data is being used
to fit the model and assess its error Akaike1974 ; hastie2001elements . To
account for this optimism, a correction term is added which is based on the
sum of covariances between the observed datapoint and the newly generated
datapoint for each input variable $x_{i}$. It can be shown that the sum of
covariances is proportional to the number of degrees of freedom in the
$N\to\infty$ limit, resulting in the following expression for AIC:

$\text{AIC}\equiv-2\log p(y|X,\hat{\theta},K_{i})+2d,$ (7)

where $d$ is the number of degrees of freedom in the model. Thus AIC provides
an estimate of the log-likelihood that would have resulted if another dataset
was independently generated at the same values of input variables (an in-
sample estimate). In the case of GP regression, $d$ needs to be replaced by
$d_{\mathit{eff}}$ in Eq. (7), where $d_{\mathit{eff}}$ is the effective
number of degrees of freedom for the GP regression with a given kernel type,
which captures the amount of smoothing induced by the GP fit gpml ; gam :

$d_{\mathit{eff}}(\hat{\theta})=\mathrm{tr}[K(\hat{\theta})(K(\hat{\theta})+\Sigma)^{-1}],$
(8)

where the dependence of the Gram matrix on the optimal kernel hyperparameters
$\hat{\theta}$ is made explicit for clarity. Note that similar to BIC, lower
AIC values are preferable; however, unlike BIC, AIC is a non-Bayesian measure
and thus provides an alternative approach to model selection.

Choosing the appropriate kernel type is crucial to the success of GP
regression, since different kernels emphasize different correlation structures
in the data. In practice, kernels are often constructed manually using simple
comparison metrics such as marginal likelihood or BIC${}^{\text{naive}}$. In
some cases, composite kernels are constructed automatically using kernel
engineering techniques (see e.g. Ref. duvenaud2013structure ). Here, we
propose a kernel selection technique which is based on the consensus between
AIC and BIC measures of model complexity. This framework allows us to compare
models with different kernels and choose a specific kernel type on the basis
of both a Bayesian approach to model selection, which emphasizes the
complexity of the kernel in terms of the number of kernel parameters, and a
non-Bayesian approach, which is based on the amount of smoothing introduced by
the GP fit.

###  2.3 Poisson likelihood

Since our data consists of integer counts in $N$ bins, we have also employed a
Poisson-type model to generate integer predictions in each bin. Specifically,
we assume that the mean of the predictive probability $f(x_{i})$ (Eq. (3))
provides the rate for the Poisson process in each bin frate2017modeling :

$\log{L_{\cal
P}}=\sum_{i=1}^{N}\left[y_{i}-f(x_{i})-y_{i}\log{\left(\frac{y_{i}}{f(x_{i})}\right)}\right].$
(9)

Note that $f(x_{i})$ implicitly depends on the kernel type and the optimized
hyperparameter values $\hat{\theta}$. Eq. (9) can be used both to generate
integer counts and compute the log-likelihood of the observed counts. Poisson
log-likelihood can also be used instead of the GP marginal log-likelihood to
compute BIC (Eq. (6)), $\text{BIC}^{\text{naive}}$ (Eq. (5)), and AIC (Eq.
(7)).

###  2.4 Gaussian process kernels

We used the GP package from scikit-learn (https://scikit-learn.org/stable/),
augmenting the implementation to include custom kernels and kernel extraction
features. In this paper, we have explored three kernels to model the continuum
background distribution: the Radial Basis Function kernel (RBF), the Matérn
kernel with $\nu=5/2$ (Matern), and the second-order polynomial kernel
(Poly2). The three kernel functions are defined below:

$k_{\text{RBF}}(x,x^{\prime})=\sigma_{0}\exp\left[-\frac{\left(x-x^{\prime}\right)^{2}}{2l^{2}}\right],$
(10)

where $\sigma_{0}$ is the amplitude and $l$ is the length scale of the
covariance function;

$k_{\text{Matern}}(x,x^{\prime})=\sigma_{0}\Bigg{[}1+\frac{\sqrt{5}}{l}d(x,x^{\prime})+\frac{5}{3l}d(x,x^{\prime})^{2}\Bigg{]}\exp\Bigg{[}-\frac{\sqrt{5}}{l}d(x,x^{\prime})\Bigg{]},$
(11)

where $\sigma_{0}$ is the covariance amplitude as in $k_{\text{RBF}}$, $l$ is
a positive parameter characterizing the covariance, and $d(x,x^{\prime})$ is
the Euclidean distance between datapoints $x$ and $x^{\prime}$;

$k_{\text{Poly2}}(x,x^{\prime})=(\sigma_{0}^{2}+x\cdot x^{\prime})^{2},$ (12)

where $\sigma_{0}$ sets the magnitude of the zeroth-order term in the
polynomial expansion. Thus the RBF, Matern and Poly2 kernels depend on 2, 2
and 1 hyperparameter, respectively.

###  2.5 Functional fit

For comparison, we also employ a fourth-order parametric polynomial fit with
explicit basis functions, which is typically used to model the background
distribution ATL-OREACH-PUB-2020-001 :

$f(x_{i})=\sum_{p=0}^{4}w_{p}x_{i}^{p}+m(x_{i}),$ (13)

where $w_{p}$ are the fitting coefficients and $m(x_{i})$ is either set to $0$
for background-only fits with the signal window masked out, or given by Eq.
(4) for signal+background fits on the entire dataset. The fits were carried
out using ROOT data analysis software Brun1997 , by maximizing the Poisson
log-likelihood in Eq. (9). For the background-only fit, the values of the
fitting coefficients are $w_{0}=1.84\times 10^{5}\pm 1.60\times 10^{2}$,
$w_{1}=-4.49\times 10^{3}\pm 1.80\times 10^{0}$, $w_{2}=4.22\times 10^{1}\pm
1.00\times 10^{-2}$, $w_{3}=-1.79\times 10^{-1}\pm 9.00\times 10^{-4}$,
$w_{4}=2.84\times 10^{-4}\pm 4.0\times 10^{-7}$. For the background+signal
fit, the fitting coefficients are $w_{0}=1.64\times 10^{5}\pm 6.16\times
10^{4}$, $w_{1}=-3.87\times 10^{3}\pm 1.86\times 10^{3}$, $w_{2}=3.52\times
10^{1}\pm 2.10\times 10^{1}$, $w_{3}=-1.43\times 10^{-1}\pm 1.05\times
10^{-1}$, $w_{4}=2.19\times 10^{-4}\pm 1.94\times 10^{-4}$, while the signal
contribution is described using $\\{A,\mu,\sigma\\}=\\{443\pm 199,124.5\pm
0.8,2.3\pm 0.9\\}$. All the parameter uncertainties have been estimated via
Hessian analysis available in ROOT. The value of $A$ and its uncertainty have
been rounded to correspond to the integer number of events. The datasets on
which the fits have been performed are described in more detail below.

##  3 Datasets

We use the di-photon sample from the open dataset made available by the ATLAS
collaboration at LHC ATL-OREACH-PUB-2020-001 . We use the selection criteria
as documented in Ref. ATL-OREACH-PUB-2020-001  to create a di-photon invariant
mass distribution, $m_{\gamma\gamma}$, that shows the Higgs decay. The Higgs
decay is a localized bump on top of the smooth background distribution,
traditionally modeled by a polynomial frate2017modeling . The di-photon
distribution consists of integer event counts $y_{i}$ in $N=30$ bins. Since in
this work we focus on the datasets in which we expect to find a localized
signal whose location is approximately known, we first mask out the region
containing the signal. We expect the signal to be localized with a
characteristic width that is small compared to the characteristic length scale
describing the background shape frate2017modeling . In new resonance searches,
we typically scan for the signal at multiple points within the full range of
the dataset, with a prior expectation for the signal width.

To search for a signal in a specific window, we use the entire range of data
with the signal window masked out to determine the optimal parameters for the
background-only GP regression fit. In new resonance searches, this process
could be repeated for multiple masked-out signal windows. Here, we model the
signal using a simple Gaussian whose mean $\mu$ and width $\sigma$ are
approximately known to be 125 GeV and 2.5 GeV, respectively ATL-OREACH-
PUB-2020-001 ; Aad_2012 . Thus, in the background-only fits we mask out a
signal window $\pm 2\sigma$ around the signal mean $\mu$; all data outside of
this window are assumed to belong to the background distribution and are
therefore fit using GP regression with $m(x_{i})=0$. Specifically, we optimize
the parameters $\theta$ of a given kernel by maximizing the marginal log-
likelihood (Eq. (2)), which yields the optimal set of hyperparameters
$\hat{\theta}$. Given $\hat{\theta}$, the predictive distribution is then
provided by Eq. (3) with $m(x_{i})=0$. The resulting hyperparameter values are
$\\{\sigma_{0},l\\}=\\{5.39\times 10^{7},46.0\\}$ for the RBF kernel,
$\\{\sigma_{0},l\\}=\\{1.12\times 10^{8},124.0\\}$ for the Matern kernel, and
$\sigma_{0}=29.3$ for the Poly2 kernel.

##  4 Model selection for background-only fits

To determine which kernel type best represents our data, we have carried out
model selection using BIC and AIC model comparison measures, summarized in
Table 1. We find that the marginal log-likelihood is much worse for Poly2 than
for either RBF or Matern. This disadvantage is too substantial to be offset by
the fact that Poly2 uses one less hyperparameter. As a result, both
$\text{BIC}^{\text{naive}}$ (Eq. (5)) and BIC (Eq. (6)) rank the kernels in
the same way, giving a slight edge to RBF over Matern. Note that this rank is
the same as with the marginal log-likelihood without BIC corrections. However,
Matern is slightly favored over RBF when considering Poisson log-likelihoods
with rates provided by the mean of the GP predictive distribution (Eq. (9)).
This preference for Matern holds when the Poisson log-likelihoods are
augmented with complexity corrections to produce $\text{BIC}^{\text{naive}}$
scores, while Func4 becomes strongly disfavored due to its larger number of
fitting parameters.

To investigate this matter further, we have considered the effect of the AIC
penalty, which effectively accounts for the amount of smoothing effected by
each kernel type gam ; Buja1989 . We observe that RBF is favored over Matern
when the AIC penalty is taken into account (Table 1). Furthermore, Poisson
log-likelihoods slightly favor Func4 over RBF or Matern GP fits. However, this
slight advantage disappears when the AIC correction is taken into account,
with the best score assigned to GP regression with the RBF kernel (Table 1).
Overall, we conclude that with Poisson log-likelihoods, there is a slight
advantage for RBF over Func4 on the basis of AIC and a distinct advantage for
RBF or Matern over Func4 on the basis of $\text{BIC}^{\text{naive}}$.
Considering all the evidence together, it appears that GP regression with the
RBF kernel is the best way to model our data, although the preference of RBF
over Matern is fairly slight.

In addition to the AIC and BIC-based model selection, we have carried out
visual comparisons of the four different models, by plotting the mean
predictive distribution of the GP regression with RBF, Matern and Poly2
kernels ($f(x_{i})$ in Eq. (3) with $m(x_{i})=0$), and the maximum-likelihood
(ML) Func4 fit (Eq. (13)) to the event counts outside of the signal window
(Fig. 1). It is clear from the upper panel of Fig. 1 that RBF, Matern and
Func4 produce very similar fits, whereas Poly2 tends to underfit the data.
This is also clear from the residuals (Fig. 1, lower panel), which are
consistently larger for Poly2 than for the other three models. Interestingly,
RBF, Matern and Func4 all show a slight spurious bump where the models have
been extrapolated across the signal window; outside of the signal window,
deviations from zero are almost always within the error bars. Thus, visual
inspection rules out Poly2 but cannot be used to differentiate between RBF,
Mattern, and Func4.

Table 1: Model comparison for background-only fits. $|H|$ is the determinant of the Hessian of the marginal log-likelihood (Eq. (2)) evaluated at $\hat{\theta}$; $n$ is the number of model parameters; $d$ is the number of model degrees of freedom ($d=d_{\mathit{eff}}$ for GP regression, $d=n$ for Func4); $\log(\text{PL})$ is the Poisson likelihood (Eq. (9)) computed using either $f(x_{i})$ evaluated at the optimal values of hyperparameters $\hat{\theta}$ from GP regression (Poly2,RBF,Matern) or found by non-linear maximization with $f(x_{i})$ given by Eq. (13) (Func4); $\log(\text{GL})$ is the marginal log-likelihood (Eq. (2)). BIC (Eq. (6)) and $\text{BIC}^{\text{naive}}$ (Eq. (5)) were computed using $\log(\text{GL})$ for BIC and $\log(\text{GL})$/$\log(\text{PL})$ for $\text{BIC}^{\text{naive}}$. AIC values (Eq. (7)) were computed using $\log(\text{PL})$.  Model | $\log|H|$ | $n$ | $d$ | -$\log(\text{PL})$ | -$\log(\text{GL})$ | $\text{BIC}^{\text{naive}}_{\text{GL}}$ | $\text{BIC}_{\text{GL}}$ | $\text{AIC}_{\text{PL}}$ | $\text{BIC}^{\text{naive}}_{\text{PL}}$
---|---|---|---|---|---|---|---|---|---
Poly2 | -0.531 | 1 | 2.99 | 38.02 | 87.52 | 89.22 | 87.25 | 82.02 | 39.72
RBF | 0.417 | 2 | 4.68 | 8.95 | 72.15 | 75.55 | 72.36 | 27.26 | 12.35
Matern | 2.906 | 2 | 5.67 | 8.69 | 72.30 | 75.70 | 73.75 | 28.72 | 12.09
Func4 | – | 5 | 5 | 8.65 | – | – | – | 27.30 | 17.15
Figure 1: Background distribution fits. Upper panel: Shown are the mean
predictive distributions predicted by GP regression with RBF, Matern and Poly2
kernels, as well as the ML fit with the fourth-order polynomial (Func4). All
models were fit on the background-only data outside of the $125\pm 5\text{
GeV}$ window where the signal is located. Event counts $y_{i}$ are shown with
error bars constructed using Garwood intervals (black markers: background
only, grey markers: background+signal). Lower panel: Residuals
$y_{i}-f(x_{i})$ for each of the four models in the upper panels, with $y_{i}$
error bars from the upper panel. Residual values are connected by lines of the
same color and type as in the upper panel, to guide the eye. In both panels,
vertical dashed lines indicate the boundaries of the signal window.

##  5 Signal extraction

In order to extract the signal superimposed on top of the background
distribution, we have carried out GP regression with the RBF kernel and
$m(x_{i})$ given by Eq. (4) using the entire dataset (Fig. 2). Importantly,
the kernel parameters were kept at the values $\hat{\theta}$ obtained via the
previous fit to the background distribution, with the signal window masked
out. Thus, the kernel hyperparameters are responsible for modeling the
background, while the Gaussian parameters in Eq. (4) are responsible for
describing the signal. The resulting parameters are
$\\{A_{\text{RBF}},\mu_{\text{RBF}},\sigma_{\text{RBF}}\\}=\\{473\pm
123,124.7\pm 0.6,2.4\pm 0.4\\}$. For comparison, we have also refit the Func4
model (Eq. (13)) on the entire dataset, by adding a Gaussian function (Eq.
(4)) to capture the signal contribution. As mentioned above, the Gaussian
fitting parameters are
$\\{A_{\text{Func4}},\mu_{\text{Func4}},\sigma_{\text{Func4}}\\}=\\{443\pm
199,124.5\pm 0.8,2.3\pm 0.9\\}$ in this case. We observe that both RBF and
Func4 are capable of capturing the approximately Gaussian signal which remains
after subtracting the background distribution. On the basis of the fitted
Gaussian parameters, the correspondence between both models and between each
model and the Higgs simulation (described in Ref. ATL-OREACH-PUB-2020-001 ) is
overall very high (Fig. 2). However, we note that the GP approach with the RBF
kernel extracts a clear Higgs signature consisting of $A_{\text{RBF}}=473\pm
123$ events above the background, compared to $A_{\text{Func4}}=443\pm 199$
events from the Func4 fit. Thus, the mean number of predicted events is higher
and the uncertainty is significantly lower with the GP RBF fit. The lower
uncertainty of the prediction indicates that GP RBF is preferable to Func4.

Figure 2: Extracting the signal from the background distribution. Upper panel:
RBF[Bkg] and Func4[Bkg] are the background fits with the masked-out signal
window (Fig. 1). RBF[Bkg+signal] and Func4[Bkg+signal] are the fits on the
entire dataset, with kernel hyperparameters (RBF) and polynomial parameters
(Func4) kept at their background-only values. Event counts $y_{i}$ are shown
with error bars constructed using Garwood intervals (black markers). Middle
panel: Gaussian signal predicted using Func4 background, residuals
$y_{i}-f(x_{i})$ with respect to the Func4 background-only model, with $y_{i}$
error bars, and the results of the Higgs simulation described in Ref. ATL-
OREACH-PUB-2020-001 . Lower panel: same as the middle panel but with RBF
instead of Func4. In all panels, vertical dashed lines indicate the boundaries
of the signal window.

###  5.1 Synthetic datasets for testing statistical significance of signal
extraction

To investigate the statistical significance of the observed signal, we have
created 500 toy datasets based on the GP fit with the RBF kernel and
$m(x_{i})=0$ to the background-only data. This fit has generated an effective
integer number of counts due to the background only,
$N_{\text{eff}}=[\sum_{i=1}^{N}f(x_{i})]$, where the square brackets indicate
the rounding operation. Next, we sampled from the GP predictive probability
${\cal N}(\tilde{y}_{i}|f(x_{i}),V(x_{i}))$, producing real-valued background
“counts” $\tilde{y}_{i}$ in each bin $i$. Finally, we used
$\tilde{y}_{i}/\sum_{i=1}^{N}\tilde{y}_{i}$ as probabilities in a multinomial
sampling process, generating a synthetic histogram of integer event counts.
Each synthetic histogram was constrained to have $N_{\text{eff}}$ counts,
equal to the total number of events inferred to be due to the background. Note
that our toy datasets include both the uncertainty inherent in GP regression
and the uncertainty related to generating integer event counts from the
underlying model. In order to create a full background+signal test set, we
have added signal counts from the Higgs simulations ATL-OREACH-PUB-2020-001
to each of the 500 background datasets. Thus the signal component is fixed,
while the background component varies from dataset to dataset according to the
background model uncertainties.

###  5.2 Test for biases in signal extraction

To test the robustness of the fit, we check for potential biases in our
background estimation procedure. Namely, for each of the 500 background+signal
toy datasets described above, we carry out a GP fit with the Gaussian mean
function (Eq. (4)) on the entire dataset, while keeping the kernel
hyperparameter values fixed at $\hat{\theta}$, values found by the previously
described fit to the background distribution, with the signal window masked
out. This procedure generates a set of predicted signal strength values,
$\\{A_{\text{pred}}\\}$, which can be compared with the corresponding exact
value, $A_{\text{true}}$, the sum of the event counts added to the background-
only counts in order to create the combined background+signal toy datasets.
Specifically, we compute a Z-score like measure:

$Z_{W}=\frac{A_{\text{pred}}-A_{\text{true}}}{\sigma_{A}},$ (14)

where $A_{\text{pred}}$ and $A_{\text{true}}$ are defined above and
$\sigma_{A}$ is the standard deviation of the $\\{A_{\text{pred}}\\}$ values.

Fig. 3 shows the resulting distribution of the $Z_{W}$ scores. We observe that
the empirical distribution is well described by a Gaussian with $\mu=0.19$ and
$\sigma=1.00$ (the latter is expected due to the normalization in Eq. (14)).
The near-zero value of $\mu$ indicates that there are no substantial biases in
our two-step background+signal reconstruction procedure. We conclude that the
signal contribution can be deconvoluted correctly from the underlying smooth
background distribution.

Figure 3: Test for potential biases in signal extraction. Shown is a
normalized histogram of $Z_{W}$ scores (Eq. (14)) obtained by generating 500
toy datasets and carrying out GP regression as described in the text (blue
bars). Orange curve: a Gaussian fit to the histogram, which yields $\mu=0.19$,
$\sigma=1.00$.

###  5.3 Posterior distributions of signal parameters and significance
analysis

We have investigated the posterior distributions of signal-characterizing
parameters by carrying out Markov Chain Monte Carlo (MCMC) sampling hasting
of the Poisson log-likelihood (Eq. (9)). Routinely employed in Bayesian
analysis, MCMC sampling of posterior probabilities is conceptually similar to
studying model parameter sensitivity and estimating confidence intervals in
frequentist statistics Feldman1998 ; Read_2002 . Poisson rates $f(x_{i})$
depend on the hyperparameter values $\hat{\theta}$ obtained via the previously
described background-only fit and on the mean function $m(x_{i})$, whose
parameters $\\{A,\mu,\sigma\\}$ were sampled from the following priors: the
prior for $A$ is uniform in the $[0,+\infty)$ range, while the prior for $\mu$
is Gaussian, with the $124.7$ GeV mean and $0.02\times 124.7$ GeV standard
deviation. The $\sigma$ prior is also Gaussian, with the $2.4$ GeV mean and
$0.1\times 2.4$ GeV standard deviation. The mean values are consistent with
the Higgs simulation ATL-OREACH-PUB-2020-001  and with the fits presented in
Fig. 2. The $0.02$ and $0.1$ scaling factors in the priors are motivated by
the ATLAS studies Aad_2012 . MCMC was implemented using the Emcee package
emcee  (https://emcee.readthedocs.io), with $10^{4}$ samples in each of 12
independent MC trajectories.

(a)

(b)

Figure 4: Posterior distributions of signal parameters obtained by MCMC
sampling. Shown are posterior probabilities in the corner plot format
(https://corner.readthedocs.io/en/latest/) corner . Diagonal entries:
marginalized posterior probabilities $P(A)$, $P(\mu)$, $P(\sigma)$, off-
diagonal entries: joint posterior probabilities $P(\mu,A)$, $P(\sigma,A)$,
$P(\mu,\sigma)$, with the contours indicating 2D sigma levels. Blue lines mark
the median values in the posterior distributions, yellow lines mark the 95%
quantile values. Panel (a): MCMC analysis of a randomly chosen synthetic
dataset (out of 500 toy datasets with background counts only). Panel (b): MCMC
analysis of the observed counts, which include both background and signal.

Fig. 4 shows MCMC posterior distributions of the three parameters
characterizing the signal: overall signal strength $A$, the mean position of
the signal peak $\mu$, and the width of the signal peak $\sigma$. In Fig. 4a
MCMC sampling was based on a synthetic dataset without any signal added, which
was randomly chosen among the 500 background-only test sets described above.
As expected, $P(A)$, the marginalized posterior probability for signal
strength, is highest when $A$ is close to zero and falls off rapidly as $A$
increases, while $P(\mu)$ and $P(\sigma)$ appear Gaussian. Moreover, the
correlations between all 3 parameter pairs appear to be weak. In contrast,
when the real data is analyzed which contains both the background counts and
Higgs events, the maximum posterior probability value of $A$ is located around
500 counts, consistent with the earlier Hessian analysis of GP regression with
the RBF kernel (Fig. 4b). Indeed, a Gaussian fit of $P(A)$ in Fig. 4b has
yielded $485\pm 121$ Higgs events, very close to the $473\pm 123$ Higgs events
obtained earlier using the GP regression framework. Thus there is a clear
signature of Higgs counts in the real data. Interestingly, the joint
probability $P(\sigma,A)$ reveals a correlation between signal strength and
signal width, with stronger signals tending to have larger widths.

To provide a more quantitative estimate of the statistical significance of the
signal strength observed in real data, we have plotted a histogram of the 95%
confidence levels for $A$ for all 500 background-only toy datasets (Fig. 5).
The value observed with the actual data is 3.02 $\tilde{\sigma}$ above the
median, where $\tilde{\sigma}$ is the distance between the median and the 84%
quantile, and is larger than 99.4% of the values empirically observed in the
histogram.

Figure 5: Distribution of the signal strength $A$ in background-only datasets.
Shown is the histogram of the 95% quantiles (confidence levels, or CL) for $A$
inferred from 500 synthetic datasets with background-only counts using MCMC
sampling. Black dashed line indicates the 50% quantile, or the median, of the
95% CL distribution (i.e., the median number of Higgs events observed above
background, reported at 95% CL). From left to right, the yellow lines show
2.5% and 97.5% quantiles and the green lines show 16% and 84% quantiles,
respectively. The red arrow indicates the 95% CL value obtained from the data
(cf. Fig. 4b).

##  6 Summary

In this work, we have developed a procedure for using Gaussian Process (GP)
regression to extract localized signals from smooth background distributions.
Although this procedure is of interest in many areas of science, including
astrophysics and crystallography, here we focus on extracting Higgs events
from an ATLAS open dataset which consists of binned event counts. Despite its
relatively small size, this is a challenging dataset since the putative signal
is masked by the background and the inference procedure used to analyze the
data affects the statistical significance of the predicted signal.
Traditionally, the background distribution is modeled using a polynomial fit,
onto which a Gaussian signal is superimposed (Eq. (13)) ATL-OREACH-
PUB-2020-001 . Here we propose an alternative framework in which GP regression
is used to model the background and the signal is modeled via the mean
function, which affects both the marginal likelihood of the observed data (Eq.
(2)) and the predictive probability, the conditional probability of a new
datapoint given the previously observed data (Eq. (3)). As with the functional
fits, the mean function is represented by a Gaussian with three free
parameters (Eq. (4)), one of which, $A$, is especially relevant to us since it
represents the total number of signal events found in the dataset.

The GP framework is more flexible than more standard approaches which employ a
fixed set of basis functions such as polynomials or Gaussians
bishop2006pattern ; gpml . This flexibility comes from the focus of the GP-
based approach on the correlation structures in the dataset, which are modeled
using kernel functions. Although GP methods are not limited by the prior
choice of the finite basis (indeed, some popular kernels correspond to
infinite basis sets) and the GP approach is in principle fully Bayesian, most
kernel functions depend on one or several hyperparameters, such as the
characteristic length scale in the RBF kernel. A fully Bayesian treatment of
the dependence of the model evidence on hyperparameters is usually
impractical; however, simply maximizing the evidence with respect to
hyperparameters may lead to overfitting for more complex kernels. In order to
provide a more principled approach to the selection of the kernel type, we
have considered two independent methodologies.

One of them, BIC, is based on evaluating the model evidence under the Laplace
approximation and the assumption that the effects of hyperparameter priors are
negligible (Eq. (6)). With several additional approximations, notably the
assumption that the Hessian matrix has full rank, BIC yields a simple
correction which penalizes model complexity (Eq. (5)). The other approach,
AIC, is non-Bayesian. Instead of concentrating on the model evidence, it
focuses on the degree of smoothing that results from applying a given kernel
to the dataset. Thus, the AIC and BIC approaches are complementary and reflect
different kernel properties (amount of data smoothing vs. the shape of the
log-likelihood landscape as a function of hyperparameters). Using both
criteria holistically, we have chosen a well-known RBF kernel for our GP
regression models, although the results with the Matern kernel are only
slightly worse.

We note that AIC yields approximately equal scores for the GP RBF fit and the
traditional fit, which models the background using a fourth-order polynomial
(Table 1). The results of the two fits are visually very similar when the GP
mean predictive probability is compared with the ML curve produced by the
functional fit, and both approaches are close to the Higgs simulations
predictions (Figs. 1,2). However, the total area $A$ under the signal bump is
somewhat higher with the GP RBF fit compared to the functional fit, with $473$
and $443$ Higgs events, respectively. More critically, Hessian analysis
reveals that the standard deviation is much smaller with the GP prediction,
$123$ vs. $199$ in the functional fit. Thus, the GP approach is preferable
since it leads to the higher signal strength prediction with considerably less
uncertainty.

After ascertaining that our signal extraction procedure is not biased (Fig.
3), we have proceeded to investigate the posterior probabilities of model
parameters by MCMC sampling (Fig. 4). This computational approach is necessary
since we have focused on the Poisson log-likelihood (Eq. (9)), which is more
appropriate for modeling integer event counts. The Poisson log-likelihood
depends on the kernel hyperparameters, which were kept fixed to their values
$\hat{\theta}$ obtained by fitting to the background-only data (Fig. 1), and
on the signal strength, mean and width, which were sampled from prior
distributions. The prior for signal strength $A$ was uninformative, assigning
equal weights to any non-negative value. The priors for the mean and the width
were informative, modeled by Gaussians whose parameters were constrained by
Higgs simulations (Fig. 2) and by the studies of instrumental errors in the
ATLAS detector Aad_2012 . The resulting posterior probability for signal
strength shows a clear Higgs signature, with $485\pm 121$ Higgs events (Fig.
4b). These numbers are consistent with the previous estimate obtained by
Hessian analysis of the signal parameters in GP regression, which yielded
$473\pm 123$ Higgs events.

When the MCMC sampling procedure is applied to synthetic datasets where no
contributions from the signal are expected, the posterior distribution for
signal strength is centered on zero and the typical predicted values are much
smaller (Fig. 4a). The latter is clearly seen by combining the data from $500$
independently generated background-only synthetic datasets into a histogram of
95% confidence levels for signal strength $A$ (Fig. 5). The corresponding
confidence level obtained from the real dataset is larger than $99.4\%$ of the
histogram values and corresponds to $3.02~{}\tilde{\sigma}$, where
$\tilde{\sigma}$ is the distance between the median and the 84% quantile. Thus
our signal strength prediction is also highly significant within the MCMC
framework.

In summary, we have developed a novel GP regression framework for extracting
localized signals from smooth background distributions of unknown functional
form. This problem appears in many areas of science where a weak signal of
interest is masked by background events due to light scattering, extraneous
emission sources, etc. The location and the width of the signal can sometimes
be guessed based on physical considerations; in other cases, consideration of
multiple putative signal windows is necessary, as in LHC anomaly detection
searches ad1 ; ad2 ; ad3 ; ad4 ; ad5 ; ad6 . In both scenarios, only rough
estimates of the position and the width of the signal window are required.
Data outside of the signal window is assumed to belong to the background and a
GP model without the signal contribution is fitted to it. We carry out model
selection using both BIC and AIC considerations, including an in-depth
analysis of the BIC assumptions. The extrapolation of the model across the
signal window then provides an estimate of the background, from which the
signal can now be separated in a second GP fit where only the signal
parameters are allowed to vary, while all the background parameters remain
fixed. This two-step procedure allows us to deconvolute the signal from the
background in a robust and reproducible manner. An application of our approach
to the open Higgs boson dataset from the ATLAS detector (known as the ATLAS
open dataset) yields a highly significant prediction of the Higgs boson
signature, outperforming the traditional approach based on fitting a
polynomial function to the background distribution.

###### Acknowledgements.

This work was supported by the National Science Foundation through grants NSF-
PHY-1607096 (to A.L.) and NSF-MCB-1920914 (to A.V.M.). This manuscript has
been authored by Fermi Research Alliance, LLC under Contract No. DE-
AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office
of High Energy Physics.

## References

  * (1) D.S. Sivia and W.I.F. David, _A Bayesian approach to extracting structure-factor amplitudes from powder diffraction data_ , _Acta Crystallographica Section A_ 50 (1994) 703.
  * (2) W.I.F. David and D.S. Sivia, _Background estimation using a robust Bayesian analysis_ , _Journal of Applied Crystallography_ 34 (2001) 318.
  * (3) T.A. Gordon, E. Agol and D. Foreman-Mackey, _A fast, two-dimensional Gaussian process method based on celerite: Applications to transiting exoplanet discovery and characterization_ , _The Astronomical Journal_ 160 (2020) 240.
  * (4) D. Foreman-Mackey, E. Agol, S. Ambikasaran and R. Angus, _Fast and scalable Gaussian process modeling with applications to astronomical time series_ , _The Astronomical Journal_ 154 (2017) 220.
  * (5) G. Aad, T. Abajyan, B. Abbott, J. Abdallah, S. Abdel Khalek, A. Abdelalim et al., _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Physics Letters B_ 716 (2012) 1.
  * (6) Particle Data Group collaboration, _Review of Particle Physics_ , _Progress of Theoretical and Experimental Physics_ 2020 (2020) 083C01.
  * (7) C. Bishop, _Pattern Recognition and Machine Learning_ , Springer (January, 2006).
  * (8) P. Mehta, M. Bukov, C.-H. Wang, A.G. Day, C. Richardson, C.K. Fisher et al., _A high-bias, low-variance introduction to machine learning for physicists_ , _Physics Reports_ 810 (2019) 1.
  * (9) J. Rocks and P. Mehta, _Memorizing without overfitting: Bias, variance, and interpolation in over-parameterized models_ , arXiv:2010.13933, 2020.
  * (10) S.S. Wilks, _The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses_ , _The Annals of Mathematical Statistics_ 9 (1938) 60.
  * (11) G. Aad, B. Abbott, J. Abdallah, O. Abdinov, B. Abeloos, R. Aben et al., _Search for new phenomena in dijet mass and angular distributions from pp collisions at s=13 TeV with the ATLAS detector_ , _Physics Letters B_ 754 (2016) 302.
  * (12) M. Titsias, _Variational learning of inducing variables in sparse Gaussian processes_ , in _Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics_ , D. van Dyk and M. Welling, eds., vol. 5 of _Proceedings of Machine Learning Research_ , (Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA), pp. 567–574, 16–18 Apr, 2009, http://proceedings.mlr.press/v5/titsias09a.html.
  * (13) C.E. Rasmussen and C.K.I. Williams, _Gaussian Processes for Machine Learning_ , The MIT Press (2005), 10.7551/mitpress/3206.001.0001.
  * (14) K. Kersting, C. Plagemann, P. Pfaff and W. Burgard, _Most likely heteroscedastic Gaussian process regression_ , in _Proceedings of the 24th International Conference on Machine Learning_ , ICML ’07, (New York, NY, USA), pp. 393–400, 2007, DOI.
  * (15) A. O’Hagan, _Curve fitting and optimal design for prediction_ , _Journal of the Royal Statistical Society: Series B (Methodological)_ 40 (1978) 1.
  * (16) K.G. Iyer, E. Gawiser, S.M. Faber, H.C. Ferguson, J. Kartaltepe, A.M. Koekemoer et al., _Nonparametric star formation history reconstruction with Gaussian processes. I. Counting major episodes of star formation_ , _The Astrophysical Journal_ 879 (2019) 116.
  * (17) C.J. Moore, C.P. Berry, A.J. Chua and J.R. Gair, _Improving gravitational wave parameter estimation using Gaussian process regression_ , _Physical Review D_ 93 (2016) 064001.
  * (18) S. Golchi and R. Lockhart, _A Bayesian search for the Higgs particle_ , arXiv:1501.02226, 2015.
  * (19) M. Frate, K. Cranmer, S. Kalia, A. Vandenberg-Rodes and D. Whiteson, _Modeling smooth backgrounds and generic localized signals with Gaussian processes_ , arXiv:1709.05681, 2017.
  * (20) L. Diosan, A. Rogozan and J. Pécuchet, _Evolving kernel functions for svms by genetic programming_ , _Sixth International Conference on Machine Learning and Applications (ICMLA 2007)_ (2007) 19.
  * (21) W. Bing, Z. Wen-qiong, C. Ling and L. Jia-hong, _A GP-based kernel construction and optimization method for RVM_ , in _2010 The 2nd International Conference on Computer and Automation Engineering (ICCAE)_ , vol. 4, pp. 419–423, 2010, DOI.
  * (22) D. Duvenaud, J.R. Lloyd, R. Grosse, J.B. Tenenbaum and Z. Ghahramani, _Structure discovery in nonparametric regression through compositional kernel search_ , arXiv:1302.4922, 2013.
  * (23) A. Wilson and R. Adams, _Gaussian process kernels for pattern discovery and extrapolation_ , in _Proceedings of the 30th International Conference on Machine Learning_ , S. Dasgupta and D. McAllester, eds., vol. 28 of _Proceedings of Machine Learning Research_ , (Atlanta, Georgia, USA), pp. 1067–1075, PMLR, 17–19 Jun, 2013, http://proceedings.mlr.press/v28/wilson13.html.
  * (24) D. Ruppert, M.P. Wand and R.J. Carroll, _Semiparametric Regression (Cambridge Series in Statistical and Probabilistic Mathematics, Series Number 12)_ , Cambridge University Press (2003).
  * (25) ATLAS Collaboration, _Review of the 13 TeV ATLAS open data release_ , Tech. Rep. ATL-OREACH-PUB-2020-001, CERN, Geneva (Jan, 2020).
  * (26) P.W. Goldberg, K.I. Williams and C.M. Bishop, _Regression with input-dependent noise: A Gaussian process treatment_ , _Advances in Neural Information Processing Systems_ 10 (1998) 493.
  * (27) F. Garwood, _Fiducial limits for the Poisson distribution_ , _Biometrika_ 28 (1936) 437.
  * (28) L. Demortier, _Interval estimation_ , in _Data Analysis in High Energy Physics: A Practical Guide to Statistical Methods_ , O. Behnke, K. Kroninger, G. Schott and T. Schorner-Sadenius, eds., (Berlin, Germany), pp. 107–151, Wiley-VCH (2013).
  * (29) G. Schwarz, _Estimating the dimension of a model_ , _The Annals of Statistics_ 6 (1978) 461.
  * (30) H. Akaike, _A new look at the statistical model identification_ , _IEEE Transactions on Automatic Control_ 19 (1974) 716.
  * (31) T. Hastie, R. Tibshirani and J. Friedman, _The Elements of Statistical Learning: Data Mining, Inference, and Prediction_ , Springer series in statistics, Springer (2001).
  * (32) T. Hastie and R. Tibshirani, _Generalized additive models_ , _Statistical Science_ 1 (1986) 297.
  * (33) R. Brun and F. Rademakers, _ROOT: An object oriented data analysis framework_ , _Nuclear Instruments and Methods in Physics A_ 389 (1997) 81.
  * (34) A. Buja, T. Hastie and R. Tibshirani, _Linear smoothers and additive models_ , _The Annals of Statistics_ 17 (1989) 453.
  * (35) W.K. Hastings, _Monte Carlo sampling methods using Markov chains and their applications_ , _Biometrika_ 57 (1970) 97.
  * (36) G.J. Feldman and R.D. Cousins, _Unified approach to the classical statistical analysis of small signals_ , _Physical Review D_ 57 (1998) 3873.
  * (37) A.L. Read, _Presentation of search results: the $CL_{s}$ technique_, _Journal of Physics G: Nuclear and Particle Physics_ 28 (2002) 2693.
  * (38) D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, _emcee: The MCMC Hammer_ , _Publications of the Astronomical Society of the Pacific_ 125 (2013) 306.
  * (39) D. Foreman-Mackey, _corner.py: Scatterplot matrices in Python_ , _The Journal of Open Source Software_ 1 (2016) 24.
  * (40) J. Collins, K. Howe and B. Nachman, _Anomaly detection for resonant new physics with machine learning_ , _Physical Review Letters_ 121 (2018) 241803.
  * (41) T. Heimel, G. Kasieczka, T. Plehn and J. Thompson, _QCD or what?_ , _SciPost Physics_ 6 (2019) 030.
  * (42) P. Jawahar, T. Aarrestad, N. Chernyavskaya, M. Pierini, K.A. Wozniak, J. Ngadiuba et al., _Improving variational autoencoders for new physics detection at the LHC with normalizing flows_ , 2021.
  * (43) O. Amram and C.M. Suarez, _Tag N’ Train: a technique to train improved classifiers on unlabeled data_ , _Journal of High Energy Physics_ 2021 (2021) 153.
  * (44) A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel et al., _Classifying Anomalies THrough Outer Density Estimation (CATHODE)_ , 2021.
  * (45) M. Farina, Y. Nakai and D. Shih, _Searching for new physics with deep autoencoders_ , _Physical Review D_ 101 (2020) 075021.


# Nuclear mass predictions based on deep neural network and finite-range
droplet model (2012)

To Chung Yiu  Department of Physics, The University of Hong Kong, Pokfulam,
999077, Hong Kong  Haozhao Liang  Department of Physics, Graduate School of
Science, The University of Tokyo, Tokyo 113-0033, Japan  Interdisciplinary
Theoretical and Mathematical Sciences Program (iTHEMS), RIKEN, Wako 351-0198,
Japan  Jenny Lee  Department of Physics, The University of Hong Kong,
Pokfulam, 999077, Hong Kong

###### Abstract

A neural network with two hidden layers is developed for nuclear mass
prediction, based on the finite-range droplet model (FRDM12). Different
hyperparameters, including the number of hidden units, the choice of
activation functions, the initializers, and the learning rates, are adjusted
explicitly and systematically. The resulting mass predictions are achieved by
averaging the predictions given by several different sets of hyperparameters
with different regularizers and seed numbers. It can provide us not only the
average values of mass predictions but also reliable estimations in the mass
prediction uncertainties. The overall root-mean-square deviations of nuclear
mass have been reduced from $0.603$ MeV for the FRDM12 model to $0.200$ MeV
and $0.232$ MeV for the training set and validation set, respectively.

nuclear mass, machine learning, deep neural network

##  I Introduction

Nuclear mass is one of the most fundamental nuclear properties [1, 2]. It not
only represents the static properties of nuclei, but also determines the
reaction energies in different nuclear processes, such as $\beta$ decay,
neutron capture, and fission [3]. All these processes play important roles in
the origin of elements and the abundance of elements in the universe [4, 3,
5].

Nowadays, there are around $2500$ nuclei with experimental masses. These
nuclei are estimated to be only $27.8\%$ of around 9000 bounded nuclei
estimated theoretically [6, 2]. In order to have a better understanding of
nuclear mass for all nuclei, theoretical mass models are required. There are
several different theoretical models for mass prediction. Some of them are the
microscopic models such as Hartree–Fock–Bogoliubov (HFB) mass models [7, 8]
and relativistic mean-field (RMF) mass models [9], while the others are the
macroscopic-microscopic models such as finite-range droplet model [6] and
Weizsäcker–Skyrme (WS) model [10]. These theoretical models give a root-mean-
square (RMS) deviation between the theoretical and experimental masses from
$0.3$ to $2.3$ MeV [11]. To reduce RMS deviation, which is critical for a
better description of final element abundance [4, 3], machine learning model
is developed.

Neural network, one of the algorithms in machine learning, has been widely
used in different research fields [12, 13, 14]. From several recent studies,
it has been proven that neural network is able to improve the accuracy of the
models for several different nuclear properties, such as the $\beta$-decay
half-lives [15, 16], neutron capture rate [17], nuclear charge radii [18], and
ground-state and excited energies [19]. Neural network has also been used in
nuclear mass prediction in several previous studies [20, 11, 21, 22]. In
particular, a significant improvement in mass predictions has been achieved by
using several neural networks [22].

Meanwhile, each neural network in these studies includes several kinds of
hyperparameters, e.g., the number of hidden units, the choice of activation
functions, the initializers, and the learning rates [23, 16]. These
hyperparameters play very important roles in both the training process and the
final performance of the neural network. In other words, it is essential to
investigate such hyperparameters in an explicit and systematic way.

In this paper, a deep neural network will be used to build a mass model
improving the current finite-range droplet model [6]. Different
hyperparameters will be in particular investigated in this study to achieve a
better neural network model. Also, several different sets of hyperparameters
will be used together to predict the nuclear masses and provide the
uncertainty of the result.

The details about the neural network used in this study will be given in Sec.
II. The results of the nuclear mass neural network model and its performance
will be discussed in Sec. III. Finally, a summary will be presented in Sec.
IV. For all the hyperparameters adjusted in this work, they will be given in
Appendices A, B, C, and D.

##  II Neural network model

Figure 1: A deep neural network with four inputs, two hidden layers, and one
output. The numbers of hidden units are $j$ and $k$ for the first and second
hidden layers, respectively. For each of the input and hidden layers, there is
a bias term shown with the yellow circle.

In this study, a deep neural network (DNN) consisting of an input layer, two
hidden layers, and an output layer is used, as shown in Fig. 1. The relation
between input, hidden, and output layers are as follows:

$\displaystyle z^{(1)}_{j}$
$\displaystyle=S(\sum_{i}w^{(1)}_{ji}x_{i}+w^{(1)}_{j0}),$ (1) $\displaystyle
z^{(2)}_{k}$
$\displaystyle=\sigma(\sum_{j}w^{(2)}_{kj}z^{(1)}_{j}+w^{(2)}_{k0}),$ (2)
$\displaystyle y_{l}$
$\displaystyle=\sum_{l}w^{(3)}_{lk}z^{(2)}_{k}+w^{(3)}_{l0},$ (3)

where $x_{i}$, $z^{(1)}_{j}$, $z^{(2)}_{k}$, and $y_{l}$ are the inputs, the
hidden units in the first hidden layer, the hidden units in the second hidden
layer, and the output of the model, respectively. Here,
$\boldsymbol{w}=\\{w^{(1)}_{ji},w^{(2)}_{kj},w^{(3)}_{lk},w^{(1)}_{j0},w^{(2)}_{k0},w^{(3)}_{l0}\\}$
are the parameters of the model, $S$ is the sigmoid function [24]

$S(x)=\frac{1}{1+e^{-x}},$ (4)

and $\sigma$ is the softmax function [24]

$\sigma(x_{k})=\frac{e^{x_{k}}}{\sum_{k}e^{x_{k}}},$ (5)

which makes the sum of $z^{(2)}_{k}$ in the second hidden equal to $1$.
Hyperbolic tangent [24] has also been tried as the activation function for the
model,

$\tanh(x)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}.$ (6)

See the key reasons for choosing sigmoid and softmax functions as the
activation functions in Appendix A.

The inputs of the neural network are given by
$\boldsymbol{x}=\\{Z,N,\delta,P\\}$, i.e., including neutron number $N$ and
proton number $Z$, together with the information about nuclear pairing
$\delta$ and shell effect $P$. In particular, the information about nuclear
pairing is characterized by

$\delta=\frac{(-1)^{N}+(-1)^{Z}}{2},$ (7)

so $\delta=1$ for even-even nuclei, $\delta=0$ for odd-$A$ nuclei, and
$\delta=-1$ for odd-odd nuclei. The information about shell effect is
characterized by

$P=\frac{v_{p}v_{n}}{v_{p}+v_{n}},$ (8)

with $v_{p},v_{n}$ being the difference between $Z$, $N$ and the closest magic
number ($2$, $8$, $20$, $28$, $50$, $82$, $126$, or $184$), respectively [25].

In this study, the finite-range droplet model (FRDM12) is used as the
theoretical mass model [6]. FRDM12 is a macroscopic–microscopic model, which
includes the finite-range liquid-drop model in Eq. (II) as the macroscopic
model,

$\displaystyle E_{\rm mac}(A,Z)=$
$\displaystyle~{}a_{v}A+a_{s}A^{2/3}+a_{3}A^{1/3}B_{k}+a_{0}A^{0}$
$\displaystyle+E_{c}-c_{2}Z^{2}A^{1/3}B_{r}-c_{4}\frac{Z^{4/3}}{A^{1/3}}$
$\displaystyle-
c_{5}Z^{2}\frac{B_{w}B_{s}}{B_{1}}+f_{0}\frac{Z^{2}}{A}-c_{a}(Z-N)$
$\displaystyle+W+E_{\rm{pairing}}-a_{el}Z^{2.39},$ (9)

and the folded-Yukawa single-particle potential as the microscopic corrections
[6]. FRDM12 has a relatively small root-mean-square (RMS) deviation between
the experimental and theoretical mass when compared with other theoretical
mass models. Its overall RMS deviation is $0.603$ MeV, while the RMS deviation
of nuclear mass for RMF model is $2.269$ MeV [11]. Besides, the FRDM model has
been used to develop different models for $\beta$-decay properties such as the
FRDM+QRPA model [3]. Therefore, the improved FRDM12 with the neural network
model can be used to improve those $\beta$-decay models.

The target we aim to achieve in this work is the mass residuals between the
theoretical FRDM12 mass $M_{l,\rm{FRDM12}}$ [6] and the experimental mass
$M_{l,\rm{exp}}$ [2],

$t_{l}=M_{l,\rm{FRDM12}}-M_{l,\rm{exp}}.$ (10)

To achieve the target mass residuals from the neural network, different
hyperparameters of the model are required to be adjusted. The mean-square-
error,

$\rm mse(\boldsymbol{x},\boldsymbol{w})=\sum_{l}(t_{l}-y_{l})^{2},$ (11)

is used as the loss function, with $y_{l}$ and $t_{l}$ obtained from Eqs. (3)
and (10), respectively. Adam algorithm [26] with a learning rate of 0.01 is
used to adjust the parameters $\boldsymbol{w}$ to achieve a smaller loss
function. The corresponding key reasons are shown in Appendix B. For the
hidden layers, each of the hidden layers consists of 22 hidden units and a
bias, see the key reasons in Appendix C. Therefore, there are $639$ parameters
in $\boldsymbol{w}$ in total for the model.

For the initial parameter $\boldsymbol{w}$ before the training, there are
three different initializers tried in this study. They include the standard
normal initializer, Glorot normal initializer, and zeros initializer [27]. The
standard normal initializer distributes the initial parameters from a standard
normal distribution with the mean being $0$ and the standard deviation being
$1$. The Glorot normal initializer distributes the initial parameters from a
normal distribution with the mean being $0$ and the standard deviation
equaling to $\sqrt{{2}/(f_{\rm{in}}+f_{\rm{out}})}$. In the above equation,
$f_{\rm in}$ and $f_{\rm out}$ are the numbers of input units and output units
of that layer, respectively [27]. The zeros initializer sets the initial
parameters equal to zero [27]. In this study, the initial parameters
$\boldsymbol{w}$ before the training are generated from the standard normal
distribution for $\\{w^{(1)}_{ji},w^{(2)}_{kj}\\}$, and Glorot normal
distribution for $\\{w^{(1)}_{j0},w^{(2)}_{k0}\\}$, see key reasons in
Appendix D.

The experimental data of this work are from the atomic mass evaluation of 2020
(AME2020) [2]. There are in total $2457$ nuclei with $Z,N\geq 8$. They are
separated into a training set and a validation set randomly. $1966$ nuclei
($80\%$) are randomly selected for the training set, while the remaining $491$
nuclei ($20\%$) are in the validation set.

To achieve a neural network with less variance in the prediction of mass, $33$
sets of hyperparameters with different regularizers and seed numbers are used
for the training. There are three different types of regularizers used in this
study. The first type is without using any regularizer. The second type is the
L2 regularizer. It includes an additional term in the loss function [28],

$L_{\rm{L2}}(\boldsymbol{w})=\lambda\times\sum(w)^{2},$ (12)

with $\lambda$ being a hyperparameter controlling the rate of regularization.
The third type is the orthogonal regularizer. It will encourage the basis of
the output space of that layer to be orthogonal to each other [28]. The
hyperparameter $\lambda$ is used to control the rate of regularization.

In the following calculations, $20000$ epochs are run for each training and
the smallest RMS deviation between $t_{l}$ and $y_{l}$ of the validation set
in different epochs is used to determine the performance of the set of
hyperparameters. Only those sets of hyperparameters giving an RMS deviation of
the validation set smaller than $0.228$ MeV (reduction of RMS deviation over
$60\%$ compared with the FRDM12 prediction) are selected. The average
prediction from the selected sets of hyperparameters is the model output $y$.
Summing the model output y and the theoretical FRDM12 mass $M_{k,\rm{FRDM12}}$
will give us the model mass prediction for the nuclei.

##  III Results and discussion

Figure 2: The RMS deviations of the training and validation sets with
different sets of hyperparameters. Only the seven sets of hyperparameters with
a reduction in RMS deviation larger than $60\%$ are selected as follows: (a)
No regularizer, seed $=3$, (b) L2 regularizer, $\lambda=0.00001$, seed $=1$,
(c) L2 regularizer, $\lambda=0.00001$, seed $=2$, (d) Orthogonal regularizer,
$\lambda=0.1$, seed $=2$, (e) Orthogonal regularizer, $\lambda=0.01$, seed
$=1$, (f) Orthogonal regularizer, $\lambda=0.001$, seed $=1$, and (g)
Orthogonal regularizer, $\lambda=0.001$, seed $=2$. The blue and orange curves
represent the RMS deviations of the training and validation sets,
respectively. The arrow points to the lowest RMS deviation of the validation
set within the $20000$ epochs.

Of the $33$ sets of hyperparameters, seven sets of them are with a reduction
in RMS deviation larger than $60\%$, as shown in Fig. 2. We can see that
different regularizers will speed up or slow down the training of the model,
compared to the model without a regularizer in panel (a). Moreover, since the
seed number will affect the random numbers generated during training and the
initial values of the parameters, different seed numbers are tried for the
model. From Table 1, the RMS derivations have average reductions from $0.603$
MeV to $0.200$ MeV and $0.232$ MeV, having improvements of $66.8\%$ and
$61.5\%$ for the training and validation sets, respectively. It shows that the
neural network can improve the accuracy of the FRDM12 model after adjusting
the hyperparameters.

Table 1: The RMS deviations of nuclear mass between the experimental data from AME2020 [2] and the model predictions. The original RMS deviation between experimental data and FRDM12 [6] is $0.603$ MeV. For each set of hyperparameters, the RMS deviations for the training and validation sets are denoted by $\sigma_{\rm training}$ and $\sigma_{\rm validation}$, respectively. The last two columns show the reduction of RMS deviations from the original deviation $0.603$ MeV to $\sigma_{\rm training}$ and $\sigma_{\rm validation}$, respectively.  Hyperparameters |  $\sigma_{\rm training}$ (MeV) |  $\sigma_{\rm validation}$ (MeV) |  $\Delta\sigma_{\rm training}$ ($\%$) |  $\Delta\sigma_{\rm validation}$ ($\%$)
---|---|---|---|---
No regularizer, seed $=3$ | $0.175$ | $0.228$ | $71.0$ | $62.2$
L2 regularizer, $\lambda=0.00001$, seed $=1$ | $0.204$ | $0.226$ | $66.2$ | $62.5$
L2 regularizer, $\lambda=0.00001$, seed $=2$ | $0.216$ | $0.235$ | $64.2$ | $61.0$
Orthogonal regularizer, $\lambda=0.1$, seed $=2$ | $0.212$ | $0.218$ | $64.8$ | $63.8$
Orthogonal regularizer, $\lambda=0.01$, seed $=1$ | $0.176$ | $0.238$ | $70.8$ | $60.5$
Orthogonal regularizer, $\lambda=0.001$, seed $=1$ | $0.219$ | $0.238$ | $63.7$ | $60.5$
Orthogonal regularizer, $\lambda=0.001$, seed $=2$ | $0.201$ | $0.238$ | $66.7$ | $60.5$
Average of the above sets | $0.200$ | $0.232$ | $66.8$ | $61.5$

For most of the isotopes and isotones, the neural network model can give a
good performance in mass prediction. For example, from Fig. 3, the average of
the selected seven sets of hyperparameters has a much better overall
performance in the mass difference compared to the FRDM12 model for the Nd
isotopes and $N=66$ isotones. The data in the training set is represented in
the gray-hatched regions in Fig. 3. It is shown that the mass prediction in
the training set is, in general, better than that in the validation set. For
the regions far away from the training set, such as $N<68$ for Nd isotopes and
$Z>59$ for $N=66$ isotones, the mass difference is still within the mass
uncertainty as shown in the yellow-hatched regions, so the neural network
model is still with a good performance in those regions. However, the error
bar for the neural network model is larger in regions far away from the
training set. It shows that the mass prediction for different sets of
hyperparameters in those regions is with large variation. Therefore, multiple
sets of hyperparameters are selected, and averaging the prediction of the
seven sets of hyperparameters can provide us with the uncertainty of nuclear
mass prediction for different nuclei. The larger the error bar for the nuclei,
the less confidence we have in the mass prediction.

Figure 3: Mass difference between the DNN mass data and the experimental data,
$t_{k}=M_{k,DNN}-M_{k,exp}$, for (a) Nd isotopes and (b) $N=66$ isotones. The
blue curves represent the mass difference between the FRDM12 model and the
experimental data. The semi-transparent curves in the background represent the
mass difference between the FRDM12+DNN mass model and the experimental data
with the seven sets of hyperparameters, respectively. The red curves with
uncertainties represent the mass difference for the average of the FRDM12+DNN
mass model. The error bars are calculated from the standard deviation of the
mass difference from the $7$ sets of hyperparameters. The gray-hatched regions
represent the training set. The yellow-hatched regions represent the mass
uncertainties by including the average RMS deviations of seven sets of
hyperparameters together with the experimental uncertainties. Figure 4: Same
as Fig. 3, but for Cr isotopes.

Furthermore, there are several nuclei in which the masses from the FRDM12
model are significantly different from the experimental data. For example, the
mass difference can be up to $2.5$ MeV. It is interesting to see whether or
not the neural network can significantly improve the mass predictions for
those specific nuclei.

For Cr isotopes shown in Fig. 4, FRDM12 cannot give a good mass prediction for
the region $N<22$. In this region, the neural network model can still provide
a better mass prediction compared with the FRDM12 model. However, the mass
difference of the two nuclei $N=17$ and $18$ is not within the yellow-hatched
region. Although the mass prediction for the two nuclei $N=17$ and $18$ is not
as good as other nuclei, the neural network model can improve the mass
difference from $-2.5$ MeV to $-2$ MeV and from $-2.1$ MeV to $-0.6$ MeV,
respectively. It shows that the neural network model can improve mass
prediction, even though the original theoretical model cannot provide an
accurate prediction.

Figure 5: Same as Fig. 3, but for Ba isotopes.

Finally, we also need to see if there are any nuclei in which the mass
prediction becomes worse after applying the neural network model. For Ba
isotopes in Fig. 5, it shows a good agreement in the mass prediction for the
training set and large-neutron-number regions. However, for $N=60$ and $62$,
the neural network gives worse predictions compared to the FRDM12 model. It
may happen because of two reasons. First, the nuclei $N=60$ and $62$ are in
the validation set, which are not involved in the training of the model.
Therefore, the overall performance of the nuclei in the validation set is not
as good as the nuclei in the training set. Second, this neural network model
still cannot tackle the odd-even staggering in light nuclei. Therefore, these
two even-even nuclei do not have a better mass prediction after using the
neural network, while the nearby even-odd nuclei in the training set have a
better performance.

##  IV Summary and perspectives

The deep neural network has been applied to study the nuclear mass based on
the FRDM12 model. Different hyperparameters, including activation functions,
learning rate, the number of hidden units, and the initializers, are adjusted
in a systematic way to achieve a better performance of the model, as shown in
Appendices A, B, C, and D. Finally, seven sets of hyperparameters with
different regularizers and seed numbers have been selected. It is important to
find that averaging the predictions given by several different sets of
hyperparameters can provide us not only the average values of mass predictions
but also reliable estimations in the mass prediction uncertainties.

With the neural network, the RMS deviations between the experimental mass and
theoretical mass have been reduced from $0.603$ MeV to $0.200$ MeV and $0.232$
MeV for the training set and validation set, respectively. For most of the
nuclei, this DNN model can give a better mass prediction compared with the
FRDM12 model. Even for those nuclei which have a poor mass prediction in the
FRDM12 model such as 41Cr and 42Cr, the DNN model can still reduce the RMS
deviation and achieve a better mass prediction. However, there are still
several nuclei with worse mass prediction compared with the FRDM12 model. It
may happen in the validation set with unsolved odd-even staggering problems.
Further studies for those nuclei are required in the future to improve the
model, such as a better description of the odd-even staggering in nuclei.

In the future, the same technique can be applied to other physics quantities
related to nuclei. With the nuclear mass predicted in this study and adjusting
different hyperparameters, DNN models for different physics quantities, such
as $\beta$-decay half-life and $\beta$-delayed neutron emission probability,
can be generated and their performance can be compared with other current
theoretical models to investigate if this algorithm can improve the present
predictions.

###### Acknowledgements.

H.L. acknowledges the JSPS Grant-in-Aid for Early-Career Scientists under
Grant No. 18K13549, the JSPS Grant-in-Aid for Scientific Research (S) under
Grant No. 20H05648, and the RIKEN Pioneering Project: Evolution of Matter in
the Universe. J.L. acknowledges the Research Grants Council (RGC) of Hong Kong
with a grant of General Research Funding (GRF-17312522).

##  Appendix A Performance with different activation function

In this appendix, we investigate explicitly the impact of the activation
functions on the performance of the training.

Figure 6: The RMS deviations of the training and validation sets in $20000$
epochs with $20$ hidden units in both hidden layers, and no activation
function. The blue and orange curves represent the RMS deviations of the
training and validation sets, respectively.

The neural network model is able to approximate non-linear relations between
inputs and outputs because of the use of the activation functions. Without
using an activation function, the model is just a linear regression that
cannot approximate non-linear relations. Figure 6 shows that, without the
activation function, the RMS deviation of the model will remain at around
$0.6$ MeV and no improvement can be made at all.

Figure 7: The RMS deviations of the training and validation sets in $20000$
epochs with $20$ hidden units in both hidden layers. The activation functions
for the two hidden layers are selected to be (a) (tanh, tanh), (b) (tanh,
sigmoid), (c) (tanh, softmax), (d) (sigmoid, tanh), (e) (sigmoid, sigmoid),
(f) (sigmoid, softmax), (g) (softmax, tanh), (h) (softmax, sigmoid), and (i)
(softmax, softmax), respectively. The blue and orange curves represent the RMS
deviations of the training and validation sets, respectively. Figure 8: Same
as Fig. 7, but the number of hidden units in both hidden layers is $22$.

Therefore, models with different activation functions are trained in this
study. From Fig. 7 with $20$ hidden units, the RMS deviations of panels (a),
(c), (h), and (i) are highly oscillating, and thus the corresponding
activation functions are not recommended. In contrast, panels (b), (d), (e),
(f), and (g) show a steady RMS deviation without large oscillation.

To find out the best activation functions for the present mass models, results
from Fig. 8 with $22$ hidden units are examined. Panel (b) is oscillating and
not recommended. Panels (e) and (g) show RMS deviations larger than $0.3$ MeV,
which is larger than the RMS deviations of panels (d) and (f), so panels (e)
and (g) are not recommended. Panels (d) and (f) show small RMS derivations and
less oscillation, so they are recommended for this mass model. In this study,
the activation functions for the two hidden layers are chosen to be (f)
sigmoid and softmax functions for all the following calculations,
respectively. If there are any future studies, the activation functions can be
chosen as (d) sigmoid and tanh functions, respectively.

##  Appendix B Performance with different learning rates

In this appendix, we investigate explicitly the impact of the learning rate in
the algorithm of the model on the performance of the training.

Figure 9: The RMS deviations of the training and validation sets in $20000$
epochs using Adam algorithm as the optimizer with the learning rate (a) $0.1$,
(b) $0.01$, (c) $0.001$, and (d) $0.0001$. The blue and orange curves
represent the RMS deviations of the training and validation sets,
respectively.

In Fig. 9, the RMS deviations of the training and validation sets using the
Adam algorithm with the learning rates $0.1$, $0.01$, $0.001$, and $0.0001$
are shown. It is seen in Fig. 9(a) that a model with a large learning rate
makes the model parameters oscillate and is unable to reach the minimum, so
its RMS deviations are highly oscillating with large values. A model with a
smaller learning rate makes the number of epochs used to reach the lowest RMS
deviation larger, and it is more time-consuming. Using the learning rate
$0.01$, the RMS deviation for the validation set reaches around $0.27$ MeV
after $20000$ epochs as shown in Fig. 9(b). In contrast, using the learning
rates $0.001$ and $0.0001$, the RMS deviations for the validation set are
still around $0.43$ and $0.54$ MeV after $20000$ epochs as shown in Figs. 9(c)
and 9(d), respectively. Therefore, the learning rate is selected as $0.01$ for
all the following calculations in this study.

##  Appendix C Performance with different number of hidden units

In this appendix, we investigate explicitly the impact of the number of hidden
units on the performance of the training.

Figure 10: The RMS deviations of the training and validation sets in $20000$
epochs with different numbers of hidden units per hidden layer within
$10$–$70$: (a) $10$, (b) $20$, (c) $30$, (d) $40$, (e) $50$, (f) $60$, and (g)
$70$. The blue and orange curves represent the RMS deviations of the training
and validation sets, respectively.

In Fig. 10, the RMS deviations of the training and validation set with
different numbers of hidden units per hidden layer are shown, where the
numbers of hidden units per hidden layer are changed from $10$ to $70$.
Generally speaking, a model with a small number of hidden units may make the
number of parameters not enough for the model and be unable to reach a small
RMS deviation. It is confirmed by Fig. 10(a) that the RMS deviation for the
validation set for $10$ hidden units per hidden layer remains around $0.37$
MeV after $10000$ epochs.

For the model with a larger number of hidden units, the number of parameters
will become larger, and the time used for training will become longer.
Moreover, since the data in the training set in general include several kinds
of noise, the large number of hidden units may lead to the overfitting problem
and cause a poor prediction outside the training set. As shown in panels (e),
(f), and (g), although the RMS deviations for the training set are still
decreasing after $10000$ epochs, the RMS deviations for the validation set are
increasing in that region, resulting in a poor prediction of mass.

On the other hand, panel (b) shows less oscillation in the RMS deviation
compared to panels (c) and (d). Meanwhile, $30$ or more hidden units show no
improvement in RMS deviation compared with $20$ hidden units. For the above
reasons, the number of hidden units around $20$ is selected.

Figure 11: The RMS deviation of the training and validation set in $20000$
epochs with different numbers of hidden units per hidden layer within
$16$–$24$: (a) $16$, (b) $17$, (c) $18$, (d) $19$, (e) $20$, (f) $21$, (g)
$22$, (h) $23$, and (i) $24$. The blue and orange curves represent the RMS
deviations of the training and validation sets, respectively.

To find a more precise number of hidden units, models with different numbers
of hidden units from $16$ to $24$ are trained and their performance is shown
in Fig. 11. It is seen that the performance of the models with different
numbers of hidden units is similar, with around $0.25$ MeV in the RMS
deviation of the validation set. However, it shows that, for models with $23$
and $24$ hidden units per hidden layer, the RMS deviations are oscillating, so
the models with $23$ and $24$ hidden units per hidden layer are not
recommended. For the remaining models with similar performance, the number of
hidden units per hidden layer is chosen to be $22$, as it can provide more
parameters without large oscillation or overfitting.

##  Appendix D Performance with different initializers

In this appendix, we investigate explicitly the impact of the initializers on
the performance of the training.

Figure 12: The RMS deviations of the training and validation sets in $20000$
epochs with different initializers for the parameters
$\\{w^{(1)}_{ji},w^{(2)}_{kj}\\}$ and bias $\\{w^{(1)}_{j0},w^{(2)}_{k0}\\}$.
The initializers adopted here include the standard normal, glorot normal, and
zeros distributions [27]. The blue and orange curves represent the RMS
deviations of the training and validation sets, respectively.

The initial values of the parameters also affect the performance of the model.
For example, when all the parameters are initialized with the same value, all
the hidden units in the same hidden layer will be identical to each other.
Such kind of model is equivalent to a model with 1 hidden unit in each hidden
layer. The small number of hidden units makes such a model with poor
performance and cannot predict the output. As shown in Fig. 12(i), the RMS
deviations are remaining at around $0.575$ MeV during the training.

Therefore, we train the models with different sets of initial parameters
($\\{w^{(1)}_{ji},w^{(2)}_{kj}\\}$, $\\{w^{(1)}_{j0},w^{(2)}_{k0}\\}$). For
different sets of initial parameters trained, panels (a),(b), (c), (f), (g),
and (h) give the RMS deviations of the validation set around $0.23$ MeV. Among
these 6 sets of initial parameters, panels (a), (b), and (f) give a steady RMS
deviation curve during the whole training. Therefore, these three sets of
initial parameters are recommended. In this work, initial parameters from
panel (b) are selected for all the following calculations. If there are any
future studies, panels (a) and (f) can be used as the initial parameters.

## References

  * Lunney  _et al._ [2003] D. Lunney, J. M. Pearson, and C. Thibault, Recent trends in the determination of nuclear masses, Reviews of Modern Physics 75, 1021 (2003).
  * Huang  _et al._ [2021] W. Huang, M. Wang, F. Kondev, G. Audi, and S. Naimi, The ame 2020 atomic mass evaluation (i). evaluation of input data, and adjustment procedures*, Chinese Physics C 45, 030002 (2021).
  * Mumpower  _et al._ [2016] M. Mumpower, R. Surman, G. McLaughlin, and A. Aprahamian, The impact of individual nuclear properties on r-process nucleosynthesis, Progress in Particle and Nuclear Physics 86, 86 (2016).
  * Burbidge  _et al._ [1957] E. M. Burbidge, G. R. Burbidge, W. A. Fowler, and F. Hoyle, Synthesis of the elements in stars, Reviews of Modern Physics 29, 547 (1957).
  * Martin  _et al._ [2016] D. Martin, A. Arcones, W. Nazarewicz, and E. Olsen, Impact of nuclear mass uncertainties on the $r$ process, Physical Review Letters 116, 121101 (2016).
  * Möller  _et al._ [2016] P. Möller, A. Sierk, T. Ichikawa, and H. Sagawa, Nuclear ground-state masses and deformations: Frdm(2012), Atomic Data and Nuclear Data Tables 109-110, 1 (2016).
  * Goriely  _et al._ [2009] S. Goriely, N. Chamel, and J. M. Pearson, Skyrme-hartree-fock-bogoliubov nuclear mass formulas: Crossing the 0.6 mev accuracy threshold with microscopically deduced pairing, Physical Review Letters 102, 152503 (2009).
  * Goriely  _et al._ [2013] S. Goriely, N. Chamel, and J. M. Pearson, Further explorations of skyrme-hartree-fock-bogoliubov mass formulas. xiii. the 2012 atomic mass evaluation and the symmetry coefficient, Physical Review C 88, 024308 (2013).
  * Zhao  _et al._ [2010] P. W. Zhao, Z. P. Li, J. M. Yao, and J. Meng, New parametrization for the nuclear covariant energy density functional with a point-coupling interaction, Physical Review C 82, 054319 (2010).
  * Wang  _et al._ [2014] N. Wang, M. Liu, X. Wu, and J. Meng, Surface diffuseness correction in global mass formula, Physics Letters B 734, 215 (2014).
  * Niu and Liang [2018] Z. Niu and H. Liang, Nuclear mass predictions based on bayesian neural network approach with pairing and shell effects, Physics Letters B 778, 48 (2018).
  * Pederson  _et al._ [2022] R. Pederson, B. Kalita, and K. Burke, Machine learning and density functional theory, Nature Reviews Physics 4, 357 (2022).
  * Kulik  _et al._ [2022] H. J. Kulik, T. Hammerschmidt, J. Schmidt, S. Botti, M. A. L. Marques, M. Boley, M. Scheffler, M. Todorović, P. Rinke, C. Oses, A. Smolyanyuk, S. Curtarolo, A. Tkatchenko, A. P. Bartók, S. Manzhos, M. Ihara, T. Carrington, J. Behler, O. Isayev, M. Veit, A. Grisafi, J. Nigam, M. Ceriotti, K. T. Schütt, J. Westermayr, M. Gastegger, R. J. Maurer, B. Kalita, K. Burke, R. Nagai, R. Akashi, O. Sugino, J. Hermann, F. Noé, S. Pilati, C. Draxl, M. Kuban, S. Rigamonti, M. Scheidgen, M. Esters, D. Hicks, C. Toher, P. V. Balachandran, I. Tamblyn, S. Whitelam, C. Bellinger, and L. M. Ghiringhelli, Roadmap on machine learning in electronic structure, Electronic Structure 4, 023004 (2022).
  * Boehnlein  _et al._ [2022] A. Boehnlein, M. Diefenthaler, N. Sato, M. Schram, V. Ziegler, C. Fanelli, M. Hjorth-Jensen, T. Horn, M. P. Kuchera, D. Lee, W. Nazarewicz, P. Ostroumov, K. Orginos, A. Poon, X.-N. Wang, A. Scheinker, M. S. Smith, and L.-G. Pang, Colloquium: Machine learning in nuclear physics, Reviews of Modern Physics 94, 031003 (2022).
  * Niu  _et al._ [2019] Z. M. Niu, H. Z. Liang, B. H. Sun, W. H. Long, and Y. F. Niu, Predictions of nuclear $\beta$-decay half-lives with machine learning and their impact on $r$-process nucleosynthesis, Physical Review C 99, 064307 (2019).
  * Minato  _et al._ [2022] F. Minato, Z. Niu, and H. Liang, Calculation of $\beta$-decay half-lives within a skyrme-hartree-fock-bogoliubov energy density functional with the proton-neutron quasiparticle random-phase approximation and isoscalar pairing strengths optimized by a bayesian method, Physical Review C 106, 024306 (2022).
  * Ma  _et al._ [2019] C. Ma, Z. Li, Z. M. Niu, and H. Z. Liang, Influence of nuclear mass uncertainties on radiative neutron-capture rates, Physical Review C 100, 024330 (2019).
  * Wu  _et al._ [2020] D. Wu, C. L. Bai, H. Sagawa, and H. Q. Zhang, Calculation of nuclear charge radii with a trained feed-forward neural network, Physical Review C 102, 054323 (2020).
  * Lasseri  _et al._ [2020] R.-D. Lasseri, D. Regnier, J.-P. Ebran, and A. Penon, Taming nuclear complexity with a committee of multilayer neural networks, Physical Review Letters 124, 162502 (2020).
  * Utama and Piekarewicz [2017] R. Utama and J. Piekarewicz, Refining mass formulas for astrophysical applications: A bayesian neural network approach, Physical Review C 96, 044308 (2017).
  * Yüksel  _et al._ [2021] E. Yüksel, D. Soydaner, and H. Bahtiyar, Nuclear binding energy predictions using neural networks: Application of the multilayer perceptron, International Journal of Modern Physics E: Nuclear Physics 30, 2150017 (2021).
  * Niu and Liang [2022] Z. M. Niu and H. Z. Liang, Nuclear mass predictions with machine learning reaching the accuracy required by $r$-process studies, Physical Review C 106, L021303 (2022).
  * Wu  _et al._ [2021] D. Wu, C. L. Bai, H. Sagawa, S. Nishimura, and H. Q. Zhang, $\beta$-delayed one-neutron emission probabilities within a neural network model, Physical Review C 104, 054303 (2021).
  * TensorFlow [2023a] TensorFlow, Module: tf.keras.activations, https://www.tensorflow.org/api_docs/python/tf/keras/activations (2023a) .
  * Kirson [2008] M. W. Kirson, Mutual influence of terms in a semi-empirical mass formula, Nuclear Physics A 798, 29 (2008).
  * TensorFlow [2023b] TensorFlow, tf.keras.optimizers.adam, https://www.tensorflow.org/api_docs/python/tf/keras/optimizers/Adam (2023b) .
  * TensorFlow [2023c] TensorFlow, Module: tf.keras.initializers, https://www.tensorflow.org/api_docs/python/tf/keras/initializers (2023c) .
  * TensorFlow [2023d] TensorFlow, Module: tf.keras.regularizers, https://www.tensorflow.org/api_docs/python/tf/keras/regularizers (2023d) .