# Unfolding with Gaussian Processes

Adam Bozson  Glen Cowan  Francesco Spanò  Department of Physics
Royal Holloway, University of London
Egham, Surrey, TW20 0EX, United Kingdom

###### Abstract

A method to perform unfolding with Gaussian processes (GPs) is presented.
Using Bayesian regression, we define an estimator for the underlying truth
distribution as the mode of the posterior. We show that in the case where the
bin contents are distributed approximately according to a Gaussian, this
estimator is equivalent to the mean function of a GP conditioned on the
maximum likelihood estimator. Regularisation is introduced via the kernel
function of the GP, which has a natural interpretation as the covariance of
the underlying distribution. This novel approach allows for the regularisation
to be informed by prior knowledge of the underlying distribution, and for it
to be varied along the spectrum. In addition, the full statistical covariance
matrix for the estimator is obtained as part of the result. The method is
applied to two examples: a double-peaked bimodal distribution and a falling
spectrum.

###### keywords:

unfolding , Gaussian process

††journal: NIM A

##  1 Introduction

Experimental measurements are distorted and biased by detector effects, due to
limitations of the measuring instrument and procedures. The need to infer the
underlying distribution using the measured data is shared by variety of
fields, from astronomy [1] and medical applications [2] to the investigation
of the parameters that describe oil well properties [3].

In most of these fields, these techniques are called _deconvolution_ or
_restoration_ [4]. They are used to solve what is defined as the inverse
problem: to infer an unknown function $f(\bm{x})$ from the measured data,
using knowledge and assumptions of the distortions.

In particle physics such techniques are known as _unfolding_ and a variety of
methods have been developed for this purpose (for some reviews see Refs. [5,
6, 7]).

In this paper, a novel Bayesian method to perform unfolding in particle
physics is proposed. We use an approach that “gives prior probability to every
possible function” via _Gaussian process regression_ [8], where higher
probabilities are assigned to functions that are considered to agree with the
observations. This approach allows greater flexibility than unfolding schemes
based on a set of parametrised functions belonging to a specific class. In
addition, it is shown to have a locally tunable regularisation scheme in terms
of the variable to be unfolded.

In Sec. 2, we define the unfolding problem and the notation for approximately
Gaussian-distributed datasets. Sec. 3 discusses the solution to the unfolding
problem based on the maximum likelihood (ML) method, and the need for
regularisation. In a Bayesian setting, the likelihood is enhanced by prior
information so that the ML solution is replaced by the mode of the posterior
distribution. Sec. 4 connects the maximum _a posteriori_ (MAP) estimator to
the solution of a regression problem which conditions prior knowledge encoded
in a Gaussian process on the ML solution extracted from data. Example
applications are provided in Sec. 5. Finally, we report the conclusions and
outlook for future exploration of this method in Sec. 6.

##  2 Definitions and notation

In particle physics, measured distributions are often reported as populations
of bins rather than continuous functions. Therefore the first step we will
take is to represent the underlying distributions with discretised bin
populations. We note that this process biases the estimated histogram away
from the true distribution.

The truth distribution is referred to as $f(\bm{x})$ and represented by a
histogram $\bm{\mu}=(\mu_{1},\ldots,\mu_{M})$ with contents
$\mu_{j}\propto\int_{\text{bin\,}j}f(\bm{x})\operatorname{d\\!}\bm{x}$,
$j=1,\ldots,M$. Observed data are contained in a histogram
$\bm{n}=(n_{1},\ldots,n_{N})$ with $N$ bins. The expectation value of $\bm{n}$
is the histogram $\bm{\nu}$.

The truth and observed distributions are related through the effects of
detector response, acceptance, and background contributions. For simplicity,
we take the background to be zero (the relaxation of this assumption is
discussed in Sec. 6). The contents of $\bm{\mu}$ and $\bm{\nu}$ are linearly
related by

$\bm{\nu}=R\bm{\mu},$ (1)

where $R$ is the $N\times M$ _response matrix_ with elements $R_{ij}$ giving
the conditional probability for an event to lie in bin $i$ of the observed
histogram, given its true value is in bin $j$ of the truth histogram.

The goal of unfolding is to construct estimators $\hat{\bm{\mu}}$ for the
truth histogram $\bm{\mu}$, along with their covariance matrix
$U_{ij}=\operatorname{cov}[\hat{\mu}_{i},\hat{\mu}_{j}]$. In an experiment,
the bin counts of the observed histogram $n_{i}$ fluctuate according to the
Poisson distribution with expectation values $\nu_{i}$ and covariance matrix
$V_{ij}=\nu_{i}\,\delta_{ij}$. The findings in this paper apply when bin
counts are approximately Gaussian, i.e., for large $\nu_{i}$.

##  3 ML solution and regularisation

Since the data $\bm{n}$ are approximately Gaussian-distributed around
$\bm{\nu}$, the likelihood is given by

$P(\bm{n}|\bm{\nu})=\left[(2\pi)^{N}\det(V)\right]^{-\frac{1}{2}}\exp\left[-\frac{1}{2}(\bm{n}-\bm{\nu})^{\mathsf{T}}V^{-1}(\bm{n}-\bm{\nu})\right]$
(2)

and hence the log-likelihood may be written

$\log
P(\bm{n}|\bm{\mu})=-\frac{1}{2}\left(\bm{n}-R\bm{\mu}\right)^{\mathsf{T}}V^{-1}\left(\bm{n}-R\bm{\mu}\right)+\ldots,$
(3)

where we have substituted $\bm{\nu}=R\bm{\mu}$ and not written terms that do
not depend on $\bm{\mu}$. It can be shown that the maximum likelihood (ML)
solution $\hat{\bm{\mu}}_{\text{ML}}$ satisfies
$\bm{n}=R\hat{\bm{\mu}}_{\text{ML}}$. We write the ML solution as
$\hat{\bm{\mu}}_{\text{ML}}=R^{-1}\bm{n}$, which may be obtained by explicit
matrix inversion for invertible $R$ when $N=M$ or by alternative methods, such
as numerically maximising Eq. (3) or singular value decomposition. The ML
covariance matrix is given by
$U_{\text{ML}}=R^{-1}V\left(R^{-1}\right)^{\textsf{T}}$ [9].

The detector response acts to smear out fine structure in the truth
distribution, so statistical fluctuations in the data can lead to a large
amount of fine structure in the unfolded result. This effect yields large
local fluctuations in the ML unfolded solution when the typical bin width is
not much larger than the detector resolution. In addition, the estimators for
neighbouring bin counts can often have strong negative correlations.

These unwanted false features are typically reduced by a technique known as
_regularisation_. Regularisation may be introduced by minimising a _cost
functional_ $\Phi(\bm{\mu})=-\alpha\log P(\bm{n}|\bm{\mu})+S(\bm{\mu})$, where
$S(\bm{\mu})$ penalises high-variance distributions, effectively constricting
the space of possible unfolded solutions. Multiple measures of smoothness may
be used, such as those based on derivatives [10, 11] or entropy [12]. The ML
solution has the minimum variance for an unbiased estimator, so any reduction
in variance must be balanced by introducing some bias. The regularisation
parameter $\alpha$ controls this bias–variance trade-off.

An unfolded distribution may alternatively be obtained by iterative techniques
[13, 14], which converge on the ML solution. Stopping after a fixed number of
iterations can yield a solution with the desired properties, although the fact
that the bias–variance trade-off is controlled by a discrete parameter, rather
than a continuous one, is seen as a disadvantage because it limits the
possibility to tune the parameter values. _Fully Bayesian unfolding_ [15]
addresses regularisation through a non-constant prior distribution, and
performs the unfolding by sampling from the posterior distribution.

##  4 Gaussian process method

The method presented in this paper builds on the ML solution from Sec. 3.
Starting from the Gaussian likelihood given by Eq. (3) and a _Gaussian
process_ (GP) prior, Bayes’ theorem is applied to obtain a posterior
distribution. We define the estimator representing the unfolded distribution
as a summary statistic of the posterior, namely the mode. We remark that while
the ML solution is a frequentist estimator, the method presented here
incorporates elements of Bayesian statistics. However the final estimator for
the unfolded distribution is a valid frequentist estimator, so this method may
be used in either fully Bayesian or hybrid analyses.

A _random process_ extends the notion of a random variable to the space of
functions of a set of indices $X$. A GP is therefore a set of indexed random
variables, any finite subset of which are distributed according to a joint
Gaussian distribution [8]. Since the Gaussian distribution is entirely defined
by its mean and covariance, a complete description of a GP requires just a
mean function $m(\bm{x})=\operatorname{E}[f(\bm{x})]$ and a _kernel function_
$k(\bm{x},\bm{x}^{\prime})=\operatorname{cov}[f(\bm{x}),f(\bm{x}^{\prime})]$.
A GP can be thought of as a probability distribution for the _latent function_
$f(\bm{x})$.

GPs may be used for regression, where one wishes to estimate the function $f$
of the (generally multidimensional) variable $\bm{x}$, given some observations
$\bm{y}$ taken at $X=(\bm{x}_{1},\allowbreak\bm{x}_{2},\ldots)$. This is done
by updating a GP prior using Bayes’ theorem to obtain a posterior GP for $f$.
The mean (or equivalently the mode) of the posterior, evaluated at
$X_{*}=(\bm{x}_{*1},\bm{x}_{*2},\ldots)$, is denoted $\bar{\bm{f}}_{*}$ and
used as the estimator for $f$. A rich treatment of using GPs for regression
may be found in Ch. 2 of Ref. [8], whose notation we follow in this paper.
Here we state the posterior mean and covariance of a GP with prior mean
function $m$ and kernel function $k$ for a vector of observations $\bm{y}$
with data covariance matrix $V$:

$\displaystyle\bar{\bm{f}}_{*}$
$\displaystyle=K_{*}^{\mathsf{T}}\left[K+V\right]^{-1}\bm{y}+\bm{m}_{*},$ (4)
$\displaystyle\operatorname{cov}(\bm{f}_{*})$
$\displaystyle=K_{**}-K_{*}^{\mathsf{T}}\left[K+V\right]^{-1}K_{*},$ (5)

with the matrices $K_{ij}=k(\bm{x}_{i},\bm{x}_{j})$,
$[K_{*}]_{ij}=k(\bm{x}_{i},\bm{x}_{*j})$,
$[K_{**}]_{ij}=k(\bm{x}_{*i},\bm{x}_{*j})$. Here $\bm{m}_{*}=m(X_{*})$.

This standard result from GP regression is used in the following section to
link the estimator for an unfolded distribution to a GP.

###  4.1 MAP estimator

We consider again the model from Sec. 3 with the likelihood
$P(\bm{n}|\bm{\mu})$ given by Eq. (3). From Bayes’ theorem, the log-posterior
probability is given by

$\log P(\bm{\mu}|\bm{n})=\log P(\bm{n}|\bm{\mu})+\log P(\bm{\mu})-\log
P(\bm{n}),$ (6)

where $P(\bm{\mu})$ is the prior probability, and the last term $P(\bm{n})$
may be ignored since it does not depend on $\bm{\mu}$.

We take the prior probability to be given by a GP with mean vector $\bm{m}$
(the _reference histogram_) and covariance matrix
$K_{ij}=k(\bm{x}_{i},\bm{x}_{j})$. The log-prior probability is then given by

$\log
P(\bm{\mu})=-\frac{1}{2}\left(\bm{\mu}-\bm{m}\right)^{\mathsf{T}}K^{-1}\left(\bm{\mu}-\bm{m}\right)+\ldots.$
(7)

Substituting the likelihood Eq. (3) and prior Eq. (7) into the expression for
the posterior Eq. (6), we obtain

$\displaystyle\log P(\bm{\mu}|\bm{n})=$
$\displaystyle-\frac{1}{2}\left(\bm{n}-R\bm{\mu}\right)^{\mathsf{T}}V^{-1}\left(\bm{n}-R\bm{\mu}\right)$
$\displaystyle-\frac{1}{2}\left(\bm{\mu}-\bm{m}\right)^{\mathsf{T}}K^{-1}\left(\bm{\mu}-\bm{m}\right)+\ldots,$
(8)

dropping terms which do not contain $\bm{\mu}$. The maximum a posteriori (MAP)
estimator $\hat{\bm{\mu}}$ is the mode of this posterior probability, and
maximises $\log P(\bm{\mu}|\bm{n})$. This summary statistic is found to be
given by

$\hat{\bm{\mu}}=K\left[K+R^{-1}V(R^{-1})^{\mathsf{T}}\right]^{-1}\left(R^{-1}\bm{n}-\bm{m}\right)+\bm{m}.$
(9)

A derivation of this is given in A.

By comparing the MAP estimator from Eq. (9) to that obtained from GP
regression in Eq. (4), we find the important result that $\hat{\bm{\mu}}$ is
the posterior mean of a GP regression whose observations are the ML solution,
which is given by $\hat{\bm{\mu}}_{\text{ML}}=R^{-1}\bm{n}$ with covariance
matrix $U_{\text{ML}}=R^{-1}V(R^{-1})^{\mathsf{T}}$. Since the posterior
distribution is a product of Gaussians, it is also Gaussian and therefore the
mode is identical to the mean. This connection allows us to write that the
covariance of the MAP estimator may be given by

$U=K-K\left[K+R^{-1}V(R^{-1})^{\mathsf{T}}\right]^{-1}K.$ (10)

Furthermore, if the observation (training) indices
$X=(\bm{x}_{1},\bm{x}_{2},\ldots)$ are different from the prediction (testing)
indices $X_{*}=(\bm{x}_{*1},\bm{x}_{*2},\ldots)$, and the reference histogram
can be obtained for bins defined by $X_{*}$, then we may use the standard
results from GP regression to generalise the MAP solution to

$\displaystyle\hat{\bm{\mu}}$
$\displaystyle=K_{*}^{\mathsf{T}}\left[K+U_{\text{ML}}\right]^{-1}\left(\hat{\bm{\mu}}_{\text{ML}}-\bm{m}\right)+\bm{m}_{*},$
(11) $\displaystyle U$
$\displaystyle=K_{**}-K_{*}^{\mathsf{T}}\left[K+U_{\text{ML}}\right]^{-1}K_{*},$
(12)

where $[K_{*}]_{ij}=k(\bm{x}_{i},\bm{x}_{*j})$,
$[K_{**}]_{ij}=k(\bm{x}_{*i},\bm{x}_{*j})$, and $\bm{m}_{*}$ is the mean
histogram at $X_{*}$.

The generalised results in Eq. (11) and Eq. (12) are simple algebraic
expressions once the ML solution is known. Therefore the unfolded estimator
and covariance are efficient to compute. This is an advantage over other, more
CPU-intensive unfolding schemes. In addition, these results are linear in
$\bm{n}$ so error propagation is simple.

###  4.2 Kernel choice and optimisation

In the proposed GP unfolding method, the regularisation is introduced via the
kernel function $k(\bm{x},\bm{x}^{\prime})$ which constricts the space of
possible solutions to those with a particular covariance. A common choice for
the kernel function is the squared-exponential:

$k(\bm{x},\bm{x}^{\prime})=A\,\exp\left(\frac{-\mathinner{\\!\left\lVert\bm{x}-\bm{x}^{\prime}\right\rVert}^{2}}{2l^{2}}\right).$
(13)

This kernel function is _stationary_ in the sense that it is a function of
only the distance between the inputs,
$\mathinner{\\!\left\lVert\bm{x}-\bm{x}^{\prime}\right\rVert}$. It is
parameterised by the amplitude $A$ and length scale $l$, referred to as the
set of _hyperparameters_ , $\bm{\theta}=\\{A,l\\}$. Various methods exist to
choose their values.

One method for this is by simulation, as is often done in particle physics
analyses. In this approach, a simulation program produces values for
$\bm{\mu}$, $\bm{\nu}$, and statistically independent $\bm{n}$. Then the
pseudo-data $\bm{n}$ are unfolded with varying hyperparameters to obtain
$\hat{\bm{\mu}}_{\bm{\theta}}$, and the agreement with $\bm{\mu}$ checked for
_closure_. The acceptable degree of closure, and its measure, are often chosen
by eye, although a more specific goodness-of-fit statistic may be used.

Another approach, taken from GP methods [8], is to maximise the marginal
likelihood,

$\displaystyle\log P(\bm{n};\bm{\theta})=$ $\displaystyle\log\left(\int
P(\bm{n}|\bm{\mu})\,P(\bm{\mu};\bm{\theta})\,\operatorname{d\\!}{\bm{\mu}}\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\left(\hat{\bm{\mu}}_{\text{ML}}-\bm{m}\right)^{\mathsf{T}}\left[K_{\bm{\theta}}+U_{\text{ML}}\right]^{-1}\left(\hat{\bm{\mu}}_{\text{ML}}-\bm{m}\right)$
$\displaystyle-\frac{1}{2}\log|{K_{\bm{\theta}}+U_{\text{ML}}}|-\frac{N}{2}\log{2\pi},$
(14)

where $K_{\bm{\theta}}\ =k_{\bm{\theta}}(X,X)$ is the kernel function
evaluated at $X$ with the hyperparameters set to $\bm{\theta}$. This is a
Bayesian approach, marginalising over the latent distribution $\bm{\mu}$. The
maximum of this marginal likelihood defines a model with a trade-off between
the fit to the data (the first term) versus model complexity (the second
term). For example, a GP using a squared-exponential kernel with very small
length scale $l$ will tend to fit to $\hat{\bm{\mu}}_{\text{ML}}$, but will be
overly complex (under-regularisation). In contrast, a large $l$ describes a
simpler model, but will fail to fit to the $\hat{\bm{\mu}}_{\text{ML}}$ (over-
regularisation). These extreme situations are penalised by the marginal
likelihood, whose maximum point may be used to choose values for
$\bm{\theta}$.

Finally, we mention the method of _cross validation_ , often employed for
hyperparameter optimisation in machine learning [16]. Various approaches for
cross validation exist, but given the relatively small number of bins $M$ in a
typical unfolding scenario, and the fast computation of the GP unfolded
result, we recommend the _leave-one-out_ variant. Here, $M$ sets of $M-1$
bins, $X_{*}^{(1)},X_{*}^{(2)},\ldots,X_{*}^{(M)}$, are produced with each
$X_{*}^{(i)}$ missing the $i$th bin. Then the prediction for $\mu_{i}$ is
compared to its true value via a _loss function_ , most often the squared
error. The set of $M$ losses for some hyperparameters $\bm{\theta}$ can then
be used to choose their values.

Other kernel functions may be more suitable for describing the truth
distribution. An attractive feature of the approach presented here is that one
may encode knowledge of the underlying physical process to derive a
physically-motivated kernel [17] which may better describe the truth
distribution.

The mathematics of _reproducing kernel Hilbert spaces_ formalises the link
between the kernel and the traditional regularisation approach used in some
particle physics results. For example, a _thin plate covariance_ [18] leads to
a solution equivalent to that of spline regularisation, known as Tikhonov
regularisation in particle physics [9, 10, 11, 19]. In one dimension, this
stationary kernel may be written $k(r)=A(2r^{3}-3Rr^{2}+R^{3})$, where
$r=|x-x^{\prime}|$ and $R$ is determined by boundary conditions. This kernel
contains a single parameter $A$, which controls the global strength of the
regularisation, as is the case with Tikhonov regularisation in its usual
implementation. In contrast, an advantage of the GP approach presented in this
paper is that the regularisation may be varied locally along the spectrum by
using a non-stationary kernel function. We provide an example of this in Sec.
5.2.

##  5 Example applications

Python code for the following examples may be found in Ref. [20]. We consider
the case of a bimodal distribution in Sec. 5.1, and a falling spectrum in Sec.
5.2.

###  5.1 Bimodal distribution

A set of 20 000 toy truth events is obtained by sampling from two Gaussian
distributions for $x$ with mean values 0.3 and 0.7, both with standard
deviation 0.1. These truth events are histogrammed in $\bm{\mu}$. The truth
events are then smeared with a Gaussian resolution of $\sigma=0.075$ to
generate the histogram $\bm{\nu}$. Events are accepted in the region $0<x<1$
and both the $\bm{\mu}$ and $\bm{\nu}$ histograms use 20 constant-width bins.
The truth and smeared events are used to determine the response matrix $R$
from a normalised 2D histogram. Finally, the observed histogram $\bm{n}$ is
generated by Poisson fluctuations around $\bm{\nu}$. The three histograms are
shown in Fig. 2.

Figure 1: Truth ($\bm{\mu}$), expected ($\bm{\nu}$), and observed ($\bm{n}$)
histograms for the two-peak unfolding example. The histogram definitions are
reported in the text. The error bars on $\bm{n}$ represent their Poisson
uncertainties.

Figure 2: Truth ($\bm{\mu}$) and unfolded truth estimators ($\hat{\bm{\mu}}$)
for the two-peak example. The error bars on $\hat{\bm{\mu}}$ represent the
standard deviations obtained from the covariance matrix as defined by Eq.
(10).

We use a GP with the squared-exponential kernel function given by Eq. (13).
The values for the two hyperparameters $A$ and $l$ are chosen to be those that
maximise the marginal likelihood given by Eq. (14). The maximum and contours
of the marginal likelihood are shown in Fig. 3.

Figure 3: Contours of the log-marginal likelihood Eq. (14) for the two-peak
example as a function of the parameters for the squared-exponential kernel,
$A$ and $l$. The cross indicates the point of maximum marginal likelihood. The
contour labels are the depth of the contour below the maximum.

The mean histogram for the GP, $\bm{m}$, is taken to be 0 for all bins since
it is found to have little impact on the final result. The unfolded estimator
for the truth, $\hat{\bm{\mu}}$ given by Eq. (9), is shown in Fig. 2. The
covariance matrix $U$ is defined by Eq. (10), and the correlation matrix
$\rho_{ij}=U_{ij}/\sqrt{U_{ii}U_{jj}}$ is shown in Fig. 4.

Figure 4: Correlation matrix for the unfolded truth estimators
$\hat{\bm{\mu}}$ for the two-peak example.

###  5.2 Falling spectrum

1 000 truth events are sampled from an exponential distribution $f(x)=e^{-x}$
in the region $1<x<5$ and accumulated in 20 bins of equal width ($\bm{\mu}$).
These events are smeared according to a Gaussian with resolution
$0.2\sqrt{x}$. The smeared events are placed in a histogram $\bm{\nu}$ with 30
bins of equal width in the region $0.5<x<5$, and the observed histogram
$\bm{n}$ is generated from a Poisson distribution around $\bm{\nu}$. These
three histograms are shown in Fig. 6.

Figure 5: Truth ($\bm{\mu}$), expected ($\bm{\nu}$), and observed ($\bm{n}$)
histograms for the falling spectrum example, as defined in the text. The error
bars on $\bm{n}$ represent their Poisson uncertainties.

Figure 6: Truth ($\bm{\mu}$) and unfolded truth estimators ($\hat{\bm{\mu}}$)
for the falling spectrum example.

In this example, $N>M$ so while the problem is well-constrained, the $N\times
M$ response matrix $R$ is not directly invertible. To mitigate this, we use
the ML estimator as the starting point and then regularise with the MAP
prescription detailed in Sec. 4. Specifically, we numerically maximise the
Gaussian likelihood in Eq. (3) using MINUIT [21] via the iminuit [22] Python
interface to obtain the ML estimator $\hat{\bm{\mu}}_{\text{ML}}$. The
covariance matrix for the ML estimator, $U_{\text{ML}}$ is obtained by
inverting the Hessian matrix with the HESSE subroutine. These results are then
substituted in Eq. (9) with
$R^{-1}\bm{n}\rightarrow\hat{\bm{\mu}}_{\text{ML}}$ and
$R^{-1}V(R^{-1})^{\mathsf{T}}\rightarrow U_{\text{ML}}$.

Since the (_a priori_ known) detector resolution increases proportionally with
$\sqrt{x}$, a kernel with constant length scale, such as the squared-
exponential Eq. (13) is unsuitable in this case. Therefore we choose a kernel
function with variable length scale, the Gibbs kernel [8, 23] in 1D,

$k(x,x^{\prime})=A\,\sqrt{\frac{2l(x)l(x^{\prime})}{l^{2}(x)+l^{2}(x^{\prime})}}\,\exp\left(-\frac{(x-x^{\prime})^{2}}{l^{2}(x)+l^{2}(x^{\prime})}\right),$
(15)

where $l(x)$ is an arbitrary positive function of $x$, here chosen to be
$l(x)=bx+c$. This allows for a linearly-changing length scale. The increased
flexibility afforded by this kernel function is realised by introducing more
regularisation parameters, $\bm{\theta}=\\{A,b,c\\}$. We remark that for a
large number of parameters, it becomes increasingly difficult to choose the
optimal point.

The unfolded estimators are shown against the truth histogram in Fig. 6. Here
the parameters $\bm{\theta}$ are chosen with the maximum marginal likelihood
prescription given in Sec. 4.2. As expected, $b>0$ so the length scale
increases with $x$.

##  6 Conclusions and Outlook

In this paper we have presented how GPs may be applied to the unfolding
problem. It is shown that conditioning a GP prior on the ML solution is
equivalent to constructing the MAP estimator. In this application, the use of
a GP regressor may be thought of as a method of regularising the ML solution.

The GP is entirely described by mean and kernel functions. While the mean
function is found to have little impact on the result, the kernel function
prescribes the covariance of unfolded estimation for the truth distribution.
By choosing an appropriate kernel function, the smoothness in the unfolded
estimator can be controlled. Furthermore, the kernel function has a direct
interpretation and may be motivated by knowledge of the underlying physics.
This means that, in contrast to other unfolding schemes, the regularisation is
a natural product of this approach.

For $N=M$, where the bins for the truth and observed histograms are equal, the
ML solution is simply given by inverting the response matrix,
$\hat{\bm{\mu}}_{\text{ML}}=R^{-1}\bm{n}$. However, generally $N\neq M$ and
this method may not be used. For the general case we envisage two
possibilities. First, the $N\times M$ response matrix $R$ is constructed from
Monte Carlo simulation, and the ML solution is found numerically by maximising
the likelihood given by Eq. (3). This is the approach taken in Sec. 5.2 for
the falling spectrum example. Alternatively, the square $N\times N$ matrix
$R^{\prime}$ could be constructed using the same binning for the truth and
observed histograms. Then the predictive GP mean function Eq. (11) is
evaluated at $X_{*}$, the $M$ centers of the bins for the desired truth
histograms.

In Sec. 2, we take the background contribution to be equal to zero for
simplicity. Background contributions, in the form of the $N$-dimensional
vector $\bm{\beta}$, are simple to include by modifying the folding equation
Eq. (1) to $\bm{\nu}=R\bm{\mu}+\bm{\beta}$. Then for the estimators used in
the method presented in this paper, one simply substitutes the data for the
background-subtracted data, $\bm{n}\rightarrow\bm{n}-\bm{\beta}$.

This paper assumes throughout that the data may be approximated as distributed
according to a Gaussian, and we note that this is not universally the case in
particle physics. However, the choice of unfolding method depends on the
analysis being done and should be tested against simulation in any case.
Therefore we recommend that for histograms with small bin populations, the
unfolding is tested to ensure it acceptably meets the requirements of the
analysis under consideration.

The treatment of systematic uncertainties is postponed to future work in this
area. We remark that that approximate variational approaches, as used in
published particle physics analyses [24, 25], may still be employed in this
case. We envisage further research into the applications of Student-$t$
processes [26] in unfolding in particle physics, as an extension to this work.

GPs have been introduced to a number of scientific fields to improve their
statistical procedures [3, 4]. They have not, however, traditionally been used
in particle physics, although recent developments in this area have shown
promise [17]. In this paper, we have introduced GPs into the important problem
of unfolding. We show that the method is generally applicable to problems of
different shapes and sizes, that the regularisation can be controlled
naturally, and that the result – including the unfolded covariance matrix –
can be obtained conveniently.

## Acknowledgements

We would like to extend our gratitude to our colleagues in ATLAS and RHUL for
their support. In particular, we thank Pedro Teixeira-Dias, and Lewis Wilkins
for their proofreading and helpful comments on this manuscript. We also thank
Veronique Boisvert and Pim Verschuuren for insightful discussions. This work
was supported by the UK Science and Technology Facilities Council.

##  Appendix A Derivation of MAP estimator

With reference to Eq. (8), we wish to find the value $\hat{\bm{\mu}}$ that
maximises the expression

$\displaystyle-\frac{1}{2}\left(\bm{n}-R\bm{\mu}\right)^{\mathsf{T}}V^{-1}\left(\bm{n}-R\bm{\mu}\right)-\frac{1}{2}\left(\bm{\mu}-\bm{m}\right)^{\mathsf{T}}K^{-1}\left(\bm{\mu}-\bm{m}\right).$
(16)

The derivative for each term is given by

$\displaystyle\frac{\partial}{\partial\bm{\mu}}\left[-\frac{1}{2}\left(\bm{n}-R\bm{\mu}\right)^{\mathsf{T}}V^{-1}\left(\bm{n}-R\bm{\mu}\right)\right]$
$\displaystyle=\left(\bm{n}-R\bm{\mu}\right)^{\mathsf{T}}V^{-1}R,$ (17)
$\displaystyle\frac{\partial}{\partial\bm{\mu}}\left[-\frac{1}{2}\left(\bm{\mu}-\bm{m}\right)^{\mathsf{T}}K^{-1}\left(\bm{\mu}-\bm{m}\right)\right]$
$\displaystyle=-\left(\bm{\mu}-\bm{m}\right)^{\mathsf{T}}K^{-1}.$ (18)

Combining these and taking the transpose ($V^{-1}$ and $K^{-1}$ are
symmetric), we therefore require that $\hat{\bm{\mu}}$ satisfies

$\displaystyle\bm{0}$
$\displaystyle=R^{\mathsf{T}}V^{-1}\left(\bm{n}-R\hat{\bm{\mu}}\right)-K^{-1}\left(\hat{\bm{\mu}}-\bm{m}\right)$
(19)
$\displaystyle=R^{\mathsf{T}}V^{-1}\bm{n}-\left[R^{\mathsf{T}}V^{-1}R+K^{-1}\right]\hat{\bm{\mu}}+K^{-1}\bm{m}.$
(20)

Now we use that the covariance of the ML solution from Sec. 3 is given by
$U_{\text{ML}}=R^{-1}V(R^{-1})^{\mathsf{T}}$ and therefore that
$R^{\mathsf{T}}V^{-1}R=U_{\text{ML}}^{-1}$. Substituting into Eq. (20) and
rearranging for $\hat{\bm{\mu}}$,

$\displaystyle\hat{\bm{\mu}}$
$\displaystyle=\left[K^{-1}+U_{\text{ML}}^{-1}\right]^{-1}\left(U_{\text{ML}}^{-1}R^{-1}\bm{n}+K^{-1}\bm{m}\right)$
(21)
$\displaystyle=K\left[K+U_{\text{ML}}\right]^{-1}R^{-1}\bm{n}+U_{\text{ML}}\left[K+U_{\text{ML}}\right]^{-1}\bm{m}$
(22)
$\displaystyle=K\left[K+U_{\text{ML}}\right]^{-1}\left(R^{-1}\bm{n}-\bm{m}\right)+\bm{m},$
(23)

where from Eq. (21) to Eq. (22) we use that
$\left[A^{-1}+B^{-1}\right]^{-1}B^{-1}\equiv A\left[A+B\right]^{-1}$ for
invertible matrices $A$ and $B$.

## References

  * [1] D. Foreman-Mackey, E. Agol, S. Ambikasaran, R. Angus, Fast and scalable Gaussian process modeling with applications to astronomical time series, Astron. J. 154 (6) (2017) 220.  arXiv:1703.09710, doi:10.3847/1538-3881/aa9332.
  * [2] I. Andersen, A. Szymkowiak, C. Rasmussen, L. Hanson, J. Marstrand, H. Larsson, L. Hansen, Perfusion quantification using Gaussian process deconvolution, Magn. Reson. Med. 48 (2) 351–361.  doi:10.1002/mrm.10213.
  * [3] J. A. Christen, B. Sansó, M. Santana-Cibrian, J. X. Velasco-Hernández, Bayesian deconvolution of oil well test data using Gaussian processes, J. Appl. Stat. 43 (4) (2016) 721–737.  doi:10.1080/02664763.2015.1077374.
  * [4] B. Hunt, Bayesian methods in nonlinear digital image restoration, IEEE Trans. Comput. C-26 (3) (1977) 219–229.  doi:10.1109/TC.1977.1674810.
  * [5] G. Cowan, A survey of unfolding methods for particle physics, in: Proc. Conference on Advanced Statistical Techniques in Particle Physics, Durham, England, 2002, pp. 248–257.
URL
http://www.ippp.dur.ac.uk/old/Workshops/02/statistics/proceedings/cowan.pdf

  * [6] V. Blobel, Unfolding methods in particle physics, in: Proc. PHYSTAT 2011 Workshop on Statistical Issues Related to Discovery Claims in Search Experiments and Unfolding, CERN, Geneva, Switzerland, 2011, pp. 240–251.  doi:10.5170/CERN-2011-006.240.
  * [7] F. Spanò, Unfolding in particle physics: a window on solving inverse problems, EPJ Web Conf. 55 (2013) 03002.  doi:10.1051/epjconf/20135503002.
  * [8] C. E. Rasmussen, C. K. I. Williams, Gaussian Processes for Machine Learning, The MIT Press, 2006.
  * [9] G. Cowan, Statistical Data Analysis, Oxford University Press, 1998.
  * [10] A. Tikhonov, On the solution of ill-posed problems and the method of regularization, Mat. Sb. 151 (3) (1963) 501–504.
  * [11] D. L. Phillips, A technique for the numerical solution of certain integral equations of the first kind, J. ACM 9 (1) (1962) 84–97.  doi:10.1145/321105.321114.
  * [12] M. Schmelling, The method of reduced cross-entropy: A general approach to unfold probability distributions, Nucl. Inst. Methods Phys. Res. A 340 (2) (1994) 400–412.  doi:10.1016/0168-9002(94)90119-8.
  * [13] G. Zech, Iterative unfolding with the Richardson-Lucy algorithm, Nucl. Inst. Methods Phys. Res. A 716 (2013) 1–9.  arXiv:1210.5177, doi:10.1016/j.nima.2013.03.026.
  * [14] G. D’Agostini, A multidimensional unfolding method based on Bayes’ theorem, Nucl. Inst. Methods Phys. Res. A 362 (1995) 487–498.  doi:10.1016/0168-9002(95)00274-X.
  * [15] G. Choudalakis, Fully Bayesian unfolding (2012).  arXiv:1201.4612.
  * [16] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine learning in Python, Journal of Machine Learning Research 12 (2011) 2825–2830.
  * [17] M. Frate, K. Cranmer, S. Kalia, A. Vandenberg-Rodes, D. Whiteson, Modeling smooth backgrounds and generic localized signals with Gaussian processes (2017).  arXiv:1709.05681.
  * [18] O. Williams, A. Fitzgibbon, Gaussian process implicit surfaces, in: Proc. Gaussian Processes in Practice, 2007, pp. 1–4.
URL https://www.microsoft.com/en-us/research/publication/gaussian-process-
implicit-surfaces-2/

  * [19] G. Wahba, Spline Models for Observational Data, Society for Industrial and Applied Mathematics, 1990.  doi:10.1137/1.9781611970128.
  * [20] https://github.com/adambozson/gp-unfold.
  * [21] F. James, M. Roos, Minuit – a system for function minimization and analysis of the parameter errors and correlations, Comput. Phys. Commun. 10 (1975) 343–367.  doi:10.1016/0010-4655(75)90039-9.
  * [22] iminuit – A Python interface to Minuit.
URL https://github.com/iminuit/iminuit

  * [23] M. N. Gibbs, Bayesian Gaussian processes for regression and classification, Ph.D. thesis, Univ. Cambridge (1997).
  * [24] G. Aad et al. (ATLAS Collaboration), Measurements of top-quark pair differential cross-sections in the lepton+jets channel in pp collisions at $\sqrt{s}=8\text{ TeV}$ using the ATLAS detector, Eur. Phys. J. C 76 (2016) 538\.  doi:10.1140/epjc/s10052-016-4366-4.
  * [25] G. Aad et al. (ATLAS Collaboration), Measurement of the differential cross-section of highly boosted top quarks as a function of their transverse momentum in $\sqrt{s}=8\text{ TeV}$ proton-proton collisions using the ATLAS detector, Phys. Rev. D 93 (2016) 1–34.  doi:10.1103/PhysRevD.93.032009.
  * [26] A. Shah, A. Wilson, Z. Ghahramani, Student-t processes as alternatives to Gaussian processes, in: Proc. Artificial Intelligence and Statistics, 2014, pp. 877–885.