# Implications of a matter-antimatter mass asymmetry in Penning-trap
experiments

Ting Cheng<EMAIL_ADDRESS>Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  Manfred Lindner
<EMAIL_ADDRESS>Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  Manibrata Sen
<EMAIL_ADDRESS>Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany

###### Abstract

The Standard Model (SM) of particle physics, being a local, unitary and
Lorentz-invariant quantum field theory, remains symmetric under the combined
action of Charge, Parity, and Time Reversal (CPT) symmetry. This automatically
implies that fundamental properties of particles and antiparticles should be
equal in magnitude. These fundamental tenets of the CPT principle have been
put to stringent tests in recent Penning-trap experiments, where the matter-
antimatter mass asymmetry has been measured. In light of these recent
advances, we compare the bounds arising on CPT invariance from kaon systems
with those from Penning-trap experiments. Using a simple yet powerful argument
of mass decomposition of hadrons, we show that bounds on quark-antiquark mass
differences from kaon oscillations are way beyond the reach of Penning-trap
experiments. We lay out a roadmap to discuss possible reformulations of our
understanding of the SM in the case of a discovery of CPT violation by these
precision experiments.

CPT violation, precision experiments

##  I Introduction

The Standard Model (SM) of Particle Physics is a local, Lorentz invariant,
Hermitian quantum field theory (QFT). As explored in a series of celebrated
papers Pauli (1940); Luders (1954); Jost (1957), one of the fundamental tenets
of such a local, Lorentz invariant theory is the conservation of CPT symmetry,
that is, invariance under the combined operations of charge conjugation,
parity inversion, and time reversal. The conservation of CPT guarantees that
physical properties of matter and antimatter are related, for example, their
masses should be identical, their charges, if any, should be equal and
opposite. In fact, the requirement that a Hermitian QFT is causal
automatically warrants the existence of antimatter which should have the exact
same mass of the corresponding matter field. Therefore, a test of whether
there exists a mass asymmetry between matter and antimatter, aptly dubbed here
as the matter-antimatter mass asymmetry (MAMA), automatically translates to a
test of the sacred principle of CPT invariance, and in turn, the foundations
of the SM.

Theoretically, a number of motivations exist for CPT symmetry to be exact,
relating the properties of matter and antimatter. However, the baryon
asymmetry of the Universe implies a matter dominated Universe. This indicates
that some form of asymmetry between matter and antimatter must have been
introduced through a new _yet-unknown_ mechanism in the early Universe. While
models of successful baryogenesis usually follow a CPT symmetric approach,
focusing on the Sakharov conditions Sakharov (1967), a baryon asymmetry could
also arise in thermal equilibrium in the presence of CPT violation, and baryon
number violation Bertolami  _et al._ (1997). Additionally, extensions of the
SM to incorporate a quantum theory of gravity often induces CPT violation
Tsujikawa (2013). Phenomenological motivations include the search for
violation of Lorentz invariance (LI), or violation of locality (L), leading to
CPT violation Colladay and Kostelecky (1997, 1998); Kostelecky and Lehnert
(2001); Coleman and Glashow (1997); Charlton  _et al._ (2020); Chaichian  _et
al._ (2011); Dolgov and Novikov (2012); Chaichian  _et al._ (2013a); Fujikawa
and Tureanu (2017).

Experimental tests of CPT invariance can be twofold: testing the properties of
particles and antiparticles directly, or probing the indirect impact of CPT
violation on other processes. In the absence of a specific model, it is
difficult to compare different experimental results at the same footing.
Therefore, in this _letter_ , we aim to bridge different mass measuring
experiments in a model independent manner, such that any new physics resulting
in CPT breaking can be investigated through a bottom-up approach. By treating
the constraints from each experimental result using dimensionless parameters,
the strongest constraint on CPT symmetry itself is currently from the kaon
oscillation experiments. The MAMA parameter is the difference between the two
diagonal terms of the Hamiltonian of $(K^{0},\bar{K}^{0})$ in flavor space,
and is tested using the Bell-Steinberger relation constructed under
assumptions of unitarity BAR (1984); Workman and Others (2022). Note that
neutral kaon-antikaon oscillations involve a process where _strangeness_ is
violated by 2 units at one loop, through a box-diagram (see Nierste (2009),
and references therein). In such a scenario, the test of MAMA could be more
sensitive to violations of the principle of locality, the underlying process
being a loop process. On the other hand, neutrino oscillation experiments also
provide another interesting test of CPT conservation. Here, for a given
neutrino energy, and baseline of the experiment, oscillation parameters (mass-
squared difference, mixing angles) are fitted separately for the neutrino and
antineutrino spectra Colladay and Kostelecky (1997, 1998); De Gouvea (2002);
Bahcall  _et al._ (2002); Kostelecky and Mewes (2004a, b); Datta  _et al._
(2004); Minakata and Uchinami (2005); Kostelecky and Mewes (2012); Ohlsson and
Zhou (2015); Abe  _et al._ (2015); Diaz and Schwetz (2016); de Gouvêa and
Kelly (2017); Barenboim  _et al._ (2018). In this case, MAMA is measured with
respect to the dispersion relation of the propagating neutrino, and therefore,
provides a more sensitive probe of Lorentz-invariance violation Liberati
(2013). For a discussion regarding tests of non-locality using neutrino
oscillation experiments, see Antusch and Fernandez-Martinez (2008) and
references therein. In all these cases, tests of CPT conservation usually
quote the results in terms of the mass difference between the particles and
the antiparticles. It is important to emphasize that while the definition of
the measured “mass” may be different from experiment to experiment, these
systems can eventually be sensitive to multiple underlying principles behind
the breaking of CPT: non-locality, and/or Lorentz invariance violation (LI-V).
On the other hand, not all CPT breaking effects are covered in our analysis,
since we have narrowed down the arbitrariness of the origin of CPT violation
by considering only MAMA tests.

Significant progress has also been achieved by precision experiments in
testing the tenets of the CPT principle. The ALPHA experiment at CERN uses
trapped antihydrogen to study its charge-neutrality Amole  _et al._ (2014);
Ahmadi  _et al._ (2016), the ratio of gravitational mass to inertial mass
Bertsche (2018), as well as a measurement for the hyperfine splitting in
neutral antihydrogen Ahmadi  _et al._ (2017). The experiment has, hence far,
shown that the SM is consistent with CPT conservation. Similar tests have been
performed using other species of antimatter such as antiprotons, which can be
trapped for longer times in Penning-traps Brown and Gabrielse (1986). The BASE
collaboration Borchert  _et al._ (2022) at CERN measures the charge-to-mass
ratio ($q/m$) of the proton/antiproton by comparing the cyclotron frequency
($\nu_{c}$) of a single antiproton $\bar{p}$ to those of a single negatively
charged hydrogen $H^{-}$(to avoid systemic uncertainties by having the charge
different from the antiproton). This is done using the Brown-Gabrielse
invariance theorem: $\nu_{c}^{2}=\nu_{+}^{2}+\nu_{z}^{2}+\nu_{-}^{2}$, where
$\nu_{+},\nu_{z},\nu_{-}$ are three eigenfrequencies, namely, the modified
cyclotron frequency, the axial frequency and the magnetron frequency, and
$\nu_{c}=1/(2\pi)(q/m)B_{0}$. The Penning-trap captures a $H^{-}$ and a
$\bar{p}$ produced at CERN; then by measuring the three eigenfrequencies under
a homogeneous magnetic field of $B_{0}=1.945\,{\rm T}$, one can obtain the
ratio of the proton to antiproton’s inertial mass. Since the inertial mass is
directly measured in this setup, measurement of the annual modulation of the
MAMA would relate directly to tests of extension of the weak equivalence
principle, such as having scalar-tensor theories of modified gravity Hughes
and Holzscheiter (1991); Charlton  _et al._ (2020).

With technology advancing in leaps and bounds, the sensitivities of the
Penning-trap experiments are expected to get even better with time. In light
of these facts, it is important to compare the bounds on CPT arising from kaon
systems, and neutrino systems, with those from current and upcoming precision
experiments. In this _letter_ , we present a simple yet crucial argument
relying on the mass decomposition of hadrons using the energy momentum tensor
in Quantum Chromodynamics (QCD), which allows a hadron mass to be separable
into individual quark contributions, and those coming from gluons, kinetic
terms, as well as anomaly terms. To zeroth order, this allows the MAMA in the
hadron system to be written as the mass difference between quark and
antiquarks. Using this parameterization, the existing bounds from kaon systems
translates to bounds on the mass difference between quarks and antiquarks,
which are well beyond the sensitivity of Penning-trap experiments to such
quark-antiquark mass differences. Similar outcomes with regards to the
differences in gravitational forces exerted on matter and antimatter have been
presented in Caldwell and Dvali (2021). This emphasizes that any discovery by
these precision experiments would warrant a serious reformulation of our
understanding of QFTs in order to be compatible with the results from the kaon
systems.

Our discussion is organized as follows. In the next section, we outline the
basic framework of decomposition of a hadron mass in QCD, and how it relates
to the constituent quark masses. Then, we apply it to different systems:
protons, kaons, and discuss the bounds arising on MAMA. We also compare bounds
on CPT arising from neutrino oscillations. We finally discuss the theoretical
implications of a positive measurement of MAMA from the Penning-trap
experiments, and/or other experiments, and conclude.

##  II Hadron mass decomposition in QCD

The masses of hadrons in QCD are governed primarily by the interaction amongst
the constituent quarks and gluons. In fact, while hadron masses have been
quite well measured experimentally, disentangling the mass contributions of
the individual quarks from the quark-gluon strong interaction dynamics is
still a matter of active research. The incredible complexity of the problem
arises from the fact that at low energies, QCD is non-perturbative, and hence
one needs to resort to lattice computations to get a clearer understanding.

However, it is possible to obtain a phenomenological decomposition of the mass
of a hadron using the energy-momentum tensor $(T_{\mu\nu})$ in QCD Ji (1995);
Yang  _et al._ (2015). The symmetric, conserved energy-momentum tensor in QCD
can be formally written in Euclidean space as

$T_{\mu\nu}=\frac{1}{4}\bar{\psi}\gamma_{(\mu}\overleftrightarrow{D}_{\nu)}\psi+F_{\mu\alpha}F_{\nu\alpha}-\frac{1}{4}\delta_{\mu\nu}F^{2},$
(1)

where $\overrightarrow{D}=\overrightarrow{\partial}_{\mu}+igA_{\mu}$ and
$\overleftarrow{D}=\overleftarrow{\partial}_{\mu}-igA_{\mu}$, for a gluon
field $A_{\mu}$ and coupling $g$, while the $()$ in the kinetic term denotes
symmetrization over the indices. The QCD Hamiltonian operator can be defined
through $T_{\mu\nu}$ as

$H_{\rm QCD}=-\int d^{3}x\,T_{44}(x)~{}.$ (2)

This Hamiltonian can be decomposed into four contributions, coming from the
kinetic energy of the quarks, the gluon field energy, the bare quark masses,
and the QCD anomaly respectively, given by

$H_{\rm QCD}=H_{E}+H_{g}+H_{m}+H_{a}~{},$ (3)

where

$\displaystyle H_{E}$ $\displaystyle=$ $\displaystyle\sum_{q}\int
d^{3}x\,\bar{\psi}_{q}(\overrightarrow{D}.\overrightarrow{\gamma})\psi_{q}~{},$
(4) $\displaystyle H_{g}$ $\displaystyle=$ $\displaystyle\int
d^{3}x\,\frac{1}{2}(B^{2}-E^{2})~{},$ (5) $\displaystyle H_{m}$
$\displaystyle=$ $\displaystyle\sum_{q}\int
d^{3}x\,m_{q}\bar{\psi}_{q}\psi_{q}~{},$ (6) $\displaystyle H_{a}$
$\displaystyle=$ $\displaystyle\,\int
d^{3}x\,\left[\frac{\gamma_{m}}{4}\sum_{q}m_{q}\bar{\psi}_{q}\psi_{q}-\frac{\beta(g)}{4g}(B^{2}+E^{2})\right].\,\,\,$
(7)

Here $\gamma_{m}$ is the anomalous mass dimension operator, and $\beta(g)$ is
the QCD beta function. Using this decomposition, the masses of hadrons were
estimated using lattice QCD for the first time to identify the contributions
of different components Yang  _et al._ (2015, 2018).

The rest mass of a single hadron state $|x\rangle$ is defined as

$m_{x}=\frac{\langle x|H_{\rm QCD}|x\rangle}{\langle x|x\rangle}~{}.$ (8)

We stress that by focusing on such a mass measurement, we are narrowing down
the arbitrariness of the origin of CPT violation to create a rather model
independent bridge relating different mass testing experiments. In other
words, $m_{x}$ does not include all CPT violating effects, for instance, terms
that do not appear in the rest frame are not included. For example, a CPT
transformation would transform a proton state ($|p\rangle$) onto an antiproton
state ($|\bar{p}\rangle$), and vice versa. Hence, a non-zero mass difference
between them would indicate CPT viloation. However, since the charge of the
proton is decided solely by the valence quark charges, one can define a CqPT
transformation with the charge conjugation operator only acting quark fields
and states (parity and time reversal operators still remain universal), which
would do the same trick. Therefore, terms independent of quark fields in the
mass decomposition formula, such as $H_{g}$, would only experience a PT
transformation under CqPT, which we take as invariant in the rest frame. This
is the fundamental assumption that enters our calculation.

From the mass decomposition formula, it is clear that the masses of the quarks
contribute to the hadron mass dominantly through the bare quark mass term in
Eq. 6, and the anomaly term in Eq. 7, the remaining contributions being from
the gluon field energy. The gluon energy contributions depend only on the
gluon field, and does not depend dominantly on the quarks involved. As a
result, when considering the mass difference between a hadron $x$, and its
corresponding antihadron ${\bar{x}}$, the gluon contributions can be expected
to cancel out at zeroth order. This mass difference should be dominated by
$H_{m}$, followed by subdominant contributions from $H_{a}$ and $H_{E}$. Note
that the above formalism is derived under the assumption of a local and LI
QFT; any violation from this might contribute to the quantum corrections in
$H_{a}$. However, the arguments that follow will hold as these can be absorbed
by a rescaling of the coefficients of $H_{g}$ and $H_{m}$. Namely, having the
same field structure, we can combine the $\gamma$ and $\beta$ terms in $H_{a}$
with $H_{m}$ and $H_{g}$, respectively, and redefine the phenomenological
coefficients. For instance, we can redefine $m_{q}\rightarrow
m_{q}(1+\gamma_{m})$, where $|\gamma_{m}|<<1$ Baikov  _et al._ (2014). Another
equivalent way of stating the same is the following: a hadron mass can be
assumed to be due to a combination of its bare quark masses, the gluon self-
energy, and the quantum corrections received by these terms. The gluon self-
energy term, as well as the quantum corrections are, to zeroth order,
independent of the different constituent quarks of the hadron. As a result,
when one takes the difference between the mass of a hadron and its
corresponding antihadron, the difference due to the bare quark masses are the
only terms to survive. Using the results of tests of CPT invariance, the
hadron mass decomposition can immediately be used to put very stringent
constraints on the mass difference between a quark and an antiquark. This is
the underlying principle of our _letter_. We elaborate on this in the
following section.

##  III Comparison among the proton, the kaon and the neutrino sector

The BASE collaboration Borchert  _et al._ (2022) measures the effective
charge-to-mass ratio of protons and antiprotons to obtain limits on the mass
asymmetry from their result $(q/m)_{p}/(q/m)_{\bar{p}}=-1.000000000003(16)$,

$\left|\frac{m_{\bar{p}}}{m_{p}}-1\right|<3\times 10^{-12}\,.$ (9)

Following our argument in the previous section, assuming the charges are equal
and opposite for particles and antiparticles, this bound on mass difference
can be translated onto the mass differences in the constituent valence quarks
and antiquarks. For estimation purpose, we consider three parameterizations of
the MAMA in quarks. We define the difference in quark-antiquark masses as
$\delta_{q}\equiv m_{\bar{q}}-m_{q}$, and the ratio between the quark-
antiquark masses as $r_{q}\equiv m_{\bar{q}}/m_{q}$, where $q=s,d,u$ denote
valence quarks. Here, the mass terms $m_{q}$ and $m_{\bar{q}}$ represent the
pole energies of the propagator in the rest frame, which could be shifted
differently for different spinor eigenstates. In particular, from Eq. (8), we
can see that $m_{q}$ and $m_{\bar{q}}$ represent the energy eigenvalues for
the corresponding spinor state of $x$ in the rest frame.

We can also parameterize the CPT violating term through, $\alpha$, where
$\alpha$ is defined as $m_{x}=m_{0}(1+\alpha)$ for particles and
$m_{\bar{x}}=m_{0}(1-\alpha)$, for antiparticles. In that case,

$\alpha\equiv\left|\frac{m_{\bar{x}}-m_{x}}{m_{\bar{x}}+m_{x}}\right|\simeq\left|\frac{\sum_{j}\delta_{j}}{2m_{x}}\right|\,.$
(10)

The limits of the Penning-trap system can thus be translated into
$|2\delta_{u}+\delta_{d}|$ and $r-1$ assuming $r=r_{d}=r_{u}$. The left-hand-
side of Eq. 9 can therefore be rewritten as
$\left|\sum_{q}C_{q}\delta_{q}/m_{p}\right|$ or
$\left|(r-1)\sum_{q}C_{q}m_{q}/m_{p}\right|$, where $C_{q}=\langle
P|\bar{\psi}_{q}\psi_{q}|P\rangle=\partial m_{p}/\partial m_{q}$, and would be
the same for proton and antiproton. Furthermore, since the sea quarks would
also be identical between proton and antiproton, CPT breaking should be from
the valance quarks with contributions proportional to the respective bare
quark charges, namely, $C_{u}:C_{d}=2:1$. The absolute value of $C_{q}$ can be
estimated by considering an analogy with the mass splitting between proton and
neutron ($\delta m^{\rm QCD}_{\rm N}$) induced by chiral symmetry breaking.
Namely, the isovector scalar charge $g^{u-d}_{s}\sim 1$ Aoki  _et al._ (2022)
for $g^{u-d}_{s}=\delta m^{\rm QCD}_{\rm N}/\delta m_{ud}$, where $\delta
m_{ud}$ denotes the up quark-down quark mass splitting. On the other hand,
$C_{u}$ and $C_{d}$ can also be estimated through the pion-nucleon
$\sigma$-term, $\sigma_{\pi N}=(m_{d}+m_{u})/2\langle
N|\bar{\psi_{u}}\psi_{u}+\bar{\psi_{d}}\psi_{d}|N\rangle$, from pion-nucleon
scattering experiments or through lattice calculations, see Aoki  _et al._
(2022); Gupta  _et al._ (2021) for a review. With the purpose of order-of-
magnitude comparison between different systems, Table 1 is listed
conservatively to a range such that inputs from the chiral symmetry breaking
analogy, the $\sigma$-term and $\langle H_{m}\rangle/m_{p}$ from the mass
decomposition Yang  _et al._ (2018), can all be included. Nonetheless,
$\alpha$ is directly related to the experimental result, hence it would be
independent of non-perturbative estimations, and is set as $\alpha<1.5\times
10^{-12}$ by Eq. 9.

MAMA | Proton | Kaon | Neutrino
---|---|---|---
$|\sum_{j}\delta_{j}|$ (MeV) | $\mathcal{O}(10^{-10}-10^{-9})$ | $\mathcal{O}(10^{-16})$ | $\mathcal{O}(10^{-9})$
$\delta$ (MeV) | $\mathcal{O}(10^{-10}-10^{-9})$ | trivial | $\mathcal{O}(10^{-9})$
$r-1$ | $\mathcal{O}(10^{-11}-10^{-10})$ | $\mathcal{O}(10^{-18})$ | $\mathcal{O}(10^{-1})$
$\alpha$ | $\mathcal{O}(10^{-12})$ | $\mathcal{O}(10^{-19})$ | $\mathcal{O}(10^{-2})$
Table 1: Theoretical limits on the different parameters of CPT violation.
$|\sum_{j}\delta_{j}|$ denotes the sensitivity of MAMA for each experiment,
which is $|2\delta_{u}+\delta_{d}|$, $|\delta_{s}-\delta_{d}|$ and
$(|\delta_{2}|,|\delta_{3}|)$ from the Penning-trap, neutral kaon oscillation
and neutrino oscillation experiments, respectively. $\delta$ assumes all
$\delta_{j}$s are identical as such value. The other quantities are as defined
in the text.

The neutral kaon oscillation system offers a stronger limit on such MAMA. The
Bell-Steinberger relation Bell and Steinberger (1966) allows a direct
connection between the Hamiltoninan eigenstates, and the decay amplitudes in
the $K^{0}-\bar{K}^{0}$ system. From this, the limits on the mass difference
between the kaon and the antikaon can be expressed using the diagonal elements
of the Hamiltonian as Workman and Others (2022)

$|m_{K^{0}}-m_{\bar{K}^{0}}|<4\times 10^{-16}\,{\rm MeV}\,.$ (11)

Following the same line of reasoning as that for the proton system, the
entries for the kaon system in Table 1 can be estimated for only bare quark
contributions ($C=1$), or all valence quark contribution ($C\sim 2.1$
considering Aguilar  _et al._ (2019)) or even including sea quarks ($C\sim
2.5$ considering Aguilar  _et al._ (2019); Yang  _et al._ (2015)). Here, we
assume $C=C_{s}=C_{d}$. Note that Eq. 11 would be trivially satisfied when
$\delta_{s}=\delta_{d}$. A rather engineered way to escape stringent bounds
set by the neutral kaon oscillation experiments is to have an additional
symmetry which sets $\delta_{s}$ (nearly) identical to $\delta_{d}$. On the
other hand, $\alpha$ would be independent of $C_{q}$ and would be limited to
$\alpha<4.0\times 10^{-19}$. Clearly, the bounds arising from neutral kaon
oscillation systems are roughly seven orders of magnitude stronger than that
of the Penning-trap experiment. This is due to a combination of the lower mass
of the kaon, as well as the stronger experimental bounds existing on this
system from tests of CPT conservation. Note that the above results should not
be considered as precise calculations, but more as estimates. What is
important to note, however, is that the naive bounds from kaon systems are
_orders of magnitude_ above proton systems, and hence clearly set a much more
stringent bound on MAMA than what can be possibly achieved by such Penning-
trap experiments in the foreseeable future.

While the discussion has been centered around hadronic systems so far, it is
also worthwhile to consider the impact of CPT invariance on MAMA for the
neutrino sector. The main advantage in considering neutrinos is that they are
fundamental particles, and hence are not plagued by the uncertainties existing
in QCD. Furthermore, neutrino masses are sensitive to new physics at higher
mass scales, where it is possible to imagine that the underlying theory might
be non-local and/or LI-V. Neutrino oscillation experiments can test CPT
invariance by its effect on the mass-squared differences Colladay and
Kostelecky (1997, 1998); De Gouvea (2002); Bahcall  _et al._ (2002);
Kostelecky and Mewes (2004a, b); Datta  _et al._ (2004); Minakata and Uchinami
(2005); Kostelecky and Mewes (2012); Ohlsson and Zhou (2015); Abe  _et al._
(2015); Diaz and Schwetz (2016); de Gouvêa and Kelly (2017); Barenboim  _et
al._ (2018). A dedicated analysis, using a combination of solar neutrino data,
and KamLAND reactor antineutrino data gives Barenboim  _et al._ (2018);
Tórtola  _et al._ (2020)

$\Delta m^{2}_{21}-\Delta\bar{m}^{2}_{21}<4.7\times 10^{-5}\,{\rm eV}^{2}\,,$
(12)

while long-baseline, and short-baseline experiments can set

$\Delta m^{2}_{31}-\Delta\bar{m}^{2}_{31}<2.5\times 10^{-4}\,{\rm eV}^{2}\,.$
(13)

For a given value of the lightest neutrino mass from zero to the upper limit
from cosmological bounds, this translates to bounds for MAMA parameters of the
two heavier neutrinos shown in Table 1. The order-of-magnitude estimations
hold for both neutrinos.

Note that in Barenboim  _et al._ (2018), the bounds arising on CPT violation
from Eq. (12-13) were compared to the quantity
$m^{2}_{K^{0}}-m^{2}_{\bar{K}^{0}}$, which is the quantity in the Lagrangian
that determines the oscillation in the kaon systems Murayama (2004). As a
result, the authors of Barenboim  _et al._ (2018) obtained that bounds on CPT
violation in the neutrino sector are stronger than the kaon sector. On the
other hand, we have used the parameterizations in Table 1, which rely on the
mass-difference and not the mass-difference-squared, as a fundamental
parameter that CPT violation tests. As a result, we obtain that the bounds on
CPT from neutrino system , while stringent, are not competitive to the kaon
system to constrain MAMA. Of course, it is crucial to remember that these two
bounds are fundamentally different: while the kaon is a composite particle in
QCD, and the limits quoted in Eq. 11 offer a probe of QCD dynamics, the
neutrinos are elementary fermions and hence can offer a clearer test of the
fundamental symmetry.

##  IV Theoretical Implications of CPT violation

This brings us to an important question. In the SM and in any extension based
on QFT, as long as quarks and neutrinos share the same source of CPT
violation, we do not expect to see a signal of MAMA from either the Penning-
trap experiment nor from neutrino oscillation due to the upper bound set by
kaon oscillation. However, if the current Penning-trap experiments would
observe a mass difference between protons and antiprotons, or in general, if a
nontrivial MAMA is observed in any system, what kind of new physics could this
signify? Immediately, this would indicate that CPT symmetry is broken, as
shown in Fig. 1, most likely by some new physics at high scale. The CPT
theorem states that in any _local, Lorentz invariant_ field theory, CPT is an
exact symmetry. Clearly, if any of the above conditions are violated, CPT need
not be an exact symmetry. In other words, unlike CP violation which is treated
independently in different sectors, CPT violation stems from the fundamental
principles not being exact, which is, in general, universal among different
sectors. For instance, if Lorentz invariance fails in the quark sector, there
is no obvious reason as to why it should not fail as well in the lepton
sector. Furthermore, as shown previously, one can connect MAMA in the hadronic
system to that of the fundamental particles under the premise that a subset of
the CPT transformation, $C_{q}PT$, is sufficient to reflect the hadronic state
onto its antihadronic state. In this case, since only the rest mass is under
consideration, the model independent mass decomposition formalism indicate
that other terms would either be cancelled out or absorbed into
$m_{q}/m_{\bar{q}}$ for the MAMA of hadrons. Here, we also don’t consider spin
dependence of the MAMA since we work in the rest frame, nonetheless, one can
in principle do a spin dependent measurement to check this assumption.

In this section, we outline some of the fascinating new directions which can
be explored in this context. This is intended to serve as a scaffolding for
the establishment of new physics. Our arguments can be summarized in Fig. 1. A
discovery of $m\neq\bar{m}$ would certainly imply violation of Lorentz
invariance and/or locality violation. This necessitates an extensions of the
SM, incorporating one or both of these tenets. Note that such extensions
usually violate micro-causality, and hence must be embedded in a UV complete
theory. However, for theories from spontaneous LI violation Kostelecky and
Lehnert (2001) or non-local field theories Tomboulis (2015), causality would
still be stably conserved at low energies, since the violation of micro-
causality is confined in the high energy region. This might be a way out in
case of a positive discovery without restricting LI-V and L-V contributions
for the breaking of micro-causality to cancel out, as shown in Fig. 2.

Figure 1: Flow chart for implications of a non-trivial asymmetry between the
mass of matter and antimatter ($m\neq\bar{m}$). LI/L/CPT-V(C) means that
Lorentz invariance/locality/CPT symmetry is violated (conserved).  Figure 2:
Both Lorentz invariant violation (LI-V) and locality violation (L-V) will
break micro-causality, however, having both might cancel out such breaking
Kostelecky and Lehnert (2001); Tomboulis (2015)

The Standard Model Extension (SME) Colladay and Kostelecky (1997, 1998) aims
to do this by writing down a phenomenological set of Lorentz-invariance
violating operators. Lorentz invariance can be broken either _spontaneously_
Kostelecky and Lehnert (2001), where some fundamental field in a Lorentz
invariant theory takes a vacuum expectation value (vev), and breaks the
symmetry, or _explicitly_ Coleman and Glashow (1997), where the Lagrangian
contains operators which break this symmetry from the beginning Charlton  _et
al._ (2020). In many cases, note that these vanilla scenarios predict the same
charge-to-mass ratio for particles and antiparticles; a violation of this will
make the theory non-causal Greenberg (2002). However, what matters is whether
the low energy observables are causal, without going into the details of the
high energy theory.

To demonstrate how one can link a certain model to our MAMA parameter, we take
the minimum SME (mSME) as an example. The CPT symmetry breaking terms in the
mSME Lagrangian with a $H_{m}$ field structure, namely $a_{\mu}$ and
$b_{\mu}$, satisfy the Dirac equation Colladay and Kostelecky (1997)

$(i\gamma_{\mu}D_{\mu}-m-a_{\mu}\gamma_{\mu}-b_{\mu}\gamma_{5}\gamma_{\mu})~{}\psi_{q}=0.$
(14)

The field operator can be written as

$\psi_{q}(x)=\int\frac{d^{3}\mathbb{p}}{(2\pi)^{3}}\sum_{s}\left(\frac{c^{s}_{\mathbb{p}}u_{q}^{s}(\mathbb{p})e^{-\mathrm{i}p_{u}^{s}x}}{\sqrt{2E_{u}^{s}}}+\frac{d^{s\dagger}_{\mathbb{p}}v_{q}^{s}(\mathbb{p})e^{\mathrm{i}p_{v}^{s}x}}{\sqrt{2E_{v}^{s}}}\right),\\\
$ (15)

where $c_{\mathbb{p}}^{s},d_{\mathbb{p}}^{s\dagger}$ denote corresponding
creation and annihilation operators, while $u^{s}_{q},v^{s}_{q}$ are the
spinors. The four-momenta are $p_{u}^{s}=(E_{u}^{s},\mathbb{p})$,
$p_{v}^{s}=(E_{v}^{s},\mathbb{p})$, with $E_{u}^{s}$ and $E_{v}^{s}$ depending
on $a_{\mu}$ and $b_{\mu}$ to satisfy the Dirac equation. For simplicity, we
assume only $a_{\mu}$ is non-zero, then
$E_{u}^{s}=E_{u}=\sqrt{m^{2}+(\mathbb{p}-\mathbb{a})^{2}}+a_{0}$ and
$E_{v}^{s}=E_{v}=\sqrt{m^{2}+(\mathbb{p}+\mathbb{a})^{2}}-a_{0}$. Hence, in
the rest frame, $m_{q}=E_{u}(\mathbb{p}=0)$, $m_{\bar{q}}=E_{v}(\mathbb{p}=0)$
and the MAMA parameter $\delta_{q}=2\,(a_{q})_{0}$. From Table 1, we get
$(a_{s})_{0}-(a_{d})_{0}<2\times 10^{-16}$ MeV from kaon oscillation and
$2(a_{u})_{0}+(a_{d})_{0}<2\times 10^{-10}$ MeV from the Penning trap bounds
after taking all $C_{q}$s as one. The energy eigenvalues for the case where
both $a_{u}$ and $b_{u}$ are non-zero are given in Colladay and Kostelecky
(1997), which would in general be spin dependent, giving four eigenvalues
corresponding to the four spinors. However, in the rest frame, two of them
would be degenerate with the other two, leaving the MAMA parameter still well
defined.

The other alternative is to extend the Lagrangian such that it becomes a _non-
local_ theory. For e.g., Chaichian  _et al._ (2011); Dolgov and Novikov
(2012); Chaichian  _et al._ (2013a); Fujikawa and Tureanu (2017) considered a
Lorentz invariant extension of QED, which however, has non-local operators,
thereby leading to CPT-violation. These kind of theories predict different
masses for particles and antiparticles at scales much below the non-locality
scales, and hence can have direct implications for Penning-trap like
experiments. However, these theories also suffer from non-unitarity Chaichian
_et al._ (2013b), as well as micro-causality problems, and should be viewed
strictly as an effective theory. Other possibilities which might lead to such
mass difference can arise in composite quark theories Larin (2020), where
different quarks pick up different CPT phases, or in models of extra dimension
like string theory Kostelecky and Samuel (1989a, b), or in theories of quantum
gravity Wald (1980).

##  V Conclusion

The conservation of CPT symmetry is a sacrosanct feature of any local, Lorentz
invariant theory. Experiments are underway to test this theory by probing the
charge-to-mass ratio of protons-antiprotons. The BASE collaboration at CERN
has recently measured the proton/antiproton charge-to-mass ratio to an
accuracy of parts per trillion. The collaboration plans to improve their
sensitivity further by upgrading the stability of the magnetic field, as well
as using transportable antiproton traps. In this regard, it is crucial to
examine the bounds arising on similar measurements from other experiments.

The most stringent bound on the extent of CPT violation comes from kaon-
antikaon oscillation experiments. Using a simple phenomenological modelling of
hadron masses, we have shown that the kaon bounds translate to bounds on the
CPT violation parameter that are orders of magnitude beyond the limits from
Penning-trap experiments. These bounds even supersede those from neutrino
oscillation experiments. As a result, any positive results from Penning-trap
experiments would definitely be a smoking-gun signal for fundamental new
physics. This can either point towards a possible violation of CPT, thereby
leading to a reformulation of the basic tenets of QFT, or a possible lack of
understanding in the QCD sector. Additionally, if the positive result appears
in the neutrino sector, it would point toward theories where the underlying
source of CPT violation is different for neutrinos than for quarks. We have
laid a tentative road-map for the future in the event of such a discovery.

## Acknowledgements

We would like to thank IMPRS-PTFS for organizing the retreat which led to
discussions shaping this paper and the support for TC.

## References

  * Pauli (1940) W. Pauli, Phys. Rev. 58, 716 (1940).
  * Luders (1954) G. Luders, Kong. Dan. Vid. Sel. Mat. Fys. Med. 28N5, 1 (1954).
  * Jost (1957) R. Jost, Helv. Phys. Acta 30, 409 (1957).
  * Sakharov (1967) A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5, 32 (1967).
  * Bertolami  _et al._ (1997) O. Bertolami, D. Colladay, V. A. Kostelecky, and R. Potting, Phys. Lett. B 395, 178 (1997), arXiv:hep-ph/9612437 .
  * Tsujikawa (2013) S. Tsujikawa, Class. Quant. Grav. 30, 214003 (2013), arXiv:1304.1961 [gr-qc] .
  * Colladay and Kostelecky (1997) D. Colladay and V. A. Kostelecky, Phys. Rev. D 55, 6760 (1997), arXiv:hep-ph/9703464 .
  * Colladay and Kostelecky (1998) D. Colladay and V. A. Kostelecky, Phys. Rev. D 58, 116002 (1998), arXiv:hep-ph/9809521 .
  * Kostelecky and Lehnert (2001) V. A. Kostelecky and R. Lehnert, Phys. Rev. D 63, 065008 (2001), arXiv:hep-th/0012060 .
  * Coleman and Glashow (1997) S. R. Coleman and S. L. Glashow, Phys. Lett. B 405, 249 (1997), arXiv:hep-ph/9703240 .
  * Charlton  _et al._ (2020) M. Charlton, S. Eriksson, and G. M. Shore, (2020), arXiv:2002.09348 [hep-ph] .
  * Chaichian  _et al._ (2011) M. Chaichian, A. D. Dolgov, V. A. Novikov, and A. Tureanu, Phys. Lett. B 699, 177 (2011), arXiv:1103.0168 [hep-th] .
  * Dolgov and Novikov (2012) A. D. Dolgov and V. A. Novikov, JETP Lett. 95, 594 (2012), arXiv:1204.5612 [hep-ph] .
  * Chaichian  _et al._ (2013a) M. Chaichian, K. Fujikawa, and A. Tureanu, Phys. Lett. B 718, 1500 (2013a), arXiv:1210.0208 [hep-th] .
  * Fujikawa and Tureanu (2017) K. Fujikawa and A. Tureanu, Int. J. Mod. Phys. A 32, 1741014 (2017), arXiv:1607.01409 [hep-ph] .
  * BAR (1984) Nuclear Physics B 247, 293 (1984).
  * Workman and Others (2022) R. L. Workman and Others (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Nierste (2009) U. Nierste, in  _Helmholz International Summer School on Heavy Quark Physics_ (2009) pp. 1–38, arXiv:0904.1869 [hep-ph] .
  * De Gouvea (2002) A. De Gouvea, Phys. Rev. D 66, 076005 (2002), arXiv:hep-ph/0204077 .
  * Bahcall  _et al._ (2002) J. N. Bahcall, V. Barger, and D. Marfatia, Phys. Lett. B 534, 120 (2002), arXiv:hep-ph/0201211 .
  * Kostelecky and Mewes (2004a) V. A. Kostelecky and M. Mewes, Phys. Rev. D 70, 031902 (2004a), arXiv:hep-ph/0308300 .
  * Kostelecky and Mewes (2004b) V. A. Kostelecky and M. Mewes, Phys. Rev. D 69, 016005 (2004b), arXiv:hep-ph/0309025 .
  * Datta  _et al._ (2004) A. Datta, R. Gandhi, P. Mehta, and S. U. Sankar, Phys. Lett. B 597, 356 (2004), arXiv:hep-ph/0312027 .
  * Minakata and Uchinami (2005) H. Minakata and S. Uchinami, Phys. Rev. D 72, 105007 (2005), arXiv:hep-ph/0505133 .
  * Kostelecky and Mewes (2012) A. Kostelecky and M. Mewes, Phys. Rev. D 85, 096005 (2012), arXiv:1112.6395 [hep-ph] .
  * Ohlsson and Zhou (2015) T. Ohlsson and S. Zhou, Nucl. Phys. B 893, 482 (2015), arXiv:1408.4722 [hep-ph] .
  * Abe  _et al._ (2015) K. Abe _et al._ (Super-Kamiokande), Phys. Rev. D 91, 052003 (2015), arXiv:1410.4267 [hep-ex] .
  * Diaz and Schwetz (2016) J. S. Diaz and T. Schwetz, Phys. Rev. D 93, 093004 (2016), arXiv:1603.04468 [hep-ph] .
  * de Gouvêa and Kelly (2017) A. de Gouvêa and K. J. Kelly, Phys. Rev. D 96, 095018 (2017), arXiv:1709.06090 [hep-ph] .
  * Barenboim  _et al._ (2018) G. Barenboim, C. A. Ternes, and M. Tórtola, Phys. Lett. B 780, 631 (2018), arXiv:1712.01714 [hep-ph] .
  * Liberati (2013) S. Liberati, Class. Quant. Grav. 30, 133001 (2013), arXiv:1304.5795 [gr-qc] .
  * Antusch and Fernandez-Martinez (2008) S. Antusch and E. Fernandez-Martinez, Phys. Lett. B 665, 190 (2008), arXiv:0804.2820 [hep-ph] .
  * Amole  _et al._ (2014) C. Amole _et al._ (ALPHA), Nature Commun. 5, 3955 (2014).
  * Ahmadi  _et al._ (2016) M. Ahmadi _et al._ (ALPHA), Nature 529, 373 (2016).
  * Bertsche (2018) W. A. Bertsche, Phil. Trans. Roy. Soc. Lond. A 376, 20170265 (2018).
  * Ahmadi  _et al._ (2017) M. Ahmadi _et al._ (ALPHA), Nature 548, 66 (2017).
  * Brown and Gabrielse (1986) L. S. Brown and G. Gabrielse, Rev. Mod. Phys. 58, 233 (1986).
  * Borchert  _et al._ (2022) M. J. Borchert _et al._ (BASE), Nature 601, 53 (2022).
  * Hughes and Holzscheiter (1991) R. J. Hughes and M. H. Holzscheiter, Phys. Rev. Lett. 66, 854 (1991).
  * Caldwell and Dvali (2021) A. Caldwell and G. Dvali, Fortsch. Phys. 69, 2000092 (2021), arXiv:1903.09096 [hep-ph] .
  * Ji (1995) X.-D. Ji, Phys. Rev. Lett. 74, 1071 (1995), arXiv:hep-ph/9410274 .
  * Yang  _et al._ (2015) Y.-B. Yang, Y. Chen, T. Draper, M. Gong, K.-F. Liu, Z. Liu, and J.-P. Ma, Phys. Rev. D 91, 074516 (2015), arXiv:1405.4440 [hep-ph] .
  * Yang  _et al._ (2018) Y.-B. Yang, J. Liang, Y.-J. Bi, Y. Chen, T. Draper, K.-F. Liu, and Z. Liu, Phys. Rev. Lett. 121, 212001 (2018), arXiv:1808.08677 [hep-lat] .
  * Baikov  _et al._ (2014) P. A. Baikov, K. G. Chetyrkin, and J. H. Kühn, JHEP 10, 076 (2014), arXiv:1402.6611 [hep-ph] .
  * Aoki  _et al._ (2022) Y. Aoki _et al._ (Flavour Lattice Averaging Group (FLAG)), Eur. Phys. J. C 82, 869 (2022), arXiv:2111.09849 [hep-lat] .
  * Gupta  _et al._ (2021) R. Gupta, S. Park, M. Hoferichter, E. Mereghetti, B. Yoon, and T. Bhattacharya, Phys. Rev. Lett. 127, 242002 (2021), arXiv:2105.12095 [hep-lat] .
  * Bell and Steinberger (1966) J. S. Bell and J. Steinberger, in  _Oxford International Conference on Elementary Particles_ (1966) pp. 195–222.
  * Aguilar  _et al._ (2019) A. C. Aguilar _et al._ , Eur. Phys. J. A 55, 190 (2019), arXiv:1907.08218 [nucl-ex] .
  * Tórtola  _et al._ (2020) M. A. Tórtola, G. Barenboim, and C. A. Ternes, JHEP 07, 155 (2020), arXiv:2005.05975 [hep-ph] .
  * Murayama (2004) H. Murayama, Phys. Lett. B 597, 73 (2004), arXiv:hep-ph/0307127 .
  * Tomboulis (2015) E. T. Tomboulis, Phys. Rev. D 92, 125037 (2015), arXiv:1507.00981 [hep-th] .
  * Greenberg (2002) O. W. Greenberg, Phys. Rev. Lett. 89, 231602 (2002), arXiv:hep-ph/0201258 .
  * Chaichian  _et al._ (2013b) M. Chaichian, K. Fujikawa, and A. Tureanu, Eur. Phys. J. C 73, 2349 (2013b), arXiv:1205.0152 [hep-th] .
  * Larin (2020) S. A. Larin, EPL 129, 21003 (2020), arXiv:2006.00910 [hep-ph] .
  * Kostelecky and Samuel (1989a) V. A. Kostelecky and S. Samuel, Phys. Rev. D 39, 683 (1989a).
  * Kostelecky and Samuel (1989b) V. A. Kostelecky and S. Samuel, Phys. Rev. D 40, 1886 (1989b).
  * Wald (1980) R. M. Wald, Phys. Rev. D 21, 2742 (1980).


# Fusion of Gaussian Processes Predictions with Monte Carlo Sampling

###### Abstract

In science and engineering, we often work with models designed for accurate
prediction of variables of interest. Recognizing that these models are
approximations of reality, it becomes desirable to apply multiple models to
the same data and integrate their outcomes. In this paper, we operate within
the Bayesian paradigm, relying on Gaussian processes as our models. These
models generate predictive probability density functions (pdfs), and the
objective is to integrate them systematically, employing both linear and log-
linear pooling. We introduce novel approaches for log-linear pooling,
determining input-dependent weights for the predictive pdfs of the Gaussian
processes. The aggregation of the pdfs is realized through Monte Carlo
sampling, drawing samples of weights from their posterior. The performance of
these methods, as well as those based on linear pooling, is demonstrated using
a synthetic dataset.

##  1 Introduction

In science and engineering, dealing with models that aim to accurately predict
a variable of interest is a common task. Since these models often are only
approximations of reality, it is prudent to apply multiple models to the same
data and somehow integrate their predictions. Within the Bayesian framework,
results based on models are typically presented as predictive probability
distributions — in the case of continuous variables, they are commonly in the
form of predictive probability density functions (pdfs). In such settings, it
is important to combine the produced predictive pdfs [1]. In this paper, we
focus on models defined by Gaussian processes (GPs) and propose an approach to
fuse the predictive pdfs of an ensemble of different GPs. This involves
aggregating Gaussians obtained through various instances of linear and log-
linear pooling of the Gaussian pdfs. The different instances correspond to
different weights assigned to the GPs, where the weights are drawn by Monte
Carlo sampling from the posterior of the weights.

GPs are probabilistic non-parametric models that find the function that
probabilistically maps an input variable $\mathbf{x}$ to an output variable
$y$. The density of the output $y$ is Gaussian with parameters that depend on
$\mathbf{x}$. This model is inaccurate when the output is not Gaussian. The
effectiveness of model ensembles in enhancing predictive performance has long
been recognized [2]. However, understanding the conditions under which an
ensemble of GPs is likely to yield reliable estimates is not straightforward.
Bayesian Model Averaging (BMA) [3] operates under the assumption that the true
model lies within the hypothesis class of the prior. However, its
suboptimality outside this case is shown empirically in [4], with an extreme
(synthetic) example provided in [5].

In contrast, ensembles adopt a different strategy by focusing on model
combinations. By merging multiple models, ensembles create a more robust and
powerful model. Ensembles demonstrate particular efficacy when the actual
model lies beyond the hypothesis class. In such instances, the collective
strength derived from the diverse models contributes to an overall improvement
in predictive performance.

One of the most popular ensembling methods for point estimators is stacking
[6], which employs $K$ point estimators denoted as
$f_{k}\colon\mathbf{x}\mapsto y_{k}$. The stacking process involves training
an additional point estimator whose inputs are the outputs of each $f_{k}$,
represented as $\hat{y}=g(\hat{y}_{1},\dots,\hat{y}_{K})$. Stacking is quite
general, and in particular, $g$ can be determined in many ways, such as using
a linear model with weights $w_{k}$, or more sophisticated functions like
neural networks.

_Bayesian stacking_ [7] adapts the stacking approach to probabilistic
estimators by replacing point estimates $y_{k}$ with predictive pdfs
$p_{k}(y_{*}\,\lvert\,{\bf x}_{*})$ and $\hat{y}$ with a mixture $\sum
w_{k}p_{k}(y_{*}\,\lvert\,{\bf x}_{*})$. One step further is Bayesian
hierarchical stacking [8], which incorporates a Bayesian prior on input-
dependent weights, making the final mixture $\sum\hat{w}_{k}({\bf
x}_{*})p_{k}(y_{*}\,\lvert\,{\bf x}_{*})$.

Distinct from stacking is the _mixture of experts (MoE)_ method [9], which
forms a mixture by jointly learning several experts, and a gating network that
determines the weights assigned to the experts. This method is generally more
expensive and difficult to train than stacking, but it allows for more
expressive predictive pdfs.

In Section 2, we present a brief background on GPs and Monte Carlo sampling,
thus providing both the setting and methodological tools for this paper. In
Section 3, we describe several existing strategies for the fusion of GP
predictions, namely Bayesian hierarchical stacking (BHS) and the mixture of GP
experts (MoGPE). We then introduce novel variants of BHS and MoGPE by
replacing their linear pooling with log-linear pooling, resulting in the
Product BHS (P-BHS) and the product of GP experts (PoGPE). Finally, we show
how all four models can be estimated with Monte Carlo sampling and random
Fourier feature-based GPs (denoted as RFF-GPs), and we provide numerical
comparisons between methods.

##  2 Background

We briefly review the basics of GP regression, random Fourier feature
approximation of kernel functions, and Bayesian inference using Monte Carlo
sampling.

###  2.1 Gaussian Processes

A GP is described by a stochastic process, where any finite subset of random
variables jointly conforms to a Gaussian distribution [10]. A GP is fully
specified by its mean $\mu(\mathbf{x})$ and covariance function
$\kappa\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ with hyperparameters
$\boldsymbol{\theta}$. In training a GP, we seek hyperparameters
$\boldsymbol{\theta}$ that maximize the log-marginal likelihood of the
hyperparameters defined by

$\displaystyle\log
p(\mathbf{y}\,\lvert\,\mathbf{X},\boldsymbol{\theta})=-\frac{1}{2}\left(\mathbf{y}^{T}\mathbf{K}^{-1}\mathbf{y}+\log\lvert\mathbf{K}\rvert\right),$
(1)

where ${\bf K}\in\mathbb{R}^{N\times N}$ is a covariance matrix with elements
$\left[\mathbf{K}\right]_{ij}=\kappa(\mathbf{x}_{i},\mathbf{x}_{j})$.
Alternatively, we can place a prior $p(\boldsymbol{\theta})$ on the kernel
hyperparameters and proceed under the Bayesian paradigm.

For a given set of hyperparameters $\boldsymbol{\theta}$, a training set
$\mathbf{X},\mathbf{y}$, and a test input
$\mathbf{x}_{*}\in\mathbb{R}^{d_{x}}$, the GP posterior predictive pdf of the
corresponding function value $f_{*}=f\left(\mathbf{x}_{*}\right)$ is Gaussian
with mean and variance given by:

$\displaystyle\mathbb{E}\left[f_{*}\right]$
$\displaystyle=m\left(\mathbf{x}_{*}\right)=\mathbf{k}_{*}^{T}\mathbf{K}^{-1}\mathbf{y},$
(2) $\displaystyle\operatorname{var}\left[f_{*}\right]$
$\displaystyle=\sigma^{2}\left(\mathbf{x}_{*}\right)=k_{**}-\mathbf{k}_{*}^{T}\mathbf{K}^{-1}\mathbf{k}_{*},$
(3)

respectively, where
$\mathbf{k}_{*}=\kappa\left(\mathbf{X},\mathbf{x}_{*}\right)$ and
$k_{**}=\kappa\left(\mathbf{x}_{*},\mathbf{x}_{*}\right)$.

###  2.2 Random Fourier Feature-based Gaussian Processes

By the very nature of being non-parametric, the number of parameters in a GP
grows with the data. This can easily be seen by reparameterizing the GP prior
for a fixed training set $\mathbf{X}$,

$\displaystyle f(\mathbf{X})$ $\displaystyle=\mathbf{L}\boldsymbol{u},$ (4)
$\displaystyle\boldsymbol{u}$
$\displaystyle\sim\mathcal{N}(\mathbf{0},\mathbb{I}_{N}),$ (5)

where $\mathbf{L}=\mathrm{chol}(\mathbf{K})$ is the Cholesky decomposition of
the kernel matrix and $\mathbb{I}_{N}$ is the $N\times N$ identity matrix.
Moreover, computation of the posterior predictive mean and variance requires
the inversion of an $N\times N$ matrix (cf. Eqs. (2) and (3)). The high
dimensionality makes sampling difficult for large datasets, and the matrix
inversion makes prediction expensive.

If the power spectral density $S(\omega)$ of the GP kernel exists, one way to
alleviate these issues is to take a Monte Carlo approximation of the kernel by
drawing frequencies from the normalized power spectral density of the desired
kernel and using them with the input data to form features. We do this by
sampling $M$ sets of _random frequencies_ from the power spectral density to
approximate the integral

$\kappa({\mathbf{x}},{\mathbf{x}}^{\prime})=\frac{1}{2\pi}\int_{\boldsymbol{\omega}}e^{i\boldsymbol{\omega}^{\top}({\mathbf{x}}-{\mathbf{x}}^{\prime})}S(\boldsymbol{\omega})\,d\boldsymbol{\omega},$
(6)

where $\boldsymbol{\omega}$ is a vector of frequencies of the same size as
${\bf x}$. The resulting model is a linear model with $2M$ parameters;
therefore, if $2M<N$, we can expect computational savings. For general kernel
machines, this approximation (known as the _random Fourier feature (RFF)
approximation_) was introduced by [11] and was adapted to GP regression in
[12].

If $\boldsymbol{\omega}_{1},\dots,\boldsymbol{\omega}_{M}$ are the vectors of
frequencies sampled from $S(\boldsymbol{\omega})$, we create a feature vector
according to

${\boldsymbol{\phi}}({\mathbf{x}})=[\cos(\boldsymbol{\omega}_{1}^{\top}{\mathbf{x}})\,\sin(\boldsymbol{\omega}_{1}^{\top}{\mathbf{x}})\,\cdots\,\cos(\boldsymbol{\omega}_{M}^{\top}{\mathbf{x}})\,\sin(\boldsymbol{\omega}_{M}^{\top}{\mathbf{x}})]^{\top}.$
(7)

Then using this reparameterization, the resulting model is a RFF-GP given by

$\displaystyle f({\mathbf{x}})$
$\displaystyle={\boldsymbol{\phi}}({\mathbf{x}})^{\top}\boldsymbol{\psi},$ (8)
$\displaystyle\boldsymbol{\psi}$
$\displaystyle\sim\mathcal{N}(\mathbf{0},\sigma^{2}_{\psi}\mathbb{I}_{2M}),$
(9)

where we assume that the vector of linear parameters $\boldsymbol{\psi}$ has a
zero mean Gaussian prior with a covariance matrix
$\sigma^{2}_{\psi}\mathbb{I}_{2M}$.

###  2.3 Bayesian Inference via Monte Carlo Algorithms

In Bayesian inference, the goal is to infer the parameters of a statistical
model for the observations by computing the posterior distribution. Let
${\boldsymbol{\eta}}$ be a parameter vector, and let $p({\bf
y}\,\lvert\,{\boldsymbol{\eta}},{\bf X})$ and $p({\boldsymbol{\eta}})$ be the
likelihood and the prior of ${\boldsymbol{\eta}}$, respectively. As we will
see below, in this work, ${\boldsymbol{\eta}}$ contains all the parameters of
the RFF-GPs. The posterior $p({\boldsymbol{\eta}}\,\lvert\,{\bf y},{\bf
X})\propto p({\bf y}\,\lvert\,{\boldsymbol{\eta}},{\bf
X})p({\boldsymbol{\eta}})$ is usually intractable, and we need to resort to
approximate inference methods, such as Markov chain Monte Carlo (MCMC) which
computes an approximation based on samples [13]. In our work, we use
Hamiltonian Monte Carlo (HMC) and its adaptively tuned variant “No-U-Turn
Sampler” (NUTS) [14], which are state-of-the-art gradient-based MCMC
algorithms. Given $N_{s}$ samples
$\\{{\boldsymbol{\eta}}^{(i)}\\}_{i=1}^{N_{s}}$ from
$p({\boldsymbol{\eta}}\,\lvert\,{\bf y},{\bf X})$, we can approximate the
posterior predictive pdf of a new observation $y_{*}$ as follows:

$\displaystyle p(y_{*}\,\lvert\,{\mathbf{x}}_{*},{\bf y},{\bf
X})\approx\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}p(y_{*}\,\lvert\,{\mathbf{x}}_{*},{\boldsymbol{\eta}}^{(i)}),\quad{\boldsymbol{\eta}}^{(i)}\sim
p({\boldsymbol{\eta}}\,\lvert\,{\bf y},{\bf X}).$ (10)

##  3 Fusion of Gaussian processes predictions

$y_{n}$$y_{nk}$$f_{nk}$$w_{nk}$$x_{n}$$N$$K$ (a) Stacking

$y_{n}$$y_{nk}$$f_{nk}$$w_{nk}$$x_{n}$$N$$K$ (b) Joint Learning

Fig. 1: Bayesian plate diagrams of joint learning and stacking for Gaussian
process fusion. The only difference is whether the stochastic function $f_{k}$
is treated as an input to be conditioned on (as in (a)), or a random variable
(as in (b)).

Let us consider the case where we have $K$ different GPs, each with a
predictive pdf
$\mathcal{N}\big{(}y\,\lvert\,\mu_{k}({\mathbf{x}}),\sigma^{2}_{k}({\mathbf{x}})\big{)}$
of $y$ conditioned on ${\mathbf{x}}$. We aim to combine these predictive pdfs,
$p_{k}(y\,\lvert\,{\mathbf{x}})$ for $k=1,\dots,K$ to a single, fused pdf
$p_{\text{f}}(y\,\lvert\,{\mathbf{x}})$. In this section, we review the
concepts of _Bayesian hierarchical stacking_ (Section 3.1) and the _mixture of
GP experts_ (Section 3.2), which both use input-dependent weights
$\mathbf{w}({\mathbf{x}})\in\mathbb{S}_{K}$ and a linear pooling function to
result in a fused pdf $p_{\text{f}}(y\,\lvert\,{\mathbf{x}})=\sum
w_{k}({\mathbf{x}})p_{k}(y\,\lvert\,{\mathbf{x}})$. We then present a novel
version of both stacking and mixture of GP experts by using a log-linear
pooling function (Section 3.3).

###  3.1 Bayesian Hierarchical Stacking

In _Bayesian hierarchical stacking (BHS)_ [8], one fuses $K$ different
predictive pdfs from pre-trained algorithms. We assume that these predictive
pdfs come from $K$ GPs, and they are
$p_{k}(y\,\lvert\,{\mathbf{x}})=\mathcal{N}\big{(}y\,\lvert\,\mu_{k}({\mathbf{x}}),\sigma^{2}_{k}({\mathbf{x}})\big{)},$
all obtained by using Eqs. (2) and (3). BHS creates fused predictive pdfs
using a linear pooling function, i.e.,

$p_{\text{lin}}\big{(}y\,\lvert\,{\mathbf{x}},\mathbf{w}({\mathbf{x}})\big{)}=\sum_{k=1}^{K}w_{k}({\mathbf{x}})\mathcal{N}\big{(}y\,\lvert\,\mu_{k}({\mathbf{x}}),\sigma^{2}_{k}({\mathbf{x}})\big{)},$
(11)

where $\mathbf{w}({\mathbf{x}})$ belongs to the simplex $\mathbb{S}_{K}$,
i.e., $\sum_{k=1}^{K}w_{k}({\mathbf{x}})=1$ for all ${\mathbf{x}}$.

BHS lays out a general framework to obtain these weights, but for our work, we
learn the input-dependent weights with an additional GP trained on a separate
stacking dataset $\mathcal{D}_{s}$. The full model becomes [8, Section 5.2],

$\displaystyle\widetilde{w}_{k}(\mathbf{x})$
$\displaystyle\sim\mathcal{GP}(0,\kappa_{k}(\mathbf{x},\mathbf{x}^{\prime})),\enskip\widetilde{w}_{K}({\mathbf{x}})=0,\;\;k=1,\dots,K-1,$
(12) $\displaystyle{\bf w}(\mathbf{x})$
$\displaystyle=\text{softmax}(\widetilde{{\bf w}}(\mathbf{x})),$ (13)
$\displaystyle y\,\lvert\,{\bf w}(\mathbf{x}),\mathbf{x}$
$\displaystyle\sim\sum_{k=1}^{K}w_{k}(\mathbf{x}){\mathcal{N}}\big{(}y\,\lvert\,\mu_{k}({\mathbf{x}}),\sigma_{k}^{2}({\mathbf{x}})\big{)}.$
(14)

We can perform inference with this model by sampling from the posterior
$p({\bf w}(\cdot)\,\lvert\,{\mathcal{D}}_{s})$ via HMC. The final predictive
pdf for an input ${\mathbf{x}}_{*}$ is

$\displaystyle
p_{\text{BHS}}(y_{*}\,\lvert\,{\mathbf{x}}_{*})=\sum_{k=1}^{K}\bar{w}_{k}({\mathbf{x}}_{*}){\mathcal{N}}(y\,\lvert\,\mu_{k}({\mathbf{x}}_{*}),\sigma_{k}^{2}({\mathbf{x}}_{*})),$
(15)

where $\bar{w}_{k}({\mathbf{x}}_{*})$ denotes the MMSE estimate/posterior mean
of $w_{k}(\cdot)$ at ${\mathbf{x}}_{*}$. Because of linearity of the integral
and the independence of $w_{k}$ from $y\,\lvert\,\mu_{k},\sigma_{k}$, the use
of the posterior mean $\bar{w}_{k}$ provides the same estimator as the full
approximation

$p_{\text{BHS}}(y_{*}\,\lvert\,{\mathbf{x}}_{*})=\frac{1}{N}\sum_{n=1}^{N}\sum_{k=1}^{K}w_{k,n}({\mathbf{x}}_{*}){\mathcal{N}}(y\,\lvert\,\mu_{k}({\mathbf{x}}_{*}),\sigma_{k}^{2}({\mathbf{x}}_{*})).$
(16)

###  3.2 Mixture of GP Experts

An important feature of the BHS strategy is that the $K$ GPs are fixed during
the learning of the weights. In other words, the functions
$\mu_{k}({\mathbf{x}})$ and $\sigma^{2}({\mathbf{x}})$ are computed using an
initial training dataset $\mathcal{D}_{0}$ and kept fixed from there on. To
create a more flexible model, we can instead use a single dataset
$\mathcal{D}_{\text{train}}$ to perform joint training of the $K$ GP experts
and the weights. This joint learning falls into the mixture of experts
framework [9], and was first proposed for GPs in [15].

In order to have flexible experts, [15] uses GP experts with input-dependent
noise; this can be performed by placing GP priors on both the mean
$\mu_{k}({\mathbf{x}})$ and log-scale $\log\sigma_{k}({\mathbf{x}})$ [16].
Again using the linear pooling function, we arrive at the _mixture of GP
experts_ (MoGPE) defined by

$\displaystyle\mu_{k}({\mathbf{x}})$
$\displaystyle\sim\mathcal{GP}(0,\kappa_{k}^{\mu}),\ k=1,\dots,K,$ (17)
$\displaystyle\log\sigma_{k}({\mathbf{x}})$
$\displaystyle\sim\mathcal{GP}(0,\kappa_{k}^{\sigma}),\ k=1,\dots,K,$ (18)
$\displaystyle\widetilde{w}_{k}(\mathbf{x})$
$\displaystyle\sim\mathcal{GP}(0,\kappa^{w}_{k}),\enskip
k=1,\dots,K-1;\widetilde{w}_{K}({\mathbf{x}})=0,$ (19) $\displaystyle{\bf
w}(\mathbf{x})$ $\displaystyle=\text{softmax}(\widetilde{{\bf
w}}(\mathbf{x})),$ (20) $\displaystyle y\,\lvert\,{\bf z},{\mathbf{x}}$
$\displaystyle\sim\sum_{k=1}^{K}w_{k}(\mathbf{x}){\mathcal{N}}(y\,\lvert\,\mu_{k}({\mathbf{x}}),\exp(2\log\sigma_{k}({\mathbf{x}}))),$
(21)

where ${\bf
z}=\\{w_{k}({\mathbf{x}}),\mu_{k}({\mathbf{x}}),\log\sigma_{k}({\mathbf{x}})\\}_{k=1}^{K}$.
In order to obtain the final predictive pdf, we marginalize the likelihood in
Eq. (21) with respect to the posterior of $w_{k}(\cdot)$, $\mu_{k}(\cdot)$,
and $\log\sigma_{k}(\cdot)$ for $k=1,\dots,K$.

We can again sample the unknowns of the model with HMC. However, the resulting
posterior of the unknowns is much more complex to sample than the BHS
posterior — as a result, we propose using RFF GPs for means, log-scales, and
non-normalized weights for scalable inference.

###  3.3 Log-Linear Pooling for Stacking and Joint Learning

While the BHS and MOGPE methods focused on linear pooling of GP predictions,
one can instead use a log-linear pooling rule, where a weighted average of the
log-pdfs is taken instead of the pdfs. The log-linear pooling scheme was first
analyzed in [17], with an axiomatic characterization and comparison available
in [1], and it possesses analytically different properties than linear
pooling. As the weighted average of log-pdfs is equivalent to a weighted
geometric average, we can express the resulting pdf as a weighted product,

$p_{\text{log-
lin}}(y\,\lvert\,{\mathbf{x}})=c\big{(}\mathbf{w}({\mathbf{x}})\big{)}\prod_{k=1}^{K}p_{k}(y\,\lvert\,{\mathbf{x}})^{w_{k}({\mathbf{x}})},$
(22)

where $c\big{(}\mathbf{w}(x)\big{)}$ is a normalizing constant, and
$\mathbf{w}$ is typically taken to be in the simplex $\mathbb{S}_{K}$.

When the predictive pdfs $p_{k}(\cdot)$ are Gaussian, it is well known that
the pdf resulting from Eq. (22) is proportional to another Gaussian. In the
context of fusing GPs, the resulting rule is known as the _generalized product
of experts (gPoE)_ [18], and it is defined by

$\displaystyle p_{\text{gPoE}}(y\,\lvert\,{\mathbf{x}})$
$\displaystyle=\mathcal{N}\big{(}y\,\lvert\,\mu_{\text{gPoE}},\sigma^{2}_{\text{gPoE}}\big{)},$
(23) $\displaystyle\mu_{\text{gPoE}}(\mathbf{x})$
$\displaystyle=\sigma^{2}_{\text{gPoE}}\sum_{k=1}^{K}w_{k}({\mathbf{x}})\sigma_{k}^{-2}({\mathbf{x}})\mu_{k}({\mathbf{x}}),$
(24) $\displaystyle\sigma^{-2}_{\text{gPoE}}({\mathbf{x}})$
$\displaystyle=\sum_{k=1}^{K}w_{k}({\mathbf{x}})\sigma^{-2}_{k}({\mathbf{x}}).$
(25)

Several strategies for setting the weights $\mathbf{w}$ have been proposed in
the literature [18, 19]. Though these existing strategies have chosen weights
belonging to the simplex $\mathbb{S}_{K}$, this is not necessary for the
mathematical validity of Eq. (23). Instead, it has been chosen for the
property that the resulting precision
$\sigma_{\text{gPoE}}^{-2}({\mathbf{x}})$ is a convex combination of the
expert precisions $\sigma^{-2}_{k}$.

The existing literature surrounding gPoE has not trained a separate
statistical model for input-dependent weights $\mathbf{w}(\cdot)$, nor has it
considered joint learning in such a model. To this end, we propose _product-
BHS (P-BHS)_ , which performs Bayesian hierarchical stacking with a log-linear
pooling rule, and the _product of GP experts (PoGPE)_ , which performs joint
learning.

Concretely, we propose the following hierarchical model for the input-
dependent weights $\mathbf{w}(\cdot)$:

$\displaystyle\log w_{k}(\mathbf{x})$
$\displaystyle\sim\mathcal{GP}(-\log(K),\kappa_{k}(\mathbf{x},\mathbf{x}^{\prime})),\enskip
k=1,\dots,K,$ (26) $\displaystyle y\,\,\lvert\,\,{\bf
w}(\mathbf{x}),\mathbf{x}$
$\displaystyle\sim\mathcal{N}(y\,\lvert\,\mu_{\text{gPoE}}(\mathbf{x};{\bf
w}(\mathbf{x})),\sigma^{2}_{\text{gPoE}}(\mathbf{x};{\bf w}(\mathbf{x}))),$
(27)

where we have made explicit the dependence of $\mu_{\text{gPoE}}$ and
$\sigma^{2}_{\text{gPoE}}$ on $\mathbf{w}({\mathbf{x}})$. Note that we model
the unconstrained weights $\log w_{k}({\mathbf{x}})\in\mathbb{R}$, meaning we
do not restrict weights to $\mathbb{S}_{K}$ but instead $\mathbb{R}^{+}$. This
has the effect of allowing a heteroscedastic prior predictive of our model,
with the use of stationary GP experts. Moreover, the GP prior mean of $-\log
K$ implies a prior mean of $1/K$ on $w_{k}$. Therefore, we avoid issues of
being overconfident as ${\mathbf{x}}$ falls outside the typical set of our
data.

The final model can be written as

$\displaystyle\mu_{k}({\mathbf{x}})$
$\displaystyle\sim\mathcal{GP}(0,\kappa_{k}^{\mu}),\ k=1,\dots,K,$ (28)
$\displaystyle\log\sigma_{k}({\mathbf{x}})$
$\displaystyle\sim\mathcal{GP}(0,\kappa_{k}^{\sigma}),$ (29)
$\displaystyle\log w_{k}(\mathbf{x})$
$\displaystyle\sim\mathcal{GP}(-\log(K),\kappa^{w}_{k}),\enskip k=1,\dots,K,$
(30) $\displaystyle y\,\,\lvert\,\,{\bf z},{\mathbf{x}}$
$\displaystyle\sim{\mathcal{N}}(y\,\lvert\,\mu_{\text{gPoE}}(\mathbf{x};{\bf
w}({\mathbf{x}})),\sigma^{2}_{\text{gPoE}}({\mathbf{x}};{\bf
w}({\mathbf{x}}))),$ (31)

where ${\bf z}=\\{\mu_{k}({\mathbf{x}}),\log\sigma_{k}({\mathbf{x}})\log
w_{k}({\mathbf{x}})\\}_{k=1}^{K}$. The variables $\mu_{k}({\mathbf{x}})$ and
$\ \log\sigma_{k}({\mathbf{x}})$ are given in P-BHS and inferred jointly in
POGPE. The predictive pdf and its Monte Carlo approximation can then be
written as

$\displaystyle p(y_{*}\,\lvert\,\mathbf{x}_{*})$
$\displaystyle\approx\frac{1}{J}\sum_{j=1}^{J}\mathcal{N}(y_{*}\,\lvert\,\mu^{(j)}_{\text{gPoE}},\sigma^{2,(j)}_{\text{gPoE}}),$
(32)

where $\mu^{(j)}_{\text{gPoE}}$ and $\sigma^{2,(j)}_{\text{gPoE}}$ are
obtained by computing Eq. (23) using the posterior samples of $\\{\log
w_{k}({\mathbf{x}}_{*})\\}_{k=1}^{K}$ in P-BHS or the posterior samples of
$\\{\mu_{k}({\mathbf{x}}_{*}),\log\sigma_{k}({\mathbf{x}}_{*}),\log
w_{k}({\mathbf{x}}_{*})\\}_{k=1}^{K}$ in PoGPE.

##  4 Experiments

Fig. 2: A depiction of the generative model and the point predictions of each
method. From top to bottom: the means of each component GP and the respective
weight process; sample draws from each GP with linewidth indicating
$w_{k}(x)$; the resulting dataset with the predictions of each method.

In this section, we run a one-dimensional numerical experiment where we test
the novel P-BHS and PoGPE. Results are compared to BHS, MoGPE, and a
heteroscedastic RFF GP (_het-RFF-GP_). We consider experiments where the
number of spectral frequencies $M$ is fixed and the number of experts $K$ is
varied (Section 4.1), and experiments where $K$ is fixed and $M$ is varied
(Section 4.2).

We generate data following the generative model in Section 3.2 (namely,
according to MoGPE). The weights for three functions,
$\widetilde{w}_{k}(\mathbf{x})$, $K=3$ are generated from a GP with zero mean
and RBF kernel with lengthscale $\ell=0.5$. Three functions are sampled from
another GP with RBF kernel and lengthscale $\ell=0.2$. Given
${\mathbf{x}}_{i}$, we first sample the unconstrained weights, apply softmax,
and we sample the $i$-th observation from
${\mathcal{N}}(y\,\lvert\,\mu_{k}({\mathbf{x}}_{i}),\exp(2\log\sigma_{k}({\mathbf{x}}))$
with probability $w_{k}({\mathbf{x}}_{i})$. Hence, conditional on GP draws,
the function of interest is a linear mixture. The experts, input-dependent
weights, and the resulting dataset are shown in Figure 2. We generate $1,000$
data points according to this procedure and consider an 80%-20% training-test
split, namely, $N_{\text{tr}}=800$ and $N_{\text{tst}}=200$. We compare the
methods in terms of their mean negative log-predictive density (NLPD) on test
data, averaged over 5 random splits. All models were implemented in the
probabilistic programming language NumPyro [20], whose NUTS implementation was
used for inference. In each case, $500$ samples were drawn from the posterior
of each of the four chains, resulting in a total of $J=2000$ samples.

(a)

(b)

Fig. 3: Results from the experiments described in (a) Section 4.1 and (b)
Section 4.2. The NLPD _(lower is better)_ are averaged across test data; boxes
show the median and quartiles of the average across $5$ random data splits for
each fusion method, with whiskers showing the minimum and maximum values. In
(a), the NLPD of the het-RFF-GP is denoted as a line, as it does not change
with $K$.

###  4.1 Performance vs number of experts

In the first experiment, we fix the number of spectral frequencies to $M=30$
and perform inference using each method with varying the number of experts
$K$. Recall that, while in in BHS and P-BHS, the experts are pre-trained in a
separate initial dataset, in MoGPE and PoGPE the experts are jointly inferred
with the weights. For the former pair, we further divide the training dataset
into two parts, using the first half to train the $K$ experts, and the second
half to train the stacking models. To train the experts, each expert is
randomly assigned $\lfloor\frac{N_{\text{tr}}}{K}\rfloor$ data points, with
hyperparameters being selected by type-II maximum likelihood. Note that this
means experts receive fewer data points as $K$ increases. The results are
shown in Figure 3(a).

Our results show that using more experts allows for higher flexibility in
MoGPE and PoGPE, with performance increasing as $K$ ranges from $2$ to $5$.
MoGPE and PoGPE perform comparably with $K=2$, with MoGPE gaining advantages
as more experts are available. For the stacking methods, the predictive power
decreases as the number of experts increases. This can be seen as an artifact
of how we trained experts, with a tradeoff between the number of experts and
their quality. The results show that P-BHS outperforms BHS for all tested
values of $K$, suggesting that stacking with log-linear pooling may result in
better performance. This is likely due to the property that log-linear pooling
with unconstrained weights can arbitrarily adjust the predictive variances.
Intuitively, joint learning (with linear or log-linear pooling) beats stacking
in either case.

###  4.2 Performance vs number of spectral frequencies

Since the quality of the RFF-GP approximation depends on the number of
spectral frequencies $M$, we compare the performance of each method for
different values of $M$, with a fixed number of experts $K=2$. The results are
shown in Figure 3(b). For MoGPE and PoGPE, as well as the het-RFF-GP, there is
a strong benefit to increasing the number of spectral frequencies, with the
predictive pdf improving by nearly an order of magnitude as $M$ ranges from
$10$ to $50$. On the other hand, stacking methods do not benefit much from
increasing the number of spectral frequencies; we expect this to be because
both experts are already quite accurate, making the stacking GPs rather
simple.

##  5 Conclusions

The fusion of GP predictions has received significant attention in the
literature, but primarily from the perspective of linear pooling. In the case
of log-linear pooling, there has been a lack of a principled method for
determining input-dependent weights $\mathbf{w}({\mathbf{x}})$. In this paper,
we cast existing approaches under an inclusive framework and introduce novel
approaches based on log-linear pooling. We show how the input-dependent
weights can be derived in a principled manner in a Bayesian hierarchical
model. Furthermore, we analyze the empirical performance of each method on a
synthetic example. Future directions of research include additional empirical
validation and mathematical analysis of log-linear fusion with unconstrained
weights.

## References

  * [1] G. Koliander, Y. El-Laham, P. M. Djurić, and F. Hlawatsch,  “Fusion of probability density functions,”  Proceedings of the IEEE, vol. 110, no. 4, pp. 404–453, 2022.
  * [2] T. G. Dietterich,  “Ensemble methods in machine learning,”  in International Workshop on Multiple Classifier Systems. Springer, 2000, pp. 1–15.
  * [3] J. M. Bernardo and A. F.M. Smith,  Bayesian Theory, vol. 405,  John Wiley & Sons, 2009.
  * [4] P. Domingos,  “Bayesian averaging of classifiers and the overfitting problem,”  in Proceedings of the Seventeenth International Conference on Machine Learning, 2000, pp. 223–230.
  * [5] T. P. Minka,  “Bayesian model averaging is not model combination,”  Tech. Rep., MIT Media Lab, 2000.
  * [6] D. H. Wolpert,  “Stacked generalization,”  Neural Networks, vol. 5, no. 2, pp. 241–259, 1992.
  * [7] Y. Yao, A. Vehtari, D. Simpson, and A. Gelman,  “Using stacking to average bayesian predictive distributions (with discussion),”  Bayesian Analysis, vol. 13, no. 3, pp. 917–1007, 2018.
  * [8] Y. Yao, G. Pirš, A. Vehtari, and A. Gelman,  “Bayesian hierarchical stacking: Some models are (somewhere) useful,”  Bayesian Analysis, vol. 17, no. 4, pp. 1043–1071, 2022.
  * [9] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton,  “Adaptive mixtures of local experts,”  Neural Computation, vol. 3, no. 1, pp. 79–87, 1991.
  * [10] C. E. Rasmussen and C. K. I. Williams,  Gaussian Processes for Machine Learning,  Adaptive Computation and Machine Learning series. MIT Press, 2005.
  * [11] A. Rahimi and B. Recht,  “Random features for large-scale kernel machines,”  Advances in Neural Information Processing Systems, vol. 20, 2007\.
  * [12] M. Lázaro-Gredilla, J. Quińonero-Candela, C. E. Rasmussen, and A. R. Figueiras-Vidal,  “Sparse spectrum Gaussian process regression,”  The Journal of Machine Learning Research, vol. 11, pp. 1865–1881, 2010.
  * [13] C. P. Robert and G. Casella,  Monte Carlo Statistical Methods,  Springer, 2004.
  * [14] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth,  “Hybrid Monte Carlo,”  Physics Letters B, vol. 195, no. 2, pp. 216–222, 1987.
  * [15] V. Tresp,  “Mixtures of Gaussian processes,”  Advances in Neural Nnformation Processing Systems, vol. 13, 2000\.
  * [16] P. Goldberg, C. Williams, and C. Bishop,  “Regression with input-dependent noise: A Gaussian process treatment,”  Advances in Neural Information Processing Systems, vol. 10, 1997\.
  * [17] C. Genest,  “A characterization theorem for externally Bayesian groups,”  The Annals of Statistics, pp. 1100–1105, 1984.
  * [18] Y. Cao and D. J. Fleet,  “Generalized product of experts for automatic and principled fusion of Gaussian process predictions,”  arXiv preprint arXiv:1410.7827, 2014.
  * [19] M. Deisenroth and J. W. Ng,  “Distributed Gaussian processes,”  in International Conference on Machine Learning. PMLR, 2015, pp. 1481–1490.
  * [20] D. Phan, N. Pradhan, and M. Jankowiak,  “Composable effects for flexible and accelerated probabilistic programming in NumPyro,”  arXiv preprint 1912.11554, 2019.


# Case study: Mapping potential informal settlements areas in Tegucigalpa with
machine learning to plan ground survey

Federico Bayle  Dymaxion LabsBuenos AiresArgentina<EMAIL_ADDRESS>and  Damian E. Silvani  Dymaxion LabsBuenos AiresArgentina
<EMAIL_ADDRESS>
(2020)

###### Abstract.

Data collection through censuses is conducted every 10 years on average in
Latin America, making difficult monitoring the growth and support needed by
communities living in informal settlements. Conducting a field survey requires
logistical resources to be able to do it exhaustively. The increasing
availability of spatial open-data, high-resolution satellite images, and open-
source tools, allow us to train machine learning algorithms map these areas in
different cities of Latin America. This case study shows the collaboration
between Dymaxion Labs and the NGO Techo to employ that techniques to create
the first informal settlements census of Tegucigalpa, Honduras.

Neural networks, informal settlements, transfer learning, GIS.

††copyright: acmcopyright††journalyear: 2020††conference: ACM SIGKDD
Conference on Knowledge Discovery and Data Mining ; August 24, 2020; San
Diego, California USA††booktitle: KDD Humanitarian Mapping Wokshop ’20: ACM
SIGKDD Conference on Knowledge Discovery and Data Mining August 24, 2020; San
Diego, California USA

##  1\. Introduction

Poverty is a major issue in Latin America. Slums and informal settlements grow
quickly, and there is a need for more up-to-date data for policymaker
decisions. If governments had updated information about slums and their
growth, they could give the affected families better life conditions and
change their future. Having updated information helps to improve health,
education, and security for the children who currently live there. In order to
carry out an accurate survey of informal areas, exhaustive coverage of the
territory is required. This task requires having costly logistical and
material resources, which in turn undermines periodicity and scope.

These settlements often do not have critical public services such as
sanitation, resulting in health and environmental hazards, especially for
children.

Tegucigalpa (formally Tegucigalpa, Municipality of the Central District) is
Honduras’ largest and most populous city as well as the nation’s political and
administrative center. The urban population is estimated to 1,143,373
inhabitants by 2019. The urban area is 201.5 square kilometers. There no
significant progress has been seen in this matter in recent decades about
informal settlements life conditions. The data required to formulate public
policies aimed at improving the situation of informal settlements population,
is almost non-existent through official and unofficial sources in Honduras.

For this reason, the NGO TECHO111www.techo.org planned to conduct in 2018 a
census of informal settlements in the urban district of the Central District,
aiming to provide clear, reliable and public information on the current
problem of thousands of people living in poverty. TECHO is a youth-led non-
profit organization founded in 1997 with a presence in Latin America and the
Caribbean. They seek to build a fair, integrated, and poverty-free society,
where everyone has the opportunities needed to develop their capacities and
fully exercise their rights.

Due to the lack of the previous census and public data about informal
settlements in Tegucigalpa, in Dymaxion Labs222www.dymaxionlabs.com we
partnered with TECHO to help them map the potential areas to survey. The main
objective was to reduce the survey area to be scouted by the volunteer team.
Based on our previous work on informal settlements mapping in Buenos Aires
(Argentina) and Asunción (Paraguay), the machine learning algorithm AP-
LATAM(Labs, 2017) looked promising to apply in Honduras based on a previous
survey of Guatemala City (Guatemala). Dymaxion Labs is a startup developing a
cloud-based API to detect objects in satellite imagery at scale.

##  2\. Related Work

Following Kuffer, Pfeffer, and Sliuzas (2016)(Kuffer et al., 2016), at least
87 papers were published in English about this topic. Reviewing some of them,
we can find studies for Argentina, Colombia, England, Ghana, Nigeria,
Tanzania, Uganda, Malawi, and Rwanda.

Focusing on Latin American cities, Hall, Malcolm, and Piwowar (2001)(Hall et
al., 2001) used Landsat and Radarsat imagery combined with GIS data to detect
urban poverty pockets in Rosario, Argentina. Patiño and Duque (2013)(Patino
and Duque, 2013) used very high-resolution images to estimate a slum index on
Medellín, Colombia.

For other cities around the world, Stoler et. al. (2012)(Stoler et al., 2012)
and Weeks et.al (2007)(Stow et al., 2007) used very high-resolution imagery
for Accra (Ghana) deriving features based on texture and land cover to
estimate a slums distribution index. Arribas-Bel, Patiño, and Duque
(2017)(Arribas-Bel et al., 2017) employed a similar approach to estimate the
Living Environment Deprivation index for small areas in Liverpool, England.
Jean et. al. (2016)(Jean et al., 2016) used nighttime imagery with very high-
resolution daily imagery to estimate variation in local-level economic
outcomes for Nigeria, Tanzania, Uganda, Malawi, and Rwanda. They employed
transfer learning of trained convolutional neural networks, a technique we
applied in this case study as well.

##  3\. Methodology

The open-source algorithm AP-LATAM, tries to reduce survey costs by analyzing
high-resolution satellite imagery to detect potential areas of informal
settlement growth. The final result is a geospatial dataset of areas that
could contain informal settlements. By having updated potential slums growth
locations, decision makers can prioritize the areas to survey.

Source code for building and using a classifier to create datasets has also
been released as open-source in a GitHub
repository333https://github.com/dymaxionlabs/ap-latam with a BSD-2 license.
Instructions on how to train the classifier and predicting over new images are
available there. This includes not only the steps for training the model but
also the neural net weights for using with new satellite imagery. AP-LATAM is
part of the Inter-American Development Bank’s Code for Development
platform444https://code.iadb.org/en/repository/60/ap-latam.

This tool is already calibrated to detect slums and informal settlements in
Buenos Aires, Asunción, and Montevideo. To obtain those results, Dymaxion Labs
applied computer vision and deep learning algorithms to satellite imagery and
other georeferenced data. Slum patterns such as texture and morphology of
roofs are detected. From Dymaxion Labs experience calibrating these kinds of
algorithms, such patterns have strong local dependencies. As TECHO has ground
truth polygons and imagery from a previous survey in Guatemala City
(Guatemala) in 2016, the strategy was based on training AP-LATAM on these data
and then mapping potential informal settlements areas on Tegucigalpa
(Honduras). Guatemala City is the capital and largest city of Guatemala, a
country that borders with Honduras. It has 2,750,965 inhabitants with an area
of 220 square kilometers. Taking into account the TECHO’s survey database and
the similarity of the geography of these cities, it was the best candidate to
train AP-LATAM to map Tegucigalpa.

The method consists of a binary classifier of image tiles of high-resolution
sub-meter satellite imagery. Each image is classified as whether it contains
an informal settlement or not. Some examples of these tiles area showed in
figure 1

(a) Tiles tagged as _true_

(b) Tiles tagged as _false_

Figure 1. Example of images tiles tagged as either true or false, whether they
contain a informal settlement or not.

To build the dataset for training and validation, the classifier takes a
vector file of polygons of previously-known informal settlements and takes
fixed-size tiles of images by sliding a window across the entire satellite
image. For each image tile, it checks if the tile intersects with any polygon
and tags it appropriately.

To make predictions over new images, it slides a window over the new image and
builds a new vector file of polygons of the size of each positively-tagged
image tile. The resulting dataset is post-processed to remove polygons with
small probability and dissolve them into bigger polygons.

In our tests with Guatemala, as a final post-processing step, we used
OpenStreetMap(OpenStreetMap contributors, 2017) datasets and crossed them with
the polygons file.

###  3.1. Imbalanced classes

A binary classifier of informal settlements has imbalanced classes, that is,
images tagged as positive (areas that contain settlements) are much less
frequent than images tagged as negative (areas that do not contain
settlements). To decrease bias and avoid overfitting, we under-sampled the
negatives by taking a set of images of size proportional to the size of
positives. We tried undersampling with a proportion of 4 and 8 and settled
with 4.

###  3.2. Data augmentation

To help prevent overfitting and make the model generalize better, we perform
data augmentation on the image tiles. We only applied horizontal and vertical
flipping, but there are other random transformations that we could use to
augment our dataset, like hue and brightness randomization (to account for
differences in atmosphere corrections) and rotations.

###  3.3. Fine-tuning

A large amount of data is needed to build a functional convolutional neural
network model. In practice, it is common to reuse a pre-trained network.
However, most pre-trained networks work for a different set of labels and were
not trained with satellite images, so for this use case it is necessary to
retrain some of the top layers to improve prediction.

The methodology used here was to fine-tune a ResNet-50(He et al., 2015)
network with our satellite imagery. We chose this architecture taking into
account the trade-off between precision and hardware requirements to perform
experiments. The procedure is roughly as follows:

  1. (1)

Instantiate the convolutional base of ResNet-50.

  2. (2)

Add a fully-connected model on top, with a standard SGD optimizer and
configure a binary cross-entropy loss function.

  3. (3)

Freeze the layers of the model up to the top 70 layers.

  4. (4)

Retrain the model.

The Keras library was used for data augmentation, training, and
prediction(Chollet et al., 2015).

###  3.4. Post-processing

The resulting dataset after prediction over sliding windows is a set of small
fixed-size squares, with a prediction probability associated. To refine the
results we apply the following:

  1. (1)

Median filter: remove squares with low probability and a small number of
neighbors.

  2. (2)

Dissolve overlapping squares: if the sliding window step size is smaller than
the size of the windows, it may end up with overlapping squares, so this step
dissolves them into a single polygon with mean probability values between the
values of each connected squares.

As mentioned before, we also used a dataset of blocks and calculated the
intersection between the squares and blocks, and if sufficient squares covered
a block, we picked them to form a new dataset of blocks that contain potential
informal settlements. The outcome has better prediction accuracy mainly
because roads and other areas are not considered.

##  4\. Results

The resultant methodology had a Cohen’s Kappa(Cohen, 1960) index of 0.82 at
pixel level on Guatemala validation polygons. Based on the previous work done
between both organizations analyzing the cities of Buenos Aires (Argentina),
Asunción (Paraguay) and Montevideo (Uruguay), imagery from those countries
were added to the training dataset. Due to the different geography, like the
soil color in Asunción, and the common used materials for construction, the
best results were achieved using only Guatemala data.

After being trained with the entire city, we applied it on the Tegucigalpa
imagery to generate the map of potential areas for surveying. This map (figure
2) was uploaded to our website (both as downloadable data and the online map).
These datasets have been released as public domain data, using the same
license555https://opendatacommons.org/licenses/pddl/1-0/ that OpenStreetMap
uses for its data, the Open Data Commons Public Domain Dedication and License.
From the AP-Latam website the user can download the datasets as GeoJSON files,
one for each area and image acquisition date. The user can also explore an
online map with the latest dataset generated.

Figure 2. Tegucigalpa map of potential informal settlements.

Combining this map and other data sources from a utility company and
government agencies, TECHO planned the fieldwork in two stages: territory
exploration and surveyors training. To start the survey in one of the most
violent cities in Latin America, a strict security protocol was designed. One
of the protocols taken for the safety of the volunteers on the ground
suggested not to survey an informal settlement without having previously
contacted local referents that would allow them to guarantee minimum security.
The search for local referents was carried out throughout the fieldwork. TECHO
also used the potential informal settlements map to double-check the
delimitation of the most hazardous places in the city to avoid sending
volunteers to survey those areas.

##  5\. Conclusions

The application of the transfer learning methodology trained in Guatemala City
to map potential informal settlements in Tegucigalpa collaborates with TECHO’s
plan to survey the entire city. In particular, it helps them to map informal
settlements located in dangerous areas where volunteers are not allowed to
survey.

The 2018 Informal Settlements Census(TECHO, 2018) finally took place in the
Metropolitan Area of the Central District (D.C.), specifically in the urban
area, made up of the cities of Tegucigalpa and Comayaguela. The objective was
to visit 796 neighborhoods and residential areas that make up the urban area,
with the aim of identifying informal settlements and subsequently
characterizing and georeferencing each one of them. The results were released
as open data666http://datos.techo.org/en/dataset/honduras-censo-de-
asentamientos-informales-casco-urbano-distrito-central. 40% of the territory
where informal settlements without property titles are located belong to the
state. Thereby, this map can be used by policymakers to relocate people living
in risky areas immediately.

After the execution, 161 informal settlements were identified, 124 settlement
files were raised with the help of references in informal settlements, and it
was estimated that there are 42,000 families living in this situation.

From a replication perspective, having a technical team is critical not only
for training and applying the algorithm, but also to help validate the mapping
results. With more frequent ground validations, more iterations could be made
on the algorithm to get a better accuracy in the final map. The technical
skills and GIS knowledge of TECHO’s team proved to be vital in achieving the
quality level of ground surveys. In terms of mapping new cities, the
availability of up-to-date high-resolution imagery is also essential to get
relevant results.

###### Acknowledgements.

To María Jesús Silva Rozas and Rubens C. Schrunder for giving us the chance to
apply our methodology in a real-world case and validate the results in
Honduras. To UNICEF Innovation Fund team for their support as part of the Data
Science cohort.

## References

  * (1)
  * Arribas-Bel et al. (2017) Daniel Arribas-Bel, Jorge E Patino, and Juan C Duque. 2017\.  Remote sensing-based measurement of Living Environment Deprivation: Improving classical approaches with machine learning.  _PloS one_ 12, 5 (2017), e0176684.
  * Chollet et al. (2015) François Chollet et al. 2015\.  Keras.  https://keras.io.
  * Cohen (1960) J Cohen. 1960.  A Coefficient of Agreement for Nominal Scales.  _Educational and Psychological Measurement_ 20, 1 (1960), 37–46.  https://doi.org/10.1177/001316446002000104 arXiv:1011.1669v3
  * Hall et al. (2001) G Brent Hall, Neil W Malcolm, and Joseph M Piwowar. 2001\.  Integration of remote sensing and GIS to detect pockets of urban poverty: The case of Rosario, Argentina.  _Transactions in GIS_ 5, 3 (2001), 235–253.
  * He et al. (2015) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2015\.  Deep Residual Learning for Image Recognition.  _arXiv preprint arXiv:1512.03385_ (2015).
  * Jean et al. (2016) Neal Jean, Marshall Burke, Michael Xie, W Matthew Davis, David B Lobell, and Stefano Ermon. 2016\.  Combining satellite imagery and machine learning to predict poverty.  _Science_ 353, 6301 (2016), 790–794.
  * Kuffer et al. (2016) Monika Kuffer, Karin Pfeffer, and Richard Sliuzas. 2016\.  Slums from Space - 15 Years of Slum Mapping Using Remote Sensing.  _Remote. Sens._ 8 (2016), 455.
  * Labs (2017) Dymaxion Labs. 2017\.  AP-Latam.  https://ap-latam.dymaxionlabs.com/.
  * OpenStreetMap contributors (2017) OpenStreetMap contributors. 2017\.  Planet dump retrieved from https://planet.osm.org .  https://www.openstreetmap.org.
  * Patino and Duque (2013) Jorge E. Patino and Juan Carlos Duque. 2013.  A review of regional science applications of satellite remote sensing in urban settings.  _Comput. Environ. Urban Syst._ 37 (2013), 1–17.
  * Stoler et al. (2012) Justin Stoler, Dean Daniels, John R Weeks, Douglas A Stow, Lloyd L Coulter, and Brian Karl Finch. 2012\.  Assessing the utility of satellite imagery with differing spatial resolutions for deriving proxy measures of slum presence in Accra, Ghana.  _GIScience & remote sensing_ 49, 1 (2012), 31–52.
  * Stow et al. (2007) D Stow, A Lopez, C Lippitt, S Hinton, and J Weeks. 2007.  Object-based classification of residential land use within Accra, Ghana based on QuickBird satellite data.  _International journal of remote sensing_ 28, 22 (2007), 5167–5173.
  * TECHO (2018) TECHO. 2018.  Informe: Censo de Asentamientos Informales, Casco Urbano Distrito Central 2018.  (2018).  http://datos.techo.org/dataset/8efe820f-fda2-4ff6-a0a5-6c28b57596e6/resource/5e75ac50-44c9-4450-ab8b-872ba939bbc4/download/informe-censo-de-asentamientos-informales-casco-urbano-dc-honduras-2018-interactivo-light.pdf

##  6\. About the authors

Federico Bayle is Bs. in Economics and MSc. in Data Mining and KDD. His master
thesis was based on a machine learning methodology to map informal settlements
in Buenos Aires, Argentina.

Damian Silvani is Bs. in Computer Sciences, with a specialization in Computer
Vision and Machine Learning. More than 15 years of experience in the software
development industry solving hard problems.