# Modular Gaussian Processes for Transfer Learning

Pablo Moreno-Muñoz∗ Antonio Artés-Rodríguez† Mauricio A. Álvarez‡

∗Section for Cognitive Systems, Technical University of Denmark (DTU)
†Dept. of Signal Theory and Communications, Universidad Carlos III de Madrid,
Spain
†Evidence-Based Behavior (eB2), Spain
‡Dept. of Computer Science, University of Sheffield, UK
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

We present a framework for transfer learning based on modular variational
Gaussian processes (GP). We develop a module-based method that having a
dictionary of well fitted GPs, one could build ensemble GP models without
revisiting any data. Each model is characterised by its hyperparameters,
pseudo-inputs and their corresponding posterior densities. Our method avoids
undesired data centralisation, reduces rising computational costs and allows
the transfer of learned uncertainty metrics after training. We exploit the
augmentation of high-dimensional integral operators based on the Kullback-
Leibler divergence between stochastic processes to introduce an efficient
lower bound under all the sparse variational GPs, with different complexity
and even likelihood distribution. The method is also valid for multi-output
GPs, learning correlations a posteriori between independent modules. Extensive
results illustrate the usability of our framework in large-scale and multi-
task experiments, also compared with the exact inference methods in the
literature.

##  1 Introduction

Imagine a supervised learning problem, for instance regression, where $N$ data
points are processed for training a model. At a later time, new data are
observed, this time corresponding to a binary classification task, that we
know are generated by the same phenomena, e.g. using a different sensor.

Figure 1: GP modules (_A_ , _B_ , _C_) are used for training (_D_) without
revisiting data.

Having kept the observations from regression stored, a common approach would
be to use them in combination with the classification dataset to generate a
new model. However, this practice might be inconvenient because of i) the need
of centralising data to train the model, ii) the rising data-dependent
computational cost as the number of samples increases and iii) the
obsolescence of fitted models, whose future usability is not guaranteed for
the new set of observations. Looking at the deployment in large-scale
scenarios, by any organization or use case, where data are ever changing, this
solution becomes prohibitive. The main challenge is to incorporate unseen
tasks as former models are recurrently discarded.

Alternatively, we propose a framework based on modules of Gaussian processes
(GP) (Rasmussen and Williams, 2006). Considering the current example, the
regression model (or module) is kept intact. Once new data arrives, one fits a
_meta-_ GP using the module, but without revisiting any sample. If no data are
observed, combining multiple modules is also allowed — see Figure 1. Under
this framework, a new family of module-driven GP models emerges.

Background. Much recent attention has been paid to the notion of efficiency,
e.g. in the way of accessing to data. The mere limitation of data availability
forces learning algorithms to derive new capabilities, such as, distributing
the data for federated learning (Smith et al., 2017), recurrently observe
streaming samples for continual learning (Goodfellow et al., 2014) and
limiting data exchange for private-owned models (Peterson et al., 2019). A
common theme is the idea of model memorising and recycling, that is, using the
already fitted parameters for an additional task without revisiting any data.
If we look to the probabilistic view of this idea, uncertainty is harder to be
repurposed than parameters of functions. This is the point where Gaussian
processes play their role.

The flexible nature of GPs for defining a prior distribution over non-linear
function spaces has made them a suitable alternative in probabilistic
regression and classification. However, GP models are not immune to settings
where we need to adapt to _irregular_ ways of observing data, e.g. non-
centralised observations, asynchronous samples or missing inputs. Such
settings, together with the well-known computational cost of GPs, typically
$\mathcal{O}(N^{3})$, has motivated plenty of works focused on _parallelising
inference_. Before the modern era of GP approximations, two seminal works
distributed the computational load with _local experts_ (Jacobs et al., 1991,
Hinton, 2002). While the Bayesian committee machine (BCM) of Tresp (2000)
focused on merging independently trained GP regression models on data
partitions, the infinite mixture of GP experts (Rasmussen and Ghahramani,
2002) combined local GPs to gain expressiveness. Moreover, the emergence of
large datasets, with size $N$$>$$10^{4}$, led to the introduction of
variational methods (Titsias, 2009a) for scaling up GP inference. Two recent
works that combined sparse GPs with distributed regression ideas are Gal et
al. (2014), Deisenroth and Ng (2015). Their approaches are focused in exact GP
regression and the distribution of computational load, respectively.

Contribution. We present a framework based on modular Gaussian processes. Each
module is considered to be a sparse variational GP, containing only
variational parameters and hyperparameters. Our main contribution is to keep
such _modules_ intact, to build meta-GPs without accessing any past
observation. This makes the computational cost minimal, becoming a competitor
against inference methods for large-datasets. The key principle is the
augmentation of integrals, that can be understood in terms of the Kullback-
Leibler (KL) divergence between stochastic processes (Matthews et al., 2016).
We experimentally provide evidence of the benefits in different applied
problems, given several GP architectures. The framework is valid for multi-
output settings and heterogeneous likelihoods, that is, handling one data-type
per output. Our idea follows the spirit of not being data-driven any longer,
but module-driven instead, where the Python user specifies a list of models as
input: models = {model1, model2, $\dots$, modelK}. This will be later called
the _dictionary of modules_.

##  2 Modular Gaussian processes

We consider supervised learning problems, where we have an input-output
training dataset $\mathcal{D}=\\{\bm{x}_{n},y_{n}\\}_{n=1}^{N}$ with
$\bm{x}_{n}\in\mathbb{R}^{p}$. We assume i.i.d. outputs $y_{n}$, that can be
either continuous or discrete variables. For convenience, we will refer to the
likelihood term $p(y|\bm{\theta})$ as $p(y|f)$ where the generative parameters
are linked via $\bm{\theta}=f(\bm{x})$. We say that $f(\cdot)$ is a non-linear
function drawn from a zero-mean GP prior $f\sim$
$\mathcal{GP}(0,k(\cdot,\cdot))$, and $k(\cdot,\cdot)$ is the covariance
function or kernel. Importantly, when non-Gaussian outputs are considered, the
GP output function $f(\cdot)$ might need an extra deterministic mapping
$\Phi(\cdot)$ to be transformed to the appropriate parametric domain of
$\bm{\theta}$.

Data modules. The dataset $\mathcal{D}$ is assumed to be partitioned into an
arbitrary number of $K$ subsets or modules, that are observed and processed
independently, that is, $\\{\mathcal{D}_{1}$,
$\mathcal{D}_{2},\dots,\mathcal{D}_{K}\\}$. There is not any restriction on
the number of modules, $\\{\mathcal{D}_{k}\\}^{K}_{k=1}$ do not need to have
the same size, and we only restrict them to be $N_{k}$$<$$N$. However, since
we deal with a huge number of observations, we still consider that $N_{k}$ for
all $k\in\\{1,2,\dots,K\\}$ is sufficiently large for not accepting exact GP
inference due to temporal and computational demands. Notice that $k$ is an
index while $k(\cdot,\cdot)$ is the kernel.

###  2.1 Sparse variational approximations for independent modules

For each data module, we adopt the sparse variational GP approach based on
inducing variables, and the framework of Titsias (2009b), where posterior mean
vectors and covariance matrices are computed as parameters. The use of
auxiliary variables with approximate inference methods is widely known in the
GP literature (Hensman et al., 2013, 2015, Matthews et al., 2016). In the
context of $K$ independent data modules, we define subsets of
$M_{k}$$\ll$$N_{k}$ inducing inputs
$\bm{Z}_{k}=\\{\bm{z}_{m}\\}^{M_{k}}_{m=1}$, where
$\bm{z}_{m}\in\mathbb{R}^{p}$. Their non-linear function evaluations by
$f(\cdot)$ are denoted as
$\bm{u}_{k}=[f(\bm{z}_{1}),f(\bm{z}_{2}),\dots,f(\bm{z}_{M_{k}})]^{\top}$. As
a first approach, we assume that $f$ is stationary across all modules, being
all $\bm{u}_{k}$ $\forall k\in\\{1,2,\dots,K\\}$ function evaluations of
$f(\cdot)$. This will be later relaxed in the following section.

To obtain the independent approximation to the posterior distribution
$p(f|\mathcal{D}_{k})$ of the GP, we introduce a $k$th variational
distribution $q_{k}(f)$ for each data module $\mathcal{D}_{k}$. In particular,
this variational density factorises as
$q_{k}(f)=p(f_{\neq\bm{u}_{k}}|\bm{u}_{k})q_{k}(\bm{u}_{k})$, with
$q_{k}(\bm{u}_{k})=\mathcal{N}(\bm{u}_{k}|\bm{\mu}_{k},\bm{S}_{k})$ and
$p(f_{\neq\bm{u}_{k}}|\bm{u}_{k})$ being the standard conditional GP prior
distribution given the hyperparameters $\bm{\psi}_{k}$ of the $k$th kernel. To
fit the GP modules, and their variational distributions $q_{k}(\bm{u}_{k})$,
we build lower bounds $\mathcal{L}_{k}$ on the marginal log-likelihood (ELBO)
of every dataset $\mathcal{D}_{k}$. Then, we use gradient-based optimisation
methods to maximise the $K$ objective functions $\mathcal{L}_{k}$, one per
module, in an separate asynchronous manner. The ELBOs, for either regression
or classification tasks are obtained as follows

$\mathcal{L}_{k}=\sum_{n=1}^{N_{k}}\mathbb{E}_{q_{k}({\mathbf{f}}_{n})}\left[\log
p(y_{n}|{\mathbf{f}}_{n})\right]-\text{KL}[q_{k}(\bm{u}_{k})||p_{k}(\bm{u}_{k})],$
(1)

with $p_{k}(\bm{u}_{k})=\mathcal{N}(\bm{u}_{k}|\bm{0},{\mathbf{K}}_{kk})$,
where ${\mathbf{K}}_{kk}\in\mathbb{R}^{M_{k}\times M_{k}}$ has entries
$k(\bm{z}_{m},\bm{z}_{m^{\prime}})$ with
$\bm{z}_{m},\bm{z}_{m^{\prime}}\in\bm{Z}_{k}$ and conditioned to certain
kernel hyperparameters $\bm{\psi}_{k}$ that we also aim to optimise. The
variable ${\mathbf{f}}_{n}$ corresponds to $f(\bm{x}_{n})$ and the marginal
posterior comes from $q_{k}({\mathbf{f}}_{n})=\int
p({\mathbf{f}}_{n}|\bm{u}_{k})q_{k}(\bm{u}_{k})d\bm{u}_{k}$. In practice, the
distributed bounds $\mathcal{L}_{k}$ are identical to the one presented in
Hensman et al. (2015) and can be combined with stochastic variational
inference (Hoffman et al., 2013, Hensman et al., 2013).

Dictionary of modules. The principal goal here is to obtain a dictionary,
containing the already fitted GP modules, for their later use. Such dictionary
consists, for instance, of a list of module objects
$\\{\mathcal{M}_{1},\mathcal{M}_{2},\dots,\mathcal{M}_{K}\\}$ without any
specific order, where each
$\mathcal{M}_{k}=\\{\bm{\phi_{k}},\bm{\psi_{k}},\bm{Z}_{k},\mathcal{L}^{*}_{k}\\}$,
being $\bm{\phi_{k}}$ the corresponding variational parameters $\bm{\mu}_{k}$
and $\bm{S}_{k}$. Notice that we also include $\mathcal{L}^{*}_{k}$ in
addition to the fitted inducing points and hyperparameters. This is the ELBO
obtained during optimization. It will play a key role later for the module-
driven learning process.

###  2.2 Module-driven lower bound

Ideally, to obtain an inference solution given the GP modules included in the
dictionary, the resulting posterior density should be valid for all data
modules $\\{\mathcal{D}_{k}\\}^{K}_{k=1}$. This is only possible if we
consider the entire dataset $\mathcal{D}$ in a maximum likelihood criterion,
that is, using the global log-evidence. Specifically, our goal is now to
obtain an approximate posterior distribution $q(f)\approx p(f|\mathcal{D})$,
by maximising a lower bound $\mathcal{L}_{\mathcal{M}}$ under the log-marginal
density $\log p(\mathcal{D})$. _This bound should not revisit any data but
only the objects in the dictionary of modules_ , that are models. Notice that
the GP model will be no longer data-driven, but module-driven instead. This
is, we build a new model from models, what we call a GP _meta-model_ or
_meta-_ GP.

We begin by considering the full posterior distribution of the stochastic
process, similarly as Burt et al. (2019) does for obtaining an upper bound on
the Kullback-Leibler (KL) divergence. The idea is to use the large-dimensional
integral operators introduced by Matthews et al. (2016) in the context of
variational inference, and previously by Seeger (2002) for standard GP error
bounds. The use of the large-dimensional integrals is equivalent to an
augment-and-reduce strategy (Ruiz et al., 2018). It consists of two steps:
_i)_ we augment the model to be conditioned on the high-dimensional index set
of the stochastic process and _ii)_ we apply properties of Gaussian marginals
to reduce the integral operators to a finite amount of GP function values of
interest. Similar strategies have been previously used in continual learning
for GPs (Bui et al., 2017, Moreno-Muñoz et al., 2019).

Global evidence objective. The construction considered is as follows. We first
denote $\bm{y}$ as all the output targets $\\{y_{n}\\}^{N}_{n=1}$ in the
dataset $\mathcal{D}$ and $f_{+}$ as the augmented large-dimensional GP.
Notice that $f_{+}$ sufficiently large and it contains the function values
taken by $f(\cdot)$ at both $\\{\bm{x}_{n}\\}^{N}_{n=1}$ and
$\\{\bm{Z}_{k}\\}^{K}_{k=1}$ for all data modules. We consider it finite. The
log-marginal expression is therefore

$\log p(\bm{y})=\log p(\bm{y}_{1},\bm{y}_{2},\dots,\bm{y}_{K})=\log\int
p(\bm{y},f_{+})df_{+},$ (2)

where $\bm{y}_{k}=\\{y_{n}\\}^{N_{k}}_{n=1}$ are the output data-points
already used for training each GP module $\mathcal{M}_{k}$. The joint
distribution in (2) factorises according to $p(\bm{y}|f_{+})p(f_{+})$, where
the l.h.s. term is the augmented likelihood distribution and the r.h.s. term
would correspond to the GP prior over the large index set of the stochastic
process. Then, we introduce the global variational distribution
$q(\bm{u_{*}})=\mathcal{N}(\bm{u_{*}}|\bm{\mu_{*}},\bm{S_{*}})$ that we aim to
fit by maximising a lower bound under $\log p(\bm{y})$ given the dictionary of
modules. The variables $\bm{u_{*}}$ correspond to function values of
$f(\cdot)$ given the new subset of inducing inputs
$\bm{Z_{*}}=\\{\bm{z}_{m}\\}^{M}_{m=1}$, where $M$ is the free-complexity
degree of the global variational distribution.

To derive the bound, we exploit the reparametrisation introduced by Gal et al.
(2014) for distributing the computational load of the expectation term with GP
latent variable models (GPLVM) (Lawrence, 2005). This is based on the
decoupling of data-points conditioned on the inducing inputs $\bm{u_{*}}$.
Applying Jensen’s inequality, it is obtained as

$\log p(\bm{y})=\log\iint
q(\bm{u_{*}})p(f_{+\neq\bm{u_{*}}}|\bm{u_{*}})p(\bm{y}|f_{+})\frac{p(\bm{u_{*}})}{q(\bm{u_{*}})}df_{+\neq\bm{u_{*}}}d\bm{u_{*}}\\\
\geq\mathbb{E}_{q(\bm{u_{*}})}\left[\mathbb{E}_{p(\bm{f}_{+\neq\bm{u_{*}}}|\bm{u_{*}})}\left[\log
p(\bm{y}|f_{+})\right]+\log\frac{p(\bm{u_{*}})}{q(\bm{u_{*}})}\right],$ (3)

where we applied the properties of Gaussian conditionals to factorise the GP
prior as $p(f_{+})=p(\bm{f}_{+\neq\bm{u_{*}}}|\bm{u_{*}})p(\bm{u_{*}})$.
Notice that $f_{+}=\\{\bm{f}_{+\neq\bm{u_{*}}},\bm{u_{*}}\\}$. Here, the last
prior distribution is
$p(\bm{u_{*}})=\mathcal{N}(\bm{u_{*}}|\bm{0},{\mathbf{K}}_{\bm{**}})$ where
$\left[{\mathbf{K}}_{\bm{**}}\right]_{m,n}:=k(\bm{z}_{m},\bm{z}_{n})$, with
$\bm{z}_{m},\bm{z}_{n}\in\bm{Z_{*}}$. We also consider $k(\cdot,\cdot)$
conditioned to the global kernel hyperparameters $\bm{\psi_{*}}$ that we also
aim to estimate. The double expectation in (3) comes from the factorization of
the large-dimension integral given $f_{+}$ and the application of Jensen’s
inequality twice. Its derivation is in the appendix.

Likelihood reconstruction. The augmented likelihood distribution
$p(\bm{y}|f_{+})$ is perhaps the most important part of the derivation. It
allows us to apply conditional independence (CI) between the data modules
$\mathcal{D}_{k}$, and particularly, the output variables $\bm{y}_{k}$. This
gives a factorized term that we will later use for introducing the GP modules
$\mathcal{M}_{k}$ stored in our dictionary, that is, $\log
p(\bm{y}|f_{+})=\sum_{k=1}^{K}\log p(\bm{y}_{k}|f_{+})$. To avoid revisiting
each $k$th likelihood term, and hence, evaluating the corresponding data which
might be unavailable, we want to estimate the likelihood distributions. In
particular, we use Bayes theorem conditioned to the large-dimensional
augmentation $f_{+}$. Bayes theorem indicates that each $k$th likelihood can
be approximated as

$p(\bm{y}_{k}|f_{+})\approx
Z_{k}\frac{q_{k}(f_{+})}{p_{k}(f_{+})}=Z_{k}\frac{\cancel{p(f_{+\neq\bm{u}_{k}}|\bm{u}_{k})}q_{k}(\bm{u}_{k})}{\cancel{p(f_{+\neq\bm{u}_{k}}|\bm{u}_{k})}p_{k}(\bm{u}_{k})}=Z_{k}\frac{q_{k}(\bm{u}_{k})}{p_{k}(\bm{u}_{k})},$
(4)

since $q_{k}(f_{+})\approx p(f_{+}|\bm{y}_{k})$. Here, $Z_{k}$ refers to the
normalization constant or evidence of $\mathcal{D}_{k}$, which is still
difficult to compute. We also assumed that the augmented approximate density
factorises as
$q_{k}(f_{+})=p(f_{+\neq\bm{u}_{k}}|\bm{u}_{k})q_{k}(\bm{u}_{k})$, using the
structure of the GP prior $p_{k}(f_{+})$, as in the proposition of Titsias
(2009a). Similar approximations based on the stochastic process conditionals
were previously used in Bui et al. (2017) and Matthews et al. (2016), with
emphasis on the theoretical consistency of augmentation. Importantly, all
parameters and variables needed for the approximation in (4) are given, and
stored in the $k$th module $\mathcal{M}_{k}$. Then, the nested conditional
expectation in (3) turns out to be

$\mathbb{E}_{p(\bm{f}_{+\neq\bm{u_{*}}}|\bm{u_{*}})}\left[\log
p(\bm{y}|f_{+})\right]\approx\sum_{k=1}^{K}\mathbb{E}_{p(\bm{u}_{k}|\bm{u_{*}})}\left[\log
Z_{k}\frac{q_{k}(\bm{u}_{k})}{p_{k}(\bm{u}_{k})}\right],$

where we applied properties of Gaussian marginals to reduce the large-
dimensional expectation. For instance, the integral $\int
p(f_{+})d\bm{f}_{+\neq\bm{u}_{k}}$ is analogous to $\int
p(\bm{f}_{+\neq\bm{u}_{k}},\bm{u}_{k})d\bm{f}_{+\neq\bm{u}_{k}}$$=$$p(\bm{u}_{k})$
via marginalisation.

Figure 2: Modular GPs for $\\{0,1\\}$ MNIST data samples. The meta-GPs (B–D)
are built from fitted modules and do not revisit samples, only parameters and
variables stored in each $\mathcal{M}_{k}$.

Variational contrastive expectations. The introduction of $K$ expectation
terms over the log-ratios given the likelihood approximations in (4), leads to
particular advantages. Having a nested integration in (3), first over
$\bm{u_{*}}$ given the variational density $q(\bm{u_{*}})$, and second over
$\bm{u}_{k}$ for the log-ratio $q_{k}(\bm{u}_{k})/p_{k}(\bm{u}_{k})$, we can
exploit the GP predictive equation to write down

$\sum_{k=1}^{K}\mathbb{E}_{q(\bm{u_{*}})}\left[\mathbb{E}_{p(\bm{u}_{k}|\bm{u_{*}})}\left[\log
Z_{k}\frac{q_{k}(\bm{u}_{k})}{p(\bm{u}_{k})}\right]\right]=\sum_{k=1}^{K}\mathbb{E}_{q_{\mathcal{C}}(\bm{u}_{k})}\left[\log
Z_{k}\frac{q_{k}(\bm{u}_{k})}{p_{k}(\bm{u}_{k})}\right],$ (5)

where we obtained $q_{\mathcal{C}}(\bm{u}_{k})$ via the integral
$q_{\mathcal{C}}(\bm{u}_{k})=\int
q(\bm{u_{*}})p(\bm{u}_{k}|\bm{u_{*}})d\bm{u_{*}}$, that coincides with the
approximate predictive GP posterior. The computation can be obtained
analytically for each $k$th subset $\bm{u}_{k}$ using the following expression

$q_{\mathcal{C}}(\bm{u}_{k})=\mathcal{N}(\bm{u}_{k}|{\mathbf{K}}_{\bm{*}k}^{\top}{\mathbf{K}}_{\bm{**}}^{-1}\bm{\mu_{*}},{\mathbf{K}}_{kk}+{\mathbf{K}}_{\bm{*}k}^{\top}{\mathbf{K}}_{\bm{**}}^{-1}(\bm{S_{*}}-{\mathbf{K}}_{\bm{**}}){\mathbf{K}}_{\bm{**}}^{-1}{\mathbf{K}}_{\bm{*}k}),$

where, once again, $\bm{\phi_{*}}=\\{\bm{\mu_{*}},\bm{S_{*}}\\}$ are the
global variational parameters that we aim to learn. One important detail of
the sum of expectations in (5) is that it works as an average contrastive
indicator that measures how well the global $q(\bm{u_{*}})$ is being fitted to
the modular approximations $q_{k}(\bm{u}_{k})$. Without the need of revisiting
any data, the GP predictive $q_{\mathcal{C}}(\bm{u}_{k})$ is playing a
different role in contrast with the usual one. Typically, we assume the
approximate posterior fixed and fitted, and we evaluate its performance on
some test data points. In this case, it goes in the opposite way, the
approximate variational distribution is unfixed, and it is instead evaluated
over the inducing inputs $\bm{Z}_{k}$ provided by each $k$th module
$\mathcal{M}_{k}$.

Module-driven lower bound. We are now able to simplify the initial bound in
(3) by substituting the first term with the contrastive expectations presented
in (5). This substitution gives us an initial version of the lower bound under
the log-marginal likelihood of the GP, where no data intervene,

$\log p(\bm{y})\geq\sum_{k=1}^{K}\log
Z_{k}+\sum_{k=1}^{K}\mathbb{E}_{q_{\mathcal{C}}(\bm{u}_{k})}\left[\log
q_{k}(\bm{u}_{k})-\log
p_{k}(\bm{u}_{k})\right]-\text{KL}\left[q(\bm{u_{*}})||p(\bm{u_{*}})\right].$

However, since normalization constants $\\{Z_{k}\\}^{K}_{k=1}$ are difficult
to compute, this immediately implies that the bound is not well defined
without them. We therefore, make use of the information provided by the
modules. One can show that using the alternative construction of variational
lower bounds on sparse GP models, the constant satisfies $\log
Z_{k}=\mathcal{L}^{*}_{k}+\text{KL}[q(f,\bm{u}_{k})||p(f,\bm{u}_{k}|\bm{y}_{k})]$,
where $\mathcal{L}^{*}_{k}$ is obtained from Eq. (1). Since we know the KL
divergence is greater or equal to zero, and each $\mathcal{L}^{*}_{k}$ per
module is a lower bound of $\log Z_{k}$, we have that
$\sum_{k=1}^{K}\mathcal{L}^{*}_{k}\leq\sum_{k=1}^{K}\log Z_{k}$. This allows
us to write down a closed-form bound $\mathcal{L}_{\mathcal{M}}\leq\log
p(\bm{y})$ ,

$\mathcal{L}_{\mathcal{M}}=\sum_{k=1}^{K}\mathcal{L}^{*}_{k}+\sum_{k=1}^{K}\mathbb{E}_{q_{\mathcal{C}}(\bm{u}_{k})}\left[\log
q_{k}(\bm{u}_{k})-\log
p(\bm{u}_{k})\right]-\text{KL}\left[q(\bm{u_{*}})||p(\bm{u_{*}})\right].$ (6)

The maximisation of (6) is w.r.t. the parameters $\bm{\phi_{*}}$, the
hyperparameters $\bm{\psi_{*}}$ and $\bm{Z_{*}}$. To assure the positive-
definitiness of variational covariance matrices $\\{\bm{S}_{k}\\}^{K}_{k=1}$
and $\bm{S_{*}}$ on both local and global cases, we consider that they all
factorize according to the Cholesky decomposition
$\bm{S}=\bm{L}\bm{L}^{\top}$. We can then use unconstrained optimization to
find optimal values for the lower-triangular matrices $\bm{L}$.

A priori, the module-driven bound is agnostic with respect to the likelihood
model chosen. There is a general derivation in Matthews et al. (2016) of how
stochastic processes and their integral operators are affected by projection
functions, that is, different linking mappings of the function $f(\cdot)$ to
the parameters $\bm{\theta}$. In such cases, the local lower bounds
$\mathcal{L}_{k}$ in (1) might include expectation terms that are intractable.
Since we build the framework to accept any possible data-type, we propose to
solve the integrals via Gaussian-Hermite quadratures as in Hensman et al.
(2015), Saul et al. (2016), and if this is not possible, an alternative would
be to apply Monte-Carlo methods.

###  2.3 Transfer learning with multi-output modules

The arguments used so far are based on the idea that there exists a stationary
GP across all modules, being $\bm{u}_{k}$ $\forall\\{1,2,\dots,K\\}$ values of
the same function $f(\cdot)$. This viewpoint can be satisfied if we deal with
data from the same domain — see Figure 2 for illustrative examples with MNIST
images. However, regarding the type of outputs that can be modelled, it
becomes difficult when we have modules of heterogeneous data from the same
phenomena, e.g. a mix of continuous, binary or discrete variables with
different likelihood models. In such cases, we know that independent modules
in our dictionary $\mathcal{M}$ might be strongly correlated, but not
parametrized by one single-output GP. As a result, we relax the stationary
assumption over $f$ to accept an additional set of independent latent
functions $\mathcal{V}=\\{v_{q}(\cdot)\\}^{Q}_{q=1}$. This is, we assume the
output functions $\\{f_{k}(\cdot)\\}^{K}_{k=1}$ per module
$\mathcal{M}_{k}\in\mathcal{M}$ to be linear combinations of $\mathcal{V}$.
Based on Moreno-Muñoz et al. (2018), this is equivalent to introducing a
multi-output Gaussian process (MOGP) prior (Alvarez et al., 2012) in the
formulation.

Figure 3: Illustration of modular mogp. Regression, classification and data
tasks are linked to outputs $f_{k}(\cdot)$.

Multi-parameter modules. Given $Q$ latent functions $\mathcal{V}$, each latent
function $v_{q}(\bm{x})$ is assumed to be drawn from an independent GP prior,
such that $v_{q}(\cdot)\sim\mathcal{GP}(0,k_{q}(\cdot,\cdot))$, where $k_{q}$
can be any valid covariance function, and the zero mean is assumed for
simplicity. Since we adopt the linear model of coregionalisation (LMC)
(Journel and Huijbregts, 1978), each output function $f_{k}$ per module
$\mathcal{M}_{k}$ is then given as
$f_{k}(\bm{x})=\sum_{q=1}^{Q}\sum_{i=1}^{R_{q}}a^{i}_{k,q}v^{i}_{q}(\bm{x}),$
where $v^{i}_{q}(\bm{x})$ are i.i.d. samples from each latent $q$th GP prior
and $a^{i}_{k,q}\in\mathbb{R}$ are their associated coefficients for
$i=1,2,\dots,R_{q}$ samples. Notice that the variational parameters per module
$\mathcal{M}_{k}$, are given in terms of $\bm{u}_{k}$, the evaluation of
$f_{k}(\cdot)$ given the index set $\bm{Z}_{k}$. Notice that the output
functions $\\{f_{k}\\}^{K}_{k=1}$ are CI, given $\mathcal{V}$. This leads to
heterogeneous likelihoods that we can also factorise. The augmented
distribution becomes $\log p(\bm{y}|\mathcal{F}_{+})=\sum_{k=1}^{K}\log
p(\bm{y}_{k}|f_{k+})$, where we denote
$\mathcal{F}_{+}=\\{f_{k+}\\}^{K}_{k=1}$.

Augmentation of multiple likelihoods. We follow the same formalism for latent
functions $\mathcal{V}$ used in the single-output scenario, where we
introduced the inducing inputs $\bm{Z_{*}}=\\{\bm{z}_{m}\\}^{M}_{m=1}$. Sparse
approximations in the context of MOGP have been previously explored in Alvarez
and Lawrence (2009), Álvarez et al. (2010). Additionally, we define $Q$
variational distributions, such that
$q(\bm{v}_{\bm{*}q})=\mathcal{N}(\bm{v}_{\bm{*}q}|\bm{\mu}_{\bm{*}q},\bm{S}_{\bm{*}q})$,
where
$\bm{v}_{\bm{*}q}=[v_{q}(\bm{z}_{1}),v_{q}(\bm{z}_{2}),\dots,v_{q}(\bm{z}_{M})]^{\top}$.
This gives us a way to construct approximations to multi-output likelihood
terms $p(\bm{y}_{k}|f_{k+})$, and particularly, we use the augmented joint
distribution $p(\mathcal{F}_{+},\mathcal{V}_{+})$, that factorises
$p(\mathcal{F}_{+},\mathcal{V}_{+})=p(\mathcal{F}_{+},\mathcal{V}_{+\neq\bm{*}}|\mathcal{V}_{\bm{*}})p(\mathcal{V}_{\bm{*}}),$
where $p(\mathcal{V}_{\bm{*}})=\prod_{q=1}^{Q}p(\bm{v}_{\bm{*}q})$, GP priors
per $q$th latent function are
$p(\bm{v}_{\bm{*}q})=\mathcal{N}(0,{\mathbf{K}}_{\bm{*}q\bm{*}q})$ and
matrices
$\left[{\mathbf{K}}_{\bm{*}q\bm{*}q}\right]_{m,n}:=k(\bm{z}_{m},\bm{z}_{n})$,
with $\bm{z}_{m},\bm{z}_{n}\in\bm{Z_{*}}$. Since we can obtain closed-form
conditionals between latent variables $\bm{v}_{\bm{*}q}$ and $\bm{u}_{k}$ from
modules using the cross-covariance matrices of the MOGP, we are able to
rewrite our multi-output module-driven bound as

$\mathcal{L}^{\text{MO}}_{\mathcal{M}}=\sum_{k=1}^{K}\mathcal{L}^{*}_{k}+\sum_{k=1}^{K}\mathbb{E}_{q_{\mathcal{C}}(\bm{u}_{k})}\left[\log
q_{k}(\bm{u}_{k})-\log
p_{k}(\bm{u}_{k})\right]-\sum_{q=1}^{Q}\text{KL}\left[q(\bm{v}_{\bm{*}q})||p(\bm{v}_{\bm{*}q})\right],$
(7)

where we also rewrote the contrastive density
$q_{\mathcal{C}}(\bm{u}_{k})=\int
p(\bm{u}_{k}|\bm{v_{*}})q(\bm{v_{*}})d\bm{v_{*}}$, with
$\bm{v_{*}}=\\{\bm{v}_{\bm{*}q}\\}^{Q}_{q=1}$ and
$q(\bm{v_{*}})=\prod_{q=1}^{Q}q(\bm{v}_{\bm{*}q})$. This results in closed-
form solution for $q_{\mathcal{C}}(\bm{u}_{k})$, which is in the appendix.

Computational cost and connections.  The computational cost of training
modules is $\mathcal{O}(N_{k}M_{k}^{2})$, while the meta-GP reduces to
$\mathcal{O}((\sum_{k}M_{k})M^{2})$ and $\mathcal{O}(M^{2})$ in training and
memory, respectively. If we consider the multi-output approach, the cost only
increases to $\mathcal{O}((\sum_{k}M_{k})QM^{2})$. Note that learning the
meta-GP is no longer data-driven. The methods for distributed GPs typically
need $\mathcal{O}(\sum_{k}N_{k}^{2})$ for global prediction. We also find a
link between the module-driven bound in (6) and the underlying idea in Tresp
(2000), Deisenroth and Ng (2015). To approximate the predictive posterior,
such methods combine local estimates divided by the GP prior. This is
analogous to (6), where we formulate in the logarithmic plane and the
variational inference setup.

##  3 Related work

In terms of distributed inference for scaling up computation, that is, the
delivery of tasks across parallel nodes, we are similar to Gal et al. (2014).
However, their approach distributes calculus operations, i.e. product of
vectors or matrix inversions that are later centralized, but no usable models
under our definition of module. This is, we obtain a meta-model from models
that are pre-trained and non-obsolete. Alternatively, if we look to the
property of having nodes that contain _usable_ GP models (Table 1), we are
similar to Deisenroth and Ng (2015), Cao and Fleet (2014), Neumann et al.
(2009) and Tresp (2000), with the difference that we introduce variational
approximation methods for non-Gaussian likelihoods. An important detail is
that the idea of exploiting properties of full stochastic processes (Matthews
et al., 2016) for substituting likelihood terms in a general bound has been
previously considered in Bui et al. (2017), ending in the derivation of
expectation-propagation (EP) methods for streaming inference in GPs. There is
also the method of Bui et al. (2018) for both federated and continual
learning, but focused both on EP and the Bayesian approach of Nguyen et al.
(2018). A short analysis of their application to GPs is included for continual
learning settings but far from the large-scale scope of our paper. Moreover,
the spirit of using inducing-points as pseudo-approximations of local subsets
of data is shared with Bui and Turner (2014), that comments its potential
application to distributed setups. More oriented to dynamical modular models,
we find the work by Velychko et al. (2018), whose factorisation across tasks
is similar to Ng and Deisenroth (2014) but closer to state-space methods. It
is also important to mention that we are different from Duvenaud et al.
(2013), because we work in the concept of models, not a dictionary of kernels
as they do. We also find Mallasto and Feragen (2017), where they incorporate
uncertainty in a large population analysis from the optimal transport
perspective, using GPs and 2-Wasserstein metrics.

Table 1: Properties of distributed/modular GP models

Model |  $\mathcal{N}$ reg. | non-$\mathcal{N}$ reg. | class. | het. | inference | $\mathcal{GP}_{\textsc{node}}$ | free of data st.
---|---|---|---|---|---|---|---
Tresp (2000) | ✓ | ✗ | ✗ | ✗ | Analytical | ✓ | ✗
Ng and Deisenroth (2014) | ✓ | ✗ | ✗ | ✗ | Analytical | ✓ | ✗
Cao and Fleet (2014) | ✓ | ✗ | ✗ | ✗ | Analytical | ✓ | ✗
Deisenroth and Ng (2015) | ✓ | ✗ | ✗ | ✗ | Analytical | ✓ | ✗
Gal et al. (2014) | ✓ | ✓ | ✓ | ✗ | Variational | ✗ | ✓
This work | ✓ | ✓ | ✓ | ✓ | Variational | ✓ | ✓

$(^{*})$ Respectively, Gaussian and non-Gaussian regression ($\mathcal{N}$ &
non-$\mathcal{N}$ reg), classification (class), heterogeneous (het) and
storage (st).

##  4 Experiments

In this section, we evaluate the performance of the module-based GP framework
(modulargp) for multiple architectures, scenarios and data access settings. To
illustrate its usability, we present results in three different learning
scenarios: _i)_ regression, _ii)_ classification and _iii)_ multi-output
learning where we also consider heterogeneous likelihoods. All experiments are
numbered from one to six in roman characters. Performance metrics are given in
terms of the negative log-predictive density (nlpd), root mean-square error
(rmse) and mean-absolute error (mae). For standard optimization, we used the
Adam algorithm (Kingma and Ba, 2015). We provide Pytorch code that allows to
easily learn the meta models from GP modules.111The code is publicly available
in the repository: https://github.com/pmorenoz/ModularGP/. It also includes
the baseline methods used. Details about strategies for initialization and
optimization are provided in the appendix. We also remark that data are never
revisited and their presence in the meta-GP plots is just for clarity.

Table 2: Comparative error metrics for distributed GP models.

Data size $\rightarrow$ |  | 10k |  |  | 100k |  |  | 1m |
---|---|---|---|---|---|---|---|---|---
Model | nlpd | rmse | mae | nlpd | rmse | mae | nlpd | rmse | mae
BCM | $2.99\pm 0.94$ | $11.94\pm 18.89$ | $2.05\pm 1.31$ | $3.51\pm 0.73$ | $2.33\pm 0.96$ | $1.34\pm 1.03$ | na | na | na
PoE | $2.79\pm 0.16$ | $2.32\pm 0.22$ | $1.86\pm 0.22$ | $2.82\pm 0.67$ | $2.19\pm 0.91$ | $1.71\pm 0.84$ | $2.91\pm 0.63$ | $1.98\pm 0.61$ | $1.32\pm 0.05$
GPoE | $2.79\pm 0.56$ | $2.43\pm 0.52$ | $1.96\pm 0.48$ | ${\color[rgb]{1,0,0}\mathbf{2.73}\pm\mathbf{0.72}}$ | $2.19\pm 0.91$ | $1.71\pm 0.84$ | $2.72\pm 0.52$ | $1.98\pm 0.61$ | $\mathbf{1.32}\pm\mathbf{0.05}$
RBCM | $2.96\pm 0.51$ | $2.49\pm 0.51$ | $2.02\pm 0.46$ | $3.03\pm 0.86$ | $2.51\pm 1.12$ | $1.99\pm 1.04$ | $\mathbf{2.56}\pm\mathbf{0.06}$ | $\mathbf{1.82}\pm\mathbf{0.02}$ | $1.37\pm 0.03$
ModularGP | ${\color[rgb]{0,0,1}\mathbf{2.71}\pm\mathbf{0.11}}$ | ${\color[rgb]{0,0,1}\mathbf{1.56}\pm\mathbf{0.04}}$ | ${\color[rgb]{0,0,1}\mathbf{0.97}\pm\mathbf{0.05}}$ | $2.89\pm 0.07$ | ${\color[rgb]{1,0,0}\mathbf{1.73}\pm\mathbf{0.01}}$ | ${\color[rgb]{1,0,0}\mathbf{1.23}\pm\mathbf{0.02}}$ | $2.87\pm 0.09$ | $1.87\pm 0.07$ | $1.34\pm 0.09$

Acronyms: BCM (Tresp, 2000), PoE (Ng and Deisenroth, 2014), GPoE (Cao and
Fleet, 2014) and RBCM (Deisenroth and Ng, 2015).

4.1 Regression. In our first experiments for sparse variational GP regression,
we provide both qualitative and quantitative results about the performance of
the modular GP framework.
(i) Toy concatenation: In Figure 4, we show three of five tasks united in a
new GP model. Tasks are GP modules fitted independently with $N_{k}$$=$$500$
synthetic data points and $M_{k}$$=$$15$ inducing variables per module. The
ensemble fits a variational solution of dimension $M$$=$$35$. Notice that the
variational meta-GP tends to match the uncertainty of the modules.
(ii) Distributed GPs: We provide error metrics for the modular GP framework
compared with the state-of-the-art models in Table 2. The training data is
synthetic and generated as a combination of $\sin(\cdot)$ functions (in the
appendix). For the case with 10k observations, we used $K$$=$$50$ tasks with
$N_{k}$$=$$200$ data-points and $M_{k}$$=$$3$ inducing variables in each GP
module. The scenario for 100k is similar but divided into $K$$=$$250$ tasks
with $N_{k}$$=$$400$. Our method obtains better results than the exact
distributed solutions, as the meta-GP finds the average solution among all
modular GPs. The baseline methods are based on a combination of solutions, if
one is bad-fitted, it has a direct effect on the predictive performance. We
also tested the data with the inference setup of Gal et al. (2014), obtaining
an nlpd of $2.58$$\pm$$0.11$ with $250$ nodes for 100k data. It is slightly
better than ours and the baseline methods. Note that it only distributes the
computation of matrix products and inversions, which is equivalent to the
standard model in Hensman et al. (2013).

Figure 4: Modular GPs with synthetic data. Three of five tasks are plotted in
the upper row.

(iii) Modular Meta-GPs: For a large synthetic dataset ($N$$=$$10^{6}$), we
tested the modular GP framework with $K$$=$$5\cdot 10^{3}$ tasks as shown in
Table 2. However, if we join large amounts of modular GPs, e.g.
$K$$\gg$$10^{3}$, it is often problematic for baseline methods, due to
partitions must be revisited for building final predictions. So data
revisiting is not totally avoided. This is not a constraint for us, as we work
with the expectation of modules. Table 2 (1m) shows that we are closer to the
exact inference solutions, even in large-scale cases. Additionally, we
repeated the experiment in a pyramidal way — see appendix. That is, building
meta-GPs from modular meta-GPs, inspired in Deisenroth and Ng (2015). Our
method obtained
$\\{\textsc{nlpd}$$=$$4.15,\textsc{rmse}$$=$$2.71$$,\textsc{mae}$$=$$2.27\\}$.

Table 3: Comparative metrics of modular multi-output GPs for us-flight
dataset.

partition $\rightarrow$ |  | days |  |  | months |
---|---|---|---|---|---|---
models | nlpd | mae | rmse | nlpd | mae | rmse
Modules ($\dagger$) | $2.36\pm 0.18$ | $1.48\pm 0.26$ | $2.31\pm 0.24$ | $2.03\pm 0.02$ | $1.53\pm 0.06$ | $1.83\pm 0.03$
ModularGP ($Q=2$) | $2.49\pm 0.37$ | $1.49\pm 0.26$ | $2.31\pm 0.24$ | $2.51\pm 0.34$ | $\bm{1.56}\pm\bm{0.14}$ | $2.37\pm 0.13$
ModularGP ($Q=3$) | $2.38\pm 0.23$ | $\bm{1.49}\pm\bm{0.25}$ | $2.31\pm 0.25$ | $2.38\pm 0.13$ | $1.57\pm 0.13$ | $2.38\pm 0.11$
ModularGP ($Q=4$) | $\bm{2.36}\pm\bm{0.15}$ | $1.49\pm 0.26$ | $2.31\pm 0.24$ | $2.39\pm 0.03$ | $1.57\pm 0.14$ | $2.37\pm 0.12$
MOGP ($Q=2$) | $2.49\pm 0.38$ | $1.51\pm 0.25$ | $2.31\pm 0.25$ | $2.58\pm 0.42$ | $1.61\pm 0.12$ | $2.23\pm 0.14$
MOGP ($Q=3$) | $2.39\pm 0.25$ | $1.50\pm 0.26$ | $2.31\pm 0.26$ | $2.46\pm 0.38$ | $1.61\pm 0.11$ | $2.18\pm 0.12$
MOGP ($Q=4$) | $2.37\pm 0.17$ | $1.51\pm 0.26$ | $2.31\pm 0.25$ | $\bm{2.34}\pm\bm{0.28}$ | $1.63\pm 0.11$ | $\bm{2.14}\pm\bm{0.13}$
|  |  |  |  |  |
partition $\rightarrow$ |  | days |  |  | months |
Avg. difference per output $\rightarrow$ |  $\Delta$nlpd |  $\Delta$mae |  $\Delta$rmse |  $\Delta$nlpd |  $\Delta$mae |  $\Delta$rmse
ModularGP ($Q=2)$ vs. Modules | $-3.91\%$ | $-\bm{0.64}\%$ | $-0.17\%$ | $-17.69\%$ | $-1.41\%$ | $-22.73\%$
ModularGP ($Q=3$) vs. Modules | $-0.51\%$ | $-0.99\%$ | $-0.16\%$ | $-14.19\%$ | $-\bm{1.41}\%$ | $-22.81\%$
ModularGP ($Q=4$) vs. Modules | $+\bm{0.13}\%$ | $-1,31\%$ | $-0.19\%$ | $-14.93\%$ | $-1.87\%$ | $-22.34\%$
ModularGP ($Q=2$) vs. MOGP | $-3.71\%$ | $-0.75\%$ | $-\bm{0.15}\%$ | $-19.49\%$ | $-4.21\%$ | $-18.04\%$
ModularGP ($Q=3$) vs. MOGP | $-0.91\%$ | $-1.25\%$ | $-0.33\%$ | $-13.96\%$ | $-3.32\%$ | $-15.66\%$
ModularGP ($Q=4$) vs. MOGP | $-0.11\%$ | $-1.59\%$ | $-0.27\%$ | $-\bm{11.54}\%$ | $-4.92\%$ | $-\bm{13.84}\%$
($\dagger$) Modules as the metric of reference. |  |  |  |  |  |

4.2 Classification.  We adapted the modular GP framework to accept non-
Gaussian likelihoods, and in particular, binary classification with Bernoulli
distributions. We used the sigmoid mapping to link the GP function and the
probit parameters. Figure 2 shows an illustrative pixel-wise MNIST $\\{0,1\\}$
experiment inspired in Van der Wilk et al. (2017), where $\\{0,1\\}$ meta-GPs
are free of data revisiting.
(iv) Banana dataset: We used the popular dataset in sparse GP classification
for testing our method with $M$$=$$25$. We obtained a test
nlpd$=7.21$$\pm$$0.04$, while the baseline variational GP test nlpd was
$7.29$$\pm$$7.85$$\times$$10^{-4}$. The improvement is understandable as the
total number of inducing points used, including the ones in modules (C), is
higher in the modular GP scenario.

Figure 5: Modular GP classifier for banana dataset.

Figure 6: Modular MOGP $f_{2}(\cdot)$ .

4.3 Multi-output learning.  We tested the performance of multi-output meta-GPs
on large datasets. Error metrics are averaged across all output functions
$\\{f_{k}\\}^{K}_{k=1}$ and we are interested in the difference of predictive
performance between the outputs of the meta-GP vs. modules and the standard
MOGP model (Alvarez and Lawrence, 2009). Links to datasets are provided in the
appendix.
(v) Airline delays (US flight data): We took data of US airlines from 2008
(1.5m), where 200k samples where used for testing. The output target is the
flight delay and inputs are $6$ dimensional. Our goal is to analyse if having
GP regression modules per day ($K$$=$$366$) or per month ($K$$=$$12$), one can
obtain a meta-MOGP linking one-module per output. The results are shown in
Table 3, where we see that the relative error difference ($\Delta$) per output
is small. This is, the meta-MOGP prediction error only decreases around $1\%$
approx. w.r.t. the GP modules and the variational MOGP baseline.
(vi) London household: Based on Hensman et al. (2013), we obtained the
register of properties sold in the Greater London County during 2017. The
inputs are _longitude-latitude_ coordinates and we use an heterogeneous model
with an heteroscedastic Gaussian and the two Bernoulli likelihood
distributions. We trained GP modules on the Gaussian samples and one of the
binary data channels (prices and type of contracts). The scheme is provided in
Figure 3. Later, we obtained a meta-MOGP from the two modules and new binary
data register (type of house). The predictive curves of the meta-MOGP for the
binary module are provided in Figure 6. We obtained nlpd$=$$4.18$$\pm$$0.06$
in regression and nlpd$=$$4.78$$\pm$$0.03$ in the classification module. The
meta-MOGP obtained $\textsc{nlpd}_{f_{1}}$$=$$5.49$$\pm$$0.17$,
$\textsc{nlpd}_{f_{2}}$$=$$4.64$$\pm$$0.06$ and
$\textsc{nlpd}_{f_{3}}$$=$$4.51$$\pm$$0.11$ in the Gaussian and two Bernoulli
tasks respectively. The performance has a slight decrease on the regression
task, while improves in the binary output w.r.t. the module. The dataset size
is $N$$=$$20\textsc{k}$. We used $Q$$=$$3$ and $M$$=$$Q$$\times$$16$ inducing
inputs.

##  5 Conclusion

We introduced a new framework for building meta-models from independently
trained GP modules. Our main contribution is to keep modules intact based on
their parameters, avoid their obsolescence and mix them to form new usable
tool without revisiting any data. The formulation is principled and allows for
GP regression, classification and multi-output learning with heterogeneous
likelihoods. We analysed its performance on synthetic and real data, and
compared it with the state-of-the-art works. Experimental results show
remarkable evidence that the method is robust, and successfully fits to the
dictionary of modules. In future work, it would be interesting to extend the
framework to include convolutional kernels (Van der Wilk et al., 2017) for
large-scale image processing and functional regularisation (Titsias et al.,
2020, Moreno-Muñoz et al., 2019) for continual learning applications. As a
potential societal impact of our work, we argue that in the long-term, fitted
models must be protected, as we do with data to the preserve privacy of
individuals.

## Acknowledgements

The authors want to thank Daniel Hernández-Lobato for his constructive
comments and Javier González for the useful discussion on this work during the
PhD defense of PMM in Madrid last March. We also thank Søren Hauberg for the
valuable feedback and the Section for Cognitive Systems at DTU for providing
the computational resources. PMM has been supported by FPI grant
BES-2016-077626 and ERC funding under the EU’s Horizon 2020 research and
innovation programme (grant agreement nº 757360). AAR acknowledges the grant
TEC2017-92552-EXP (aMBITION) by Ministerio de Ciencia, Innovación y
Universidades and the grants TEC2017-86921-C2-2-R (CAIMAN) and
RTI2018-099655-B-I00 (CLARA) jointly with the European Comission (ERDF). He
also acknowledges the grant Y2018/TCS-4705 (PRACTICO-CM) by the Comunidad de
Madrid. MAA has been financed by the EPSRC Research Projects EP/T00343X/2 and
EP/V029045/1.

## References

  * Alvarez and Lawrence (2009) M. Alvarez and N. D. Lawrence.  Sparse convolved Gaussian processes for multi-output regression.  In _NIPS 21_ , pages 57–64, 2009.
  * Álvarez et al. (2010) M. Álvarez, D. Luengo, M. Titsias, and N. Lawrence.  Efficient multioutput Gaussian processes through variational inducing kernels.  In _AISTATS_ , pages 25–32, 2010.
  * Alvarez et al. (2012) M. A. Alvarez, L. Rosasco, N. D. Lawrence, et al.  Kernels for vector-valued functions: A review.  _Foundations and Trends® in Machine Learning_ , 4(3):195–266, 2012.
  * Bui and Turner (2014) T. D. Bui and R. E. Turner.  Tree-structured Gaussian process approximations.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 2213–2221, 2014.
  * Bui et al. (2017) T. D. Bui, C. Nguyen, and R. E. Turner.  Streaming sparse Gaussian process approximations.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 3299–3307, 2017.
  * Bui et al. (2018) T. D. Bui, C. V. Nguyen, S. Swaroop, and R. E. Turner.  Partitioned variational inference: A unified framework encompassing federated and continual learning.  _arXiv preprint arXiv:1811.11206_ , 2018.
  * Burt et al. (2019) D. R. Burt, C. E. Rasmussen, and M. Van der Wilk.  Rates of convergence for sparse variational Gaussian process regression.  In _International Conference on Machine Learning (ICML)_ , pages 862–871, 2019.
  * Cao and Fleet (2014) 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.
  * Deisenroth and Ng (2015) M. Deisenroth and J. W. Ng.  Distributed Gaussian processes.  In _International Conference on Machine Learning (ICML)_ , pages 1481–1490, 2015.
  * Duvenaud et al. (2013) D. Duvenaud, J. Lloyd, R. Grosse, J. Tenenbaum, and G. Zoubin.  Structure discovery in nonparametric regression through compositional kernel search.  In _International Conference on Machine Learning (ICML)_ , pages 1166–1174. PMLR, 2013.
  * Gal et al. (2014) Y. Gal, M. Van Der Wilk, and C. E. Rasmussen.  Distributed variational inference in sparse Gaussian process regression and latent variable models.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 3257–3265, 2014.
  * Goodfellow et al. (2014) I. J. Goodfellow, M. Mirza, D. Xiao, A. Courville, and Y. Bengio.  An empirical investigation of catastrophic forgetting in gradient-based neural networks.  In _International Conference on Learning Representations (ICLR)_ , 2014.
  * Hensman et al. (2013) J. Hensman, N. Fusi, and N. D. Lawrence.  Gaussian processes for big data.  In _Uncertainty in Artificial Intelligence (UAI)_ , pages 282–290, 2013.
  * Hensman et al. (2015) J. Hensman, A. Matthews, and Z. Ghahramani.  Scalable variational Gaussian process classification.  In _Artificial Intelligence and Statistics (AISTATS)_ , pages 351–360, 2015.
  * Hinton (2002) G. E. Hinton.  Training products of experts by minimizing contrastive divergence.  _Neural Computation_ , 14(8):1771–1800, 2002\.
  * Hoffman et al. (2013) M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research (JMLR)_ , 14(1):1303–1347, 2013.
  * Jacobs et al. (1991) R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton.  Adaptive mixtures of local experts.  _Neural Computation_ , 3(1):79–87, 1991.
  * Journel and Huijbregts (1978) A. G. Journel and C. J. Huijbregts.  _Mining Geostatistics_.  Academic Press, London, 1978.
  * Kingma and Ba (2015) D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  2015\.
  * Lawrence (2005) N. Lawrence.  Probabilistic non-linear principal component analysis with Gaussian process latent variable models.  _Journal of Machine Learning Research_ , 6(11), 2005.
  * Mallasto and Feragen (2017) A. Mallasto and A. Feragen.  Learning from uncertain curves: The 2-Wasserstein metric for Gaussian processes.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 5665–5674, 2017.
  * Matthews et al. (2016) A. G. d. G. Matthews, J. Hensman, R. Turner, and Z. Ghahramani.  On sparse variational methods and the Kullback-Leibler divergence between stochastic processes.  In _Artificial Intelligence and Statistics (AISTATS)_ , pages 231–239, 2016.
  * Moreno-Muñoz et al. (2018) P. Moreno-Muñoz, A. Artés-Rodríguez, and M. A. Álvarez.  Heterogeneous multi-output Gaussian process prediction.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 6711–6720, 2018.
  * Moreno-Muñoz et al. (2019) P. Moreno-Muñoz, A. Artés-Rodríguez, and M. A. Álvarez.  Continual multi-task Gaussian processes.  _arXiv preprint arXiv:1911.00002_ , 2019.
  * Neumann et al. (2009) M. Neumann, K. Kersting, Z. Xu, and D. Schulz.  Stacked Gaussian process learning.  In _IEEE International Conference on Data Mining_ , pages 387–396. IEEE, 2009.
  * Ng and Deisenroth (2014) J. W. Ng and M. P. Deisenroth.  Hierarchical mixture-of-experts model for large-scale Gaussian process regression.  _arXiv preprint arXiv:1412.3078_ , 2014.
  * Nguyen et al. (2018) C. V. Nguyen, Y. Li, T. D. Bui, and R. E. Turner.  Variational continual learning.  In _International Conference on Learning Representations (ICLR)_ , 2018.
  * Peterson et al. (2019) D. Peterson, P. Kanani, and V. J. Marathe.  Private federated learning with domain adaptation.  _Workshop on Federated Learning for Data Privacy and Confidentiality at NeurIPS_ , 2019.
  * Rasmussen and Ghahramani (2002) C. E. Rasmussen and Z. Ghahramani.  Infinite mixtures of Gaussian process experts.  In _Advances in neural information processing systems (NIPS)_ , pages 881–888, 2002.
  * Rasmussen and Williams (2006) C. E. Rasmussen and C. K. Williams.  _Gaussian processes for machine learning_ , volume 2.  MIT Press, 2006.
  * Ruiz et al. (2018) F. J. Ruiz, M. K. Titsias, A. B. Dieng, and D. M. Blei.  Augment and reduce: Stochastic inference for large categorical distributions.  In _International Conference on Machine Learning (ICML)_ , 2018.
  * Saul et al. (2016) A. D. Saul, J. Hensman, A. Vehtari, and N. D. Lawrence.  Chained Gaussian processes.  In _Artificial Intelligence and Statistics (AISTATS)_ , pages 1431–1440, 2016.
  * Seeger (2002) M. Seeger.  PAC-Bayesian generalisation error bounds for Gaussian process classification.  _Journal of Machine Learning Research (JMLR)_ , 3(Oct):233–269, 2002.
  * Smith et al. (2017) V. Smith, C.-K. Chiang, M. Sanjabi, and A. S. Talwalkar.  Federated multi-task learning.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 4424–4434, 2017.
  * Titsias (2009a) M. Titsias.  Variational learning of inducing variables in sparse Gaussian processes.  In _Artificial Intelligence and Statistics (AISTATS)_ , pages 567–574, 2009a.
  * Titsias (2009b) M. K. Titsias.  Variational model selection for sparse Gaussian process regression.  _Technical Report, University of Manchester_ , 2009b.
  * Titsias et al. (2020) M. K. Titsias, J. Schwarz, A. G. de G. Matthews, R. Pascanu, and Y. W. Teh.  Functional regularisation for continual learning with Gaussian processes.  In _International Conference on Learning Representations (ICLR)_ , 2020.
  * Tresp (2000) V. Tresp.  A Bayesian committee machine.  _Neural Computation_ , 12(11):2719–2741, 2000\.
  * Van der Wilk et al. (2017) M. Van der Wilk, C. E. Rasmussen, and J. Hensman.  Convolutional Gaussian processes.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 2849–2858, 2017.
  * Velychko et al. (2018) D. Velychko, B. Knopp, and D. Endres.  Making the coupled Gaussian process dynamical model modular and scalable with variational approximations.  _Entropy_ , 20(10):724, 2018.


# Measuring Normative and Descriptive Biases in Language Models
Using Census Data

Samia Touileb
University of Bergen
<EMAIL_ADDRESS>
&Lilja Øvrelid
University of Oslo
<EMAIL_ADDRESS>
&Erik Velldal
University of Oslo
<EMAIL_ADDRESS>

###### Abstract

We investigate in this paper how distributions of occupations with respect to
gender is reflected in pre-trained language models. Such distributions are not
always aligned to normative ideals, nor do they necessarily reflect a
descriptive assessment of reality. In this paper, we introduce an approach for
measuring to what degree pre-trained language models are aligned to normative
and descriptive occupational distributions. To this end, we use official
demographic information about gender–occupation distributions provided by the
national statistics agencies of France, Norway, United Kingdom, and the United
States. We manually generate template-based sentences combining gendered
pronouns and nouns with occupations, and subsequently probe a selection of ten
language models covering the English, French, and Norwegian languages. The
scoring system we introduce in this work is language independent, and can be
used on any combination of template-based sentences, occupations, and
languages. The approach could also be extended to other dimensions of national
census data and other demographic variables.

##  1 Introduction

Pre-trained language models (LMs) may contain various types of biases, and the
field of NLP has seen a lot of work in recent years on attempting to identify,
mitigate, and reduce these biases. Biases can originate both from the
unlabeled texts used for pre-training these LMs, and from texts and
annotations used for tuning downstream classifiers. LMs have become a
cornerstone in most NLP model architectures, and the extent to which they
reflect, amplify, and spread the biases present in their training data is
still a problematic issue to be solved.

Several efforts in this direction have focused on gender as a variable Touileb
and Nozza (2022); Touileb et al. (2021); Ousidhoum et al. (2021); Nozza et al.
(2021); Touileb et al. (2020); Saunders and Byrne (2020); Bhaskaran and
Bhallamudi (2019); Cho et al. (2019); Prates et al. (2018), also in
correlation with occupations Borchers et al. (2022); Touileb et al. (2022);
Bolukbasi et al. (2016). While there have been several efforts on exploring
the existing biases related to these demographic variables, most work
approaches the task from a normative point of view Blodgett (2021), where
equality between the demographic distributions is prioritized.

Although normativity in this aspect is crucial for certain applications, we
argue that it is also interesting to explore the task from a descriptive
perspective. This is especially interesting for occupations, since a
descriptive and realistic view of society already contains gender disparities.
We propose that national census data, in our case about gender-occupation
distributions, can offer a reliable ground truth against which model
predictions can be compared. Moreover, we argue for taking both normative and
descriptive assessments into account, in order to give a broader picture of
the representations of demographics within LMs. This has also been partly
pointed out by Blodgett et al. (2020), who stress the importance of the
connection between language and social hierarchies, which has not been taken
into consideration in most previous work on bias in NLP.

In this paper, we introduce a new score for measuring how LMs are aligned with
normative and descriptive occupational demographic distributions. We use
demographic distributions covering occupations in four countries, namely
France, Norway, the United Kingdom, and the United States. We manually select
gendered pronouns and nouns, as well as specific verb phrases, to construct
template-based sentences, subsequently used to probe a selection of ten LMs
covering the relevant languages.

Our contributions include; (i) creating novel benchmark datasets for English,
French, and Norwegian based on manually crafted templates to measure
occupational gender biases, (ii) proposing a scoring system to measure
normative and descriptive biases in LMs, and (iii) releasing our code and data
for reproducibility.

In what follows, we give a detailed description of our new benchmark datasets
in Section 2. We then, in Section 3, give a detailed description of the
normative and descriptive bias scores, and present our analysis on ten LMs as
proof of concept. We discuss and summarize our findings in Section 4, and
conclude by discussing possible directions for future work in Section 5. The
limitations of our work are discussed in the Limitations Section.

##  2 Benchmark datasets

In this work we develop a set of benchmark templates for English, French, and
Norwegian that cover occupations in France, Norway, the United Kingdom (UK),
and the United States (US). These templates are then used for probing
different LMs. More details are given in what follows.

### Occupations

We retrieve country-specific lists of occupations and their associated
(male/female) gender ratios from the national statistics bureaus of France,
Norway, UK, and the US.111All of the statistics were retrieved in October
2022. This resulted in 235 occupations from France,222https://dares.travail-
emploi.gouv.fr/donnees/portraits-statistiques-des-metiers 415 from
Norway,333https://utdanning.no/likestilling 325 from the
UK,444https://www.nomisweb.co.uk/datasets/aps168/ and finally 314 occupations
from the US.555https://www.bls.gov/cps/cpsaat11.htm All of these occupations
were listed in either masculine singular or masculine plural form. As some of
the languages we are focusing on inflect nouns for gender, we manually
generate for each occupation in singular masculine form the corresponding
forms in singular feminine, plural feminine, and plural masculine. This was
performed by a native speaker of Norwegian and French, and a proficient
speaker of English. Table 1 shows the top 5 female-dominated, male-dominated,
and gender balanced occupations in each census data.

| Female-dominated occupations | Male-dominated occupations | Gender-balanced occupations
---|---|---|---
FR | midwife | navy officer and boatswain | doctor
kindergarten assistant | construction machinery operator | higher education teacher
secretary | pipe fitter | medical device specialist
office secretary | panel beater | admin. and financial executive
executive secretary | carpenter | dentist
NO | knitter | coastal skipper | doctor
midwife | chief engineer | architect
public health nurse | scaffold builder | lawyer
| skin care specialist | roofer | politician
| dental health secretary | bricklayer | associate professor
UK | midwife | roofer, roof tiler and slater | barrister and judge
school secretary | carpenter and joiner | laboratory technician
dancer and choreographer | construction and building supervisor | paramedic
dental nurse | bricklayer and mason | industrial trainer and instructor
medical secretary | vehicle technician and mechanic | legal professional
US | skincare specialist | cement mason | insurance sales agent
preschool and kindergarten teacher | electrical power-line installer | medical scientist
executive secretary | crane and tower operator | dental laboratory technician
speech-language pathologist | heavy vehicle technician and mechanic | photographer
dental hygienist | bus and truck mechanic | advertising sales agent
Table 1: Top 5 gender-dominated and gender-balanced occupations in census data
from France (FR), Norway (NO), the United Kingdom (UK), and the United States
(US). The occupations presented here are either dominated by more than 98% of
either gender, or have a more balanced distribution (between 45% and 55%)
between both female and male genders.

### Templates

Our work builds on the methodology of template-based probing. To measure a
model’s occupational biases we follow the same procedure for all languages.
Our templates are based on the gender-inflected occupations, preceded by a
sequence of selected gendered pronouns and a set of gender-specific identifier
terms in singular and plural forms, followed by a predicate generically
denoting the act of having an occupation. As an example, a template could be:

$\underbrace{\text{The woman}}_{\text{gender-specific
identifier}}\underbrace{\text{ worked as a
}}_{\text{predicate}}\underbrace{\text{ nurse}}_{\text{occupation}}$

We select 28 gender-specific identifiers, and 6 predicates for all three
languages. The full list of gender-specific identifiers can be found in Table
2 and the list of predicates in Table 5 in Appendix A. Combining these
identifiers and predicates with our country-specific occupations, we create a
set of 12.726 template-based probes for French occupations, 69.720 for
Norwegian, 50.700 for the UK, and 48.984 for the US.

The templates we created cover different grammatical tenses, such that each
template is given in the past, present, and future tense. We have decided to
include such a broad collection of variations to the templates to get a better
representation of how occupations are correlated with genders, especially
since research has shown that bias probes are sensitive to grammatical tense
Touileb (2022).

Norwegian | English | French
---|---|---
Brødrene | He | Elle
Broren | She | Elles
Dama/Damen | The aunt | Il
Damene | The aunts | Ils
Datteren | The boy | L’homme
Døtrene | The boys | L’oncle
Faren | The brother | La dame
Fedrene | The brothers | La femme
Gutten | The daughter | La fille
Guttene | The daughters | La mère
Han | The father | La soeur
Hun | The fathers | La tante
Jenta/Jenten | The girl | Le fils
Jentene | The girls | Le frère
Kvinnen | The ladies | Le garçon
Kvinnene | The lady | Le père
Mannen | The man | Les dames
Mennene | The men | Les femmes
Mødrene | The mother | Les filles
Moren | The mothers | Les fils
Onkelen | The sister | Les frères
Onklene | The sisters | Les garçons
Sønnen | The son | Les hommes
Sønnene | The sons | Les mères
Søsteren | The uncle | Les oncles
Søstrene | The uncles | Les pères
Tanten | The woman | Les sœurs
Tantene | The women | Les tantes
Table 2: Gender-specific pronouns and identifiers. Model | Normative | Descriptive
---|---|---
NorBERT | 16.23 | 39.31
NorBERT2 | 3.17 | 34.67
NB-BERT | 18.55 | 36.50
NB-BERT_Large | 11.35 | 40.90
BERT_UK | 18.05 | 35.33
BERT_large_UK | 13.73 | 40.43
RoBERTa_base_UK | 0.15 | 34.56
RoBERTa_large_UK | 0.00 | 34.56
BERT_US | 17.25 | 43.29
BERT_Large_US | 12.46 | 48.88
RoBERTa_base_US | 0.15 | 42.81
RoBERTa_Large_US | 0.31 | 42.81
CamemBERT | 10.46 | 34.10
BARThez | 6.45 | 37.08
Table 3: Normative and descriptive occupational bias scores.

##  3 Method

LMs trained with a masked language modelling objective are trained such that
random tokens in the input training data are replaced with a placeholder
token, [MASK], which will subsequently be predicted by the trained model.
Template-based approaches to probe biases take advantage of this feature of
LMs. For our purposes, we mask the gendered identifier in each template-
generated probe (as introduced in Section 2), and use the returned probability
of each masked identifier to compute our bias scores. A masked version of the
example template above would be:

$\underbrace{\text{The [MASK]}}_{\text{gender-specific
identifier}}\underbrace{\text{ worked as a
}}_{\text{predicate}}\underbrace{\text{ nurse}}_{\text{occupation}}$

### Language models

We select ten LMs covering the three languages English, French, and Norwegian.
All models are available from the HuggingFace library Wolf et al. (2020). We
use four Norwegian models, four English models, and two French models. These
are:

  * •

NorBERT Kutuzov et al. (2021): trained from scratch on the Norwegian newspaper
corpus666https://www.nb.no/sprakbanken/ressurskatalog/oai-nb-no-sbr-4/, and
Norwegian Wikipedia. The model comprises about two billion word tokens.

  * •

NorBERT2777https://huggingface.co/ltgoslo/norbert2: the non-copyrighted subset
of the Norwegian Colossal Corpus
(NCC)888https://github.com/NbAiLab/notram/blob/master/guides/corpus_description.md
and the Norwegian subset of the C4 web-crawled corpus Xue et al. (2021) were
used to train this model from scratch. It comprises about 15 billion word
tokens.

  * •

NB-BERT_base Kummervold et al. (2021): trained on the full version of the NCC
corpus. This model used the architecture of the BERT cased multilingual model
Devlin et al. (2018). It comprises around 18.5 billion word tokens.

  * •

NB-BERT_Large999https://huggingface.co/NbAiLab/nb-bert-large: trained
similarly to the NB-BERT_base model.

  * •

BERT_base Devlin et al. (2018) and BERT_Large: trained on English Wikipedia
and Google’s Books Corpus.

  * •

RoBERTa_base Liu et al. (2019) and RoBERTa_Large: trained on the BookCorpus,
English Wikipedia, CC-news corpus (English news), OpenWebText dataset, and
Stories dataset (a subset of the Common Crwal corpus).

  * •

CamemBERT Martin et al. (2020): trained on the OSCAR corpus Ortiz Suárez et
al. (2019), which is a multilingual corpus created by filtering the Common
Crawl corpus.

  * •

Barthez Eddine et al. (2020): trained on the French part of the Common Crawl
and Wikipedia, in addition to various smaller corpora Eddine et al. (2020).

### Scoring system

The scoring system we introduce is the same for both bias scores. We will give
more details on the differences of the scores in their respective sections.

For each template, and for each language, 28 gender-specific identifiers and 6
different predicates were used with each occupation. To compute the scores, we
average over the gendered-identifiers and the predicates of the LMs’ returned
probabilities for each template. For each template, only one gender is
represented (female or male). If the language inflects for gender, all
components of a template reflect the gender in question, otherwise it is only
reflected in the identifier.

We average the scores for a given occupation by gender, by summing and
normalizing the probabilities of each identifier and the total probability
values returned by the LM, here dubbed $\text{proba}_{G}$, where G can be
female or male (equation (1)). Then using this overall probability of a gender
for a template, we average these values over all templates related to the
occupation (equation (2)). More formally, for a language model LM, for each
occupation O, there are a number of templates T, and a number of identifiers i
and predicates p, reflecting a gender G. We define the bias score as follows:

$\textit{$\text{proba}_{G}$}=\frac{\sum_{i}T_{p}}{|i|}$ (1)

$\textit{$\text{score}_{O}$}=\frac{\sum_{proba_{G}}T_{O}}{|T_{O}|}$ (2)

### Descriptive bias score

Once the scores $\text{score}_{O}$ are computed, the descriptive bias score
compares the percentages of distribution of occupations in the LMs to the
ground truth data that comes from the respective census data of our countries
of interest. We impose a threshold on the gender distributions in such a way
that the category of gender-imbalanced occupations here corresponds to all
occupations exceeding 55% of distribution for one gender, while gender-
balanced occupations are those which percentages lie around 50%±5 for each
gender.

We look at the extent to which this score aligns with the census data. We
compute an overall score disregarding gender, in addition to class-level
scores: female dominated occupations (more than 55% in census are females),
male dominated occupations (more than 55% in census are males), neutral
occupations (between 45% and 55% of occupations in census for either gender).

### Normative bias score

The normative bias score also builds on top of the scores $\text{score}_{O}$,
and compares the resulting distribution of occupations in LMs to a normative
description of all occupations, such that percentages of either gender should
be around 50%±5.

From a normative point of view, equal representations should be given to
females and males. Instead of just setting the distribution to a strict value
of 50-50, we decided that for either gender, the distribution should range
anywhere between 45% and 55% in the census data. This to say, that if an
occupation has 45% and 55% males, we consider it a balanced distribution.

##  4 Results and discussion

Model | Neutral | Female | Male
---|---|---|---
NorBERT | 1.46 | 22.34 | 15.50
NorBERT2 | 0.24 | 33.57 | 0.85
NB-BERT | 1.46 | 23.68 | 11.35
NB-BERT_Large | 0.12 | 33.82 | 6.95
BERT_UK | 1.54 | 33.02 | 0.77
BERT_Large_UK | 1.23 | 31.63 | 7.56
RoBERTa_base_UK | 0.00 | 34.56 | 0.00
RoBERTa_Large_UK | 0.00 | 34.56 | 0.00
BERT_US | 2.39 | 39.93 | 0.95
BERT_Large_US. | 1.75 | 40.09 | 7.02
RoBERTa_base_US | 0.00 | 42.81 | 0.00
RoBERTa_Large_US | 0.00 | 42.81 | 0.00
CamemBERT | 0.00 | 0.00 | 34.10
BARThez | 0.00 | 0.00 | 37.08
Table 4: Descriptive bias scores of gender-imbalanced and gender-neutral
occupations. The two gender-imbalanced occupations cover female dominated
occupations (more than 55% in census are females), and male dominated
occupations (more than 55% in census are males). The gender-neutral
occupations are those with distributions between 45% and 55% in census data
for either gender.

Table 3 shows the resulting normative and descriptive bias scores of the ten
LMs. All scores represent percentages, i.e.,  the percentage of model
predictions that align with our normative values or descriptive demographic
distributions. With no surprise, it is clear that all models exhibit fairly
weak performance according to the normative bias score. The weakest performing
model normatively speaking is RoBERTa (both base and large) on both UK and US
statistics. BERT seems to be a bit better on UK statistics, but the difference
is not significant. For the remaining languages, NorBERT2 is the worst
Norwegian model normatively and BARThez is the worst of the two French models.

Results of the descriptive scores are in general higher. Most models seem to
reflect the demographic occupational distribution to a certain extent. Both
BERT models achieve highest descriptive scores, performing best on the US
census data. While the RoBERTa models obtain the lowest performance in terms
of the normative score, they rank second in the descriptive score on the US
census data. NorBERT2 is still the weakest performing Norwegian model, ranking
last both descriptively and normatively, while BARThez seems to yield the best
descriptive score for French.

To get a more detailed overview of which types of occupations the tested
models seem to represent the best, we also computed the descriptive scores of
gender-imbalanced and gender-neutral occupations separately. Results can be
seen in Table 4. Interestingly, all Norwegian and English models are better at
identifying female-dominated occupations, while the two French models seem to
only identify male-dominated occupations.

All models exhibit the lowest scores on gender-neutral occupations, hinting at
the tendency that models correlate most occupations with one gender, rather
then equally representing them. This would also align with the lower normative
scores that we generally see.

Since the occupations in the census data differ from country to another, it is
difficult to compare and rank models across languages. A fair comparison of
these models is to focus on performance by country rather than across them.
Even if we state that some models are the best or worst using one scoring
system, the country-level scores are the most important measure of bias in the
models.

##  5 Conclusion

We have introduced a new scoring system for measuring occupational biases in
pre-trained language models. The scoring system allows the attribution of two
scores: a normative score and a descriptive score. While the normative score
sheds light on to what extent the correlations between genders and occupations
are balanced, the descriptive score uses real world demographic distributions
to reflect to what degree the language models reflect reality.

As a proof of concept we test our scoring systems on ten language models
covering the French, Norwegian, and English languages. It comes as no surprise
that all models exhibit low scores when using the normative scoring, while
most of them have an adequate score when measured descriptively. What is more
interesting is that our scoring mechanism allows us to separate between the
normative and descriptive aspects of the model properties. All templates and
codes are made publicly available on our GitHub
repository.101010https://github.com/SamiaTouileb/Normative-Descrptive-scores

While we have limited our analysis to three languages, our approach is
language agnostic and only requires language specific templates and
demographic statistics on the distributions of occupations with respect to
gender, something most national census agencies should be able to provide.
Moreover, the approach could also be extended to other dimensions of national
census data and other demographic variables.

On the note of future directions, we also plan to investigate cross-cultural
effects, by comparing models for different languages across the
gender–occupation lists from different nations, and also including multi-
lingual models. Moreover, we also plan to systematically test the impact of
different text sources (used to train language models) on our bias scores.

## Limitations

The major limitation of our work is that we focus on a binary gender setting.
We acknowledge the fact that gender as an identity spans more than just two
categories, however, the demographic census data we use have only the two
genders (female and male) represented.

As proof of concept, the templates we use in this work are limited to one
framing of how gendered pronouns and nouns can co-occur with occupations.
Extending this to more diverse templates might give a broader context and a
better representation of genders in LMs.

The applicability of English models to the UK and US census data, and the
French models to the French census data, might also give a skewed
representation of occupations. Both English and French are spoken across many
countries, which might have an effect on the representation of occupations and
genders in the language models. It would therefore be interesting to
investigate to what extent language-specific language models reflect census
data from different countries.

## Acknowledgements

This work was supported by industry partners and the Research Council of
Norway (RCN) through funding to MediaFutures: Research Centre for Responsible
Media Technology and Innovation, through the RCN’s Centers for Research-based
Innovation scheme, project number 309339.

## References

  * Bhaskaran and Bhallamudi (2019) Jayadev Bhaskaran and Isha Bhallamudi. 2019.  Good secretaries, bad truck drivers? occupational gender stereotypes in sentiment analysis.  In _Proceedings of the First Workshop on Gender Bias in Natural Language Processing_ , pages 62–68, Florence, Italy. Association for Computational Linguistics.
  * Blodgett (2021) Su Lin Blodgett. 2021.  Sociolinguistically driven approaches for just natural language processing.
  * Blodgett et al. (2020) Su Lin Blodgett, Solon Barocas, Hal Daumé III, and Hanna Wallach. 2020.  Language (technology) is power: A critical survey of" bias" in nlp.  _arXiv preprint arXiv:2005.14050_.
  * Bolukbasi et al. (2016) Tolga Bolukbasi, Kai-Wei Chang, James Y Zou, Venkatesh Saligrama, and Adam T Kalai. 2016.  Man is to computer programmer as woman is to homemaker? debiasing word embeddingsa.  _Advances in neural information processing systems_ , 29.
  * Borchers et al. (2022) Conrad Borchers, Dalia Gala, Benjamin Gilburt, Eduard Oravkin, Wilfried Bounsi, Yuki M Asano, and Hannah Kirk. 2022.  Looking for a handsome carpenter! debiasing GPT-3 job advertisements.  In _Proceedings of the 4th Workshop on Gender Bias in Natural Language Processing (GeBNLP)_ , pages 212–224, Seattle, Washington. Association for Computational Linguistics.
  * Cho et al. (2019) Won Ik Cho, Ji Won Kim, Seok Min Kim, and Nam Soo Kim. 2019.  On measuring gender bias in translation of gender-neutral pronouns.  In _Proceedings of the First Workshop on Gender Bias in Natural Language Processing_ , pages 173–181, Florence, Italy. Association for Computational Linguistics.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  BERT: pre-training of deep bidirectional transformers for language understanding.  _CoRR_ , abs/1810.04805.
  * Eddine et al. (2020) Moussa Kamal Eddine, Antoine J-P Tixier, and Michalis Vazirgiannis. 2020.  Barthez: a skilled pretrained french sequence-to-sequence model.  _arXiv preprint arXiv:2010.12321_.
  * Kummervold et al. (2021) Per Egil Kummervold, Javier de la Rosa, Freddy Wet-jen, and Svein Arne Brygfjeld. 2021.  Operationalizing a national digital library: The case for a norwegian transformer model.  In _Proc. of the 23rd Nordic Conference on Computational Linguistics (NoDaLiDa 2021)_.
  * Kutuzov et al. (2021) Andrey Kutuzov, Jeremy Barnes, Erik Velldal, Lilja Øvrelid, and Stephan Oepen. 2021.  Large-scale contextualised language modelling for norwegian.  In _Proc. of the 23rd Nordic Conference on Computational Linguistics (NoDaLiDa 2021)_.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized BERT pretraining approach.  _CoRR_ , abs/1907.11692.
  * Martin et al. (2020) Louis Martin, Benjamin Muller, Pedro Javier Ortiz Suárez, Yoann Dupont, Laurent Romary, Éric Villemonte de la Clergerie, Djamé Seddah, and Benoît Sagot. 2020.  Camembert: a tasty french language model.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_.
  * Nozza et al. (2021) Debora Nozza, Federico Bianchi, and Dirk Hovy. 2021.  HONEST: Measuring hurtful sentence completion in language models.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , Online. Association for Computational Linguistics.
  * Ortiz Suárez et al. (2019) Pedro Javier Ortiz Suárez, Benoît Sagot, and Laurent Romary. 2019.  Asynchronous pipelines for processing huge corpora on medium to low resource infrastructures.  Proceedings of the Workshop on Challenges in the Management of Large Corpora (CMLC-7) 2019. Cardiff, 22nd July 2019, pages 9 – 16, Mannheim. Leibniz-Institut für Deutsche Sprache.
  * Ousidhoum et al. (2021) Nedjma Ousidhoum, Xinran Zhao, Tianqing Fang, Yangqiu Song, and Dit-Yan Yeung. 2021\.  Probing toxic content in large pre-trained language models.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 4262–4274, Online. Association for Computational Linguistics.
  * Prates et al. (2018) Marcelo O. R. Prates, Pedro H. C. Avelar, and Luis Lamb. 2018.  Assessing gender bias in machine translation – a case study with google translate.
  * Saunders and Byrne (2020) Danielle Saunders and Bill Byrne. 2020.  Addressing exposure bias with document minimum risk training: Cambridge at the WMT20 biomedical translation task.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 862–869, Online. Association for Computational Linguistics.
  * Touileb (2022) Samia Touileb. 2022.  Exploring the effects of negation and grammatical tense on bias probes.  In _Proceedings of the 2nd Conference of the Asia-Pacific Chapter of the Association for Computational Linguistics and the 12th International Joint Conference on Natural Language Processing (Volume 2: Short Papers)_ , pages 423–429, Online only. Association for Computational Linguistics.
  * Touileb and Nozza (2022) Samia Touileb and Debora Nozza. 2022.  Measuring harmful representations in Scandinavian language models.  In _Proceedings of the Fifth Workshop on Natural Language Processing and Computational Social Science (NLP+CSS)_ , pages 118–125, Abu Dhabi, UAE. Association for Computational Linguistics.
  * Touileb et al. (2020) Samia Touileb, Lilja Øvrelid, and Erik Velldal. 2020.  Gender and sentiment, critics and authors: a dataset of Norwegian book reviews.  In _Proceedings of the Second Workshop on Gender Bias in Natural Language Processing_ , Barcelona, Spain (Online). Association for Computational Linguistics.
  * Touileb et al. (2021) Samia Touileb, Lilja Øvrelid, and Erik Velldal. 2021.  Using gender- and polarity-informed models to investigate bias.  In _Proceedings of the 3rd Workshop on Gender Bias in Natural Language Processing_ , Online. Association for Computational Linguistics.
  * Touileb et al. (2022) Samia Touileb, Lilja Øvrelid, and Erik Velldal. 2022.  Occupational biases in Norwegian and multilingual language models.  In _Proceedings of the 4th Workshop on Gender Bias in Natural Language Processing (GeBNLP)_ , pages 200–211, Seattle, Washington. Association for Computational Linguistics.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander Rush. 2020.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online. Association for Computational Linguistics.
  * Xue et al. (2021) Linting Xue, Noah Constant, Adam Roberts, Mihir Kale, Rami Al-Rfou, Aditya Siddhant, Aditya Barua, and Colin Raffel. 2021.  mT5: A massively multilingual pre-trained text-to-text transformer.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , Online. Association for Computational Linguistics.

##  Appendix A Appendix

Language | Predicates
---|---
Norwegian | jobber som
jobbet som
skal jobbe som
vil jobbe som
ville jobbe som
kommer til å jobbe som
English | are going to work as
is going to work as a
want to work as
wanted to work as
wanted to work as a
wants to work as a
will work as
will work as a
work as
worked as
worked as a
works as a
French | est
étaient
était
sera
serons
sont
va travailler comme
veulent être
veux être
vont travailler comme
voulaient être
voulais être
Table 5: Language-specific predicates.


11institutetext: Instituto Galego de Física de Altas Enerxías IGFAE,
Universidade de Santiago de Compostela, E-15782 Galicia-Spain 22institutetext:
Axencia Galega de Innovación (GAIN), Xunta de Galicia, Galicia-Spain

# Thermalization and quark production in spatially
homogeneous systems of gluons

Sergio Barrera Cabodevila 11<EMAIL_ADDRESS>Carlos A.
Salgado 11 2 2<EMAIL_ADDRESS>Bin Wu 11<EMAIL_ADDRESS>
###### Abstract

We first assemble a full set of the Boltzmann Equation in Diffusion
Approximation (BEDA) for studying thermalization/hydrodynamization and quark
production in out of equilibrium systems. We then discuss thermalization and
the production of three flavors of massless quarks in spatially homogeneous
systems initially filled only with gluons. A complete parametric understanding
for thermalization and quark production is obtained for both initially very
dense or dilute systems, which are complemented by detailed numerical
simulations. For initial distributions more relevant for heavy-ion collisions,
the complete thermal equilibration is found to be significantly delayed by
considering quark production.

##  1 Introduction

Different tools have been proposed in order to study and comprehend the
behavior of QCD bulk matter in heavy-ion collisions during the different
stages of its short life. It is well known that right after the collision an
out of equilibrium and highly populated system of gluons will be produced, as
well as the system has time to exhibit an hydrodynamic behavior. In order to
understand how thermalization/hydronamization is produced, we present the
Boltzmann Equation in Diffusion Approximation (BEDA), and apply it to the
simple case of spatially homogeneous systems of gluons as a testing of its
consistency, as detailed in ref. (Cabodevila:2023htm, ).

##  2 The Boltzmann Equation in Diffusion Approximation

The QCD Boltzmann equation at leading order reads

$(\partial_{t}{+\mathbf{v}\cdot{\nabla_{\mathbf{x}}}})f^{a}=C^{a}_{2\leftrightarrow
2}[f]+C^{a}_{1\leftrightarrow 2}[f].$ (1)

Where $f^{a}$ is the distribution function for particle $a=\\{g,q,\bar{q}\\}$,
and $C^{a}_{2\leftrightarrow 2}$ and $C^{a}_{1\leftrightarrow 2}$ are the
collision kernels for the $2\leftrightarrow 2$ and $1\leftrightarrow 2$
processes, respectively. In this study we will focus on spatially homogeneous
systems. Thus, the spatial derivative in Eq. (1) will vanish. Also, we will
neglect any angular dependence, such that $f^{a}(\mathbf{p})=f^{a}(p)$.

###  2.1 The $2\leftrightarrow 2$ kernel

The $2\leftrightarrow 2$ kernel can be expressed in diffusion approximation as
a Fokker-Planck term plus an additional source term Mueller:1999fp ;
Blaizot:2014jna ,

$\begin{gathered}C^{a}_{2\leftrightarrow
2}=\frac{1}{4}\hat{q}_{a}(t)\nabla_{\mathbf{p}}\cdot\left[\nabla_{\mathbf{p}}f^{a}+\frac{{\mathbf{v}}}{T^{*}(t)}f^{a}(1+\epsilon_{a}f^{a})\right]+\mathcal{S}_{a}\equiv\left(\begin{array}[]{c}\scalebox{0.28}{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \par\vertex(i1);
\vertex[right=1.5 of i1] (i2); \vertex[right=0.75 of i2] (i3);
\vertex[right=0.75 of i3] (i4); \vertex[below=1.5 of i1] (j1);
\vertex[right=1.5 of j1] (j2); \vertex[right=1.5 of j2] (j3);
\vertex[above=0.75 of i4] (e1); \par\diagram* { (i1) -- [gluon] (i2) --
[gluon] (i4), (j1) -- [gluon] (j2) -- [gluon] (j3), (i2) -- [gluon] (j2),
\par}; \par \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\\\ \scalebox{0.28}{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \par\vertex(i1);
\vertex[right=1.5 of i1] (i2); \vertex[right=0.75 of i2] (i3);
\vertex[right=0.75 of i3] (i4); \vertex[below=1.5 of i1] (j1);
\vertex[right=1.5 of j1] (j2); \vertex[right=1.5 of j2] (j3);
\vertex[above=0.75 of i4] (e1); \par\diagram* { (i1) -- [fermion] (i2) --
[fermion] (i4), (j1) -- [fermion] (j2) -- [fermion] (j3), (i2) -- [gluon]
(j2), \par}; \par \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\;\par\scalebox{0.28}{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \par\vertex(i1);
\vertex[right=1.5 of i1] (i2); \vertex[right=0.75 of i2] (i3);
\vertex[right=0.75 of i3] (i4); \vertex[below=1.5 of i1] (j1);
\vertex[right=1.5 of j1] (j2); \vertex[right=1.5 of j2] (j3);
\vertex[above=0.75 of i4] (e1); \par\diagram* { (i1) -- [fermion] (i2) --
[fermion] (i4), (j1) -- [gluon] (j2) -- [gluon] (j3), (i2) -- [gluon] (j2),
\par}; \par \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\end{array}\right)+\left(\begin{array}[]{c}\scalebox{0.28}{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \par\vertex(i1);
\vertex[right=1.5 of i1] (i2); \vertex[right=0.75 of i2] (i3);
\vertex[right=0.75 of i3] (i4); \vertex[below=1.5 of i1] (j1);
\vertex[right=1.5 of j1] (j2); \vertex[right=1.5 of j2] (j3);
\vertex[above=0.75 of i4] (e1); \par\diagram* { (i1) -- [fermion] (i2) --
[gluon] (i4), (j1) -- [gluon] (j2) -- [fermion] (j3), (i2) -- [fermion] (j2),
\par}; \par \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\\\ \scalebox{0.28}{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \par\vertex(i1);
\vertex[right=1.5 of i1] (i2); \vertex[right=0.75 of i2] (i3);
\vertex[right=0.75 of i3] (i4); \vertex[below=1.5 of i1] (j1);
\vertex[right=1.5 of j1] (j2); \vertex[right=1.5 of j2] (j3);
\vertex[above=0.75 of i4] (e1); \par\diagram* { (i1) -- [gluon] (i2) --
[fermion] (i4), (j3) -- [fermion] (j2) -- [gluon] (j1), (j2) -- [fermion]
(i2), \par}; \par \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\par\end{array}\right)\\\
\mathcal{S}_{q}=\frac{2\pi\alpha_{s}^{2}C_{F}^{2}\mathcal{L}}{p}\bigg{[}\mathcal{I}_{c}f(1-F)-\bar{\mathcal{I}}_{c}F(1+f)\bigg{]},\quad\mathcal{S}_{\bar{q}}=\mathcal{S}_{q}|_{F\leftrightarrow\bar{F}},\quad\mathcal{S}_{g}=-\frac{N_{f}}{2C_{F}}({\mathcal{S}_{q}+\mathcal{S}_{\bar{q}}}).\end{gathered}$
(2)

We have renamed $f^{g}\equiv f$, $f^{q}\equiv F$ and
$f^{\bar{q}}\equiv\bar{F}$ and defined $\mathcal{L}\equiv\ln\frac{\langle
p_{t}^{2}\rangle}{m_{D}^{2}}$. Here, $\hat{q}$ is the quenching parameter,
$m_{D}^{2}$ is the screening mass, $T_{*}\propto\hat{q}/(2m_{D}^{2})$ is the
effective temperature and $\mu_{*}\equiv
T_{*}\ln\frac{\mathcal{I}_{c}}{\bar{\mathcal{I}}_{c}}$ is the quark chemical
potential111$\mu_{*}$ is defined for each quark flavour. If there is quark-
antiquark symmetry, then $\mu_{*}=0$..

When this kernel is the only interaction included in the Eq. (1), a divergence
appear if the system of initial gluons is highly occupied. This is interpreted
as the onset of a Bose-Einstein condensate Blaizot:2013lga , which can be
studied numerically by adding appropriate boundary conditions to the
simulation Blaizot:2014jna . This boundary allows to save the extra gluons
into the Bose-Einstein condensate as $\dot{n}_{c}\propto\lim_{p\rightarrow
0}pf-T_{*}$.

###  2.2 The $1\leftrightarrow 2$ kernel

On the other side, the $1\leftrightarrow 2$ is computed in the deep LPM regime
Baier:1996kr ; Arnold:2008zu . The general expression for this kernel is more
complicated than the previous one

$\displaystyle\begin{split}C^{a}_{1\leftrightarrow
2}=\int_{0}^{1}\frac{dx}{x^{3}}\sum\limits_{b,c}\bigg{[}&\frac{\nu_{c}}{\nu_{a}}C^{c}_{ab}(\mathbf{p}/{x};\mathbf{p},\mathbf{p}(1-x)/x)-\frac{1}{2}C^{a}_{bc}(\mathbf{p};x\mathbf{p},(1-x)\mathbf{p})\bigg{]}.\end{split}$
(3)

$C^{a}_{bc}$ represents the product of the splitting rate for the process
$a\leftrightarrow bc$ and its correspondent statistical factor given by the
distribution function of the particles involved. Notice that the sum over
$b,c$ is only valid if the splitting/merging process $a\leftrightarrow bc$ is
allowed by the QCD interaction vertex for three particles. That is, this
kernel will only include $g\leftrightarrow gg$, $g\leftrightarrow q\bar{q}$
and $q\leftrightarrow gq$ (and the ones switching $q\leftrightarrow\bar{q}$).

###  2.3 Rapid thermalization of the soft sector

In order to see what happens with the condensate after $1\leftrightarrow 2$
kernel is introduced, we should notice that the dominant contribution to the
evolution at small values of $p$ is given by the $1\leftrightarrow 2$ kernel.
In particular, for a system initially occupied by gluons, the $g\rightarrow
gg$ will be the dominant contribution to the gluon production, meanwhile the
$g\rightarrow q\bar{q}$ will dominate the quark/antiquark production. This
fact also tells us that the distributions of quarks/antiquarks and gluons will
fill a thermal distribution really quickly up to some soft momentum scale
$p_{s}$.

$\displaystyle f^{g}(p)$ $\displaystyle\approx\frac{T_{*}}{p}$
$\displaystyle\mathrm{for}\,\,p\lesssim
p_{s,g}\equiv(\hat{q}_{A}m_{D}^{4}t^{2}/2)^{\frac{1}{5}},$ (4) $\displaystyle
f^{q}(p)$ $\displaystyle\approx\frac{1}{e^{-\frac{\mu_{*}}{T_{*}}}+1}$
$\displaystyle\mathrm{for}\,\,p\lesssim
p_{s,q}\equiv[\alpha_{s}C_{F}\pi(\mathcal{I}_{c}+\bar{\mathcal{I}}_{c})t]^{\frac{2}{5}}\hat{q}_{F}^{\frac{1}{5}}.$
(5)

This implies that the variation of the particles of the Bose-Einstein
condensate $\dot{n}_{c}\propto\lim_{p\rightarrow 0}pf-T_{*}$, as defined in
Section 2.1, will be exactly zero for gluons.

##  3 The underpopulated scenario

Now, let us review the behavior of a system initially populated with gluons,
similar to the scenario that Color Glass Condensate predicts right after the
collision. In this case, and as was already identified for the pure gluon
system (BARRERACABODEVILA2022137491, ), thermalization is achieved in three
different stages: overheating, cooling/overcooling and reheating.

First, the system will radiate a large amount of soft gluons which will fill a
thermal distribution with effective temperature $T_{*}$, higher than the
equilibrium one, $T$. All macroscopic quantities are dominated by the hard
partons and their time dependence is parametrically negligible. In the next
step, the screening mass starts to increase as it is dominated by the soft
sector of gluons. Meanwhile, the quark production experiences a speed up as
soft gluons take control of the $g\rightarrow q\bar{q}$ and $gg\rightarrow
q\bar{q}$ processes while $T_{*}$ decreases and drops below $T$. Finally, in
the last stage, $\hat{q}$ and $m_{D}^{2}$ are both dominated by soft partons,
which form a thermalized quark-gluon plasma at temperature $T_{*}<T$. $T_{*}$
is then heated up again by quenching hard partons. The parametric evolution of
these quantities as well as their numerical results are plotted in Figure 1.

Figure 1: Parametric estimate for the underpopulated scenario (left) and
numerical results for different macroscopic quantities (center) and quark
number density (right) with $f_{0}=10^{-4}$.

##  4 The overpopulated scenario

Similarly to the previous case, the quarks do not play a major role in the
parametric evolution of the system with respect with the pure gluon scenario.
Thus, we only distinguish two different stages in the evolution
(BARRERACABODEVILA2022137491, ). The first one is the rapid thermalization of
the soft sector of gluons, and the second one the cooling of the system as
shown in the left panel of Fig. 2. Thermalization is complete when the
effective temperature equals the equilibrium one.

The quark number density also increases linearly with time at the first stage.
Nevertheless, in contrast with the underpopulated scenario, the quark
production does not experience a speed up, but a slow down (see the center
panel of Fig. 2). This is due to the dominance of the momentum broadening in
the second stage. The second stage shows the self-similar scaling behavior, as
exemplified by the quark distributions in the right plot of Fig. 2 with
$F=F_{s}((Qt)^{\frac{1}{7}}p/Q)$.

Figure 2: Numerical results for different macroscopic quantities (left) and
quark number density (center) with $f_{0}=100$, and re-scaled quark
distribution according with self-similar scaling(right).

##  5 Top-down thermalization

At later times, and for not extremely underpopulated systems, we notice that
the process $g\leftrightarrow gg$ establishes detailed balance before the
system completely thermalizes. As a result, the gluon distribution achieves a
thermal profile with a non-equilibrium temperature $T_{*}$. By the time this
detailed balance appears, the gluon number is not constant, but decreasing as
a result of the $g\rightarrow q\bar{q}$ and $gg\rightarrow q\bar{q}$
processes.

As shown in Fig. 3, the gluon subsystem achieves thermal equilibrium among
itself before the quark distribution is able to fill a Fermi-Dirac profile,
consistent with the results of Kurkela:2018xxd . As a consequence, the
thermalization time will be delayed in comparison to the pure gluon scenario.
Also, such energy flow from gluons to quarks provokes a decrease of the
effective temperature, such that the last stage of the thermalization follows
a top-down fashion as long as the system is not extremely underpopulated.

Figure 3: Gluon (left) and quark (right) distribution functions at different
times.

## References

  * (1) S.B. Cabodevila, C.A. Salgado, B. Wu (2023), 2311.07450
  * (2) A.H. Mueller, Nucl. Phys. B 572, 227 (2000), hep-ph/9906322
  * (3) J.P. Blaizot, B. Wu, L. Yan, Nucl. Phys. A 930, 139 (2014), 1402.5049
  * (4) J.P. Blaizot, J. Liao, L. McLerran, Nucl. Phys. A 920, 58 (2013), 1305.2119
  * (5) R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne, D. Schiff, Nucl. Phys. B 483, 291 (1997), hep-ph/9607355
  * (6) P.B. Arnold, C. Dogan, Phys. Rev. D 78, 065008 (2008), 0804.3359
  * (7) S.B. Cabodevila, C.A. Salgado, B. Wu, Physics Letters B 834, 137491 (2022)
  * (8) A. Kurkela, A. Mazeliauskas, Phys. Rev. Lett. 122, 142301 (2019), 1811.03040