#  How to use GP: Effects of the mean function and hyperparameter selection on
Gaussian Process regression

Seung-gyu Hwang*  | Benjamin L’Huillier*  | Ryan E. Keeley*  | M. James Jee  Arman Shafieloo*  |

###### Abstract

Gaussian processes have been widely used in cosmology to reconstruct
cosmological quantities in a model-independent way. However, the validity of
the adopted mean function and hyperparameters, and the dependence of the
results on the choice have not been well explored. In this paper, we study the
effects of the underlying mean function and the hyperparameter selection on
the reconstruction of the distance moduli from type Ia supernovae. We show
that the choice of an arbitrary mean function affects the reconstruction: a
zero mean function leads to unphysical distance moduli and the best-fit
$\Lambda\text{CDM}$ to biased reconstructions. We propose to marginalize over
a family of mean functions and over the hyperparameters to effectively remove
their impact on the reconstructions. We further explore the validity and
consistency of the results considering different kernel functions and show
that our method is unbiased.

##  1 Introduction

The nature of the accelerating universe is one of the biggest questions in
modern physics and cosmology. In a matter-dominated universe, the expansion
speed of the Universe should gradually decrease, but observations suggest that
the cosmic expansion is accelerated by some unknown component of the Universe,
dark energy (DE), e.g. [1, 2].

The role of DE is to create a repulsive force against gravity. DE can be
described as a perfect fluid with equation of state (EoS) $w=p/\rho$, which is
the ratio between pressure $p$ and density $\rho$ of the fluid. The simplest
model is the presence of a positive constant in the Einstein field equations,
which can be modeled as a fluid with $w=-1$.

However, as observations become more and more precise, some tensions arise
such as the Hubble tension between local measurements of the Hubble-Lemaître
constant $H_{0}$ [3, 4] and inferences of this parameter from high-redshift
experiments [5]. These tensions, in addition to our ignorance of the nature of
DE, may lead to modifications to the standard model. A possibility is that DE
is dynamical, i.e., that the EoS $w$ is an evolving function of the redshift.
For a flat Friedmann–Lemaître–Robertson–Walker metric with EoS $w(z)$, the
expansion history $h(z)=H(z)/H_{0}$ is described by

$h^{2}(z)=\Omega_{\mathrm{m0}}(1+z)^{3}+(1-\Omega_{\mathrm{m0}})\exp{\left(3\int_{0}^{z}\frac{1+w(z^{\prime})}{1+z^{\prime}}dz^{\prime}\right)},$
(1.1)

where $\Omega_{\text{m0}}$ is the matter energy density parameter at the
present time.

Two main approaches are carried out in the literature to test cosmology beyond
the concordance model. On one hand, theoreticians are building new models,
usually adding new degrees of freedom, then compare them to the data and
obtain constraints on their parameters. While this approach can give tight
constraints on the parameters, the results may be biased towards the
considered models, and may even be meaningless if the model is not a good
description of the reality. On the other hand, model-independent, data-driven
approaches can be applied to infer cosmological quantities without assuming a
model, and the resulting reconstructions are thus not biased by a particular
model assumption. These methods can be more flexible, and thus allow the
detection of features in the data that could be missed by traditional model-
fitting approaches. However, they can be prone to overfitting and the
interpretation of the results is not always straightforward.

Among these model-independent methods, Gaussian process (GP) regression [6]
has become particularly popular in the last decade or so. Refs. [7, 8, 9]
introduced GP to analyze cosmological observations phenomenologically. In
addition, the flexibility of model-independent approaches enable the discovery
of interesting trends. In a nutshell, GP modeling implicitely makes the
assumption that the data and the theoretical function to be reconstructed are
jointly Gaussian, and that their distribution can be fully described by a mean
function and a covariance kernel, itself defined by its hyperparameters. While
a large amount of work has been devoted to the choice of the kernel, few
studies have investigated the dependence on the choice of the mean function
and hyperparameters (e.g., [7, 10]). Most notably, most of the literature
simply uses zero as the mean function, and fixes the hyperparameters to the
value that maximize the marginal likelihood, claiming that zero mean function
is a general case. As we will show, the choice of this mean function is not
necessarily well motivated, in particular for non-stationary functions such as
cosmic distances or expansions. We will also demonstrate that the choice of
the hyperparameters is also dependent on that of the mean function. Moreover,
we will suggest a rigorous approach to attain a completely Bayesian analysis
with the GPR.

In this study, we focus on several aspects of GP regression, in particular,
the effect of the mean function and the choice of the hyperparameters. In
Section 2, we describe the method and the mock data used in this study.
Section 3 presents our main results on the factors that affect reconstruction
results. We summarize our findings in Section 4.

##  2 Method and Data

###  2.1 Gaussian process regression

A GP is a collection of random variables, such that every finite collection of
those random variables has a multivariate normal distribution. The
distribution of $f$ at any finite number of points
$\boldsymbol{x}=(x_{1},\dots,x_{N})$ is specified by a mean function
$\mu(\mathrm{\boldsymbol{x}})$ and a covariance matrix $\mathbf{K}$,

$\displaystyle f(\mathrm{\boldsymbol{x}})$
$\displaystyle\sim\mathcal{GP}(\mu(\mathrm{\boldsymbol{x}}),\mathbf{K}),$
(2.1) where the covariance matrix $\mathbf{K}$ is defined by
$\displaystyle[\mathbf{K}]_{ij}$ $\displaystyle=k(x_{i},x_{j}),$ (2.2)

where $k(x,x^{\prime})$ is the covariance kernel.

Specifically, the mean and covariance functions are

$\displaystyle\mu(\boldsymbol{x})=\mathbb{E}[f(\boldsymbol{x})],$ (2.3)
$\displaystyle
k(x_{i},x_{j})=\mathbb{E}[(f(x_{i})-\mu(x_{i}))(f(x_{j})-\mu(x_{j}))].$

A covariance function is a kernel function describing the correlation between
the random variables. While many shapes of kernels have been considered in the
literature, in this study, we will focus on the square exponential (SE)
kernel, defined as

$k(x,x^{\prime})=\sigma_{f}^{2}\exp{\left[-\frac{(x-x^{\prime})^{2}}{2l^{2}}\right]}$
(2.4)

The squared exponential kernel has two hyperparameters that control the
correlation length ($l$) and magnitude ($\sigma_{f}$) of the fluctuations.
Section 3.2.1 explores the role of the kernel.

Gaussian process regression (GPR) is used in cosmology to reconstruct a
desired quantity from a discrete number of observations [7, 8, 10, 11, 12,
13]. Two approaches can be used to perform a regression.

####  2.1.1 Prior approach: forward-modeling

The first approach is to model the function to be reconstructed as a GP. We
will refer to this approach as the _prior approach_ , or _untrained GP_ [7, 8,
9, 14, 12, 13]. In essence, it consists in the following steps:

  1. 1.

Choose a mean function $\boldsymbol{\mu^{*}}$, a covariance kernel $k$, and
its hyperparameters and calculate the covariance matrix
$\mathbf{K}_{**}=K(\boldsymbol{X^{*}},\boldsymbol{X^{*}})$.

  2. 2.

Generate a collection of random functions $\boldsymbol{f^{*}}$ following the
Gaussian distribution
$\boldsymbol{f^{*}}\sim\mathcal{N}(\boldsymbol{\mu^{*}},\mathbf{K}_{**})$.

  3. 3.

For each of these random functions $f^{*}$, evaluate the goodness-of-fit, for
instance by calculating the $\chi^{2}$ for the samples $\boldsymbol{f^{*}}$
evaluated against the observation $\boldsymbol{y}$.

  4. 4.

Obtain the conditional distribution $\boldsymbol{f^{*}}|\boldsymbol{y}$
weighted by the goodness-of-fit.

In other words, we generate the samples $f^{*}$ following a GP, and
individually calculate the goodness-of-fit for the observation. The data
covariance is not used to generate $f^{*}$, but is only used to evaluate the
$\chi^{2}$. In practice, this approach can be computationally expensive, since
we have to generate a large quantity of samples. However, it allows for
forward-modeling and combining data sets [7, 8, 9, 14, 12, 13].

####  2.1.2 Posterior approach: GP regression

Because of the computational cost of the prior approach, a more commonly used
approach is to model the data and the function to be reconstructed as a GP,
i.e., assume that they are jointly distributed and follow a Gaussian
distribution with a chosen mean and covariance. This approach will be referred
to as the _posterior approach_ , or _trained GP_ [10, 11]. In order to take
into account the observational error in the process, the covariance matrix of
the observation can be added to the kernel matrix,

$\begin{bmatrix}\boldsymbol{y}\\\
\boldsymbol{f^{*}}\end{bmatrix}\mathtt{\sim}~{}\mathcal{N}\left(\begin{bmatrix}\boldsymbol{\mu}\\\
\boldsymbol{\mu^{*}}\end{bmatrix},\begin{bmatrix}K(\boldsymbol{X},\boldsymbol{X})+C&K(\boldsymbol{X},\boldsymbol{X^{*}})\\\
K(\boldsymbol{X^{*}},\boldsymbol{X})&K(\boldsymbol{X^{*}},\boldsymbol{X^{*}})\end{bmatrix}\right),$
(2.5)

where $\boldsymbol{y}$ is the observation, $\boldsymbol{f^{*}}$ is the
reconstruction, $\boldsymbol{\mu^{(*)}}$ is the input mean function, $K$ is a
kernel matrix given by
$[K(\boldsymbol{X^{(*)}},\boldsymbol{X^{(*)}})]_{ij}=k(x^{(*)}_{i},x^{(*)}_{j})$,
$C$ is the covariance matrix of the observation, and * denotes the test points
where we wish to reconstruct the function. The reconstruction is done with the
conditional distribution to model a physical quantity given the observation as
a GP,

$\boldsymbol{f^{*}}|\boldsymbol{y}\sim\mathcal{N}\left(\overline{\boldsymbol{f^{*}}},\operatorname{\mathrm{cov}}(\boldsymbol{f^{*}})\right).$
(2.6)

For a given input mean function and hyperparameters, the GPR is readily
computed with the conditional distribution via the following matrix
calculation, thanks to the matrix inversion lemma:

$\displaystyle\overline{\boldsymbol{f^{*}}}=\boldsymbol{\mu^{*}}+K(\boldsymbol{X^{*}},\boldsymbol{X})[K(\boldsymbol{X},\boldsymbol{X})+C]^{-1}(\boldsymbol{y}-\boldsymbol{\mu}),$
(2.7)
$\displaystyle\operatorname{\mathrm{cov}}(\boldsymbol{f^{*}})=K(\boldsymbol{X^{*}},\boldsymbol{X^{*}})-K(\boldsymbol{X^{*}},\boldsymbol{X})[K(\boldsymbol{X},\boldsymbol{X})+C]^{-1}K(\boldsymbol{X},\boldsymbol{X^{*}}).$

We note that in the posterior approach, it is not necessary to actually draw
the individual samples $\boldsymbol{f^{*}}$, and if one is only interested in
the distribution of $\boldsymbol{f^{*}}$, all the information is contained in
$\boldsymbol{\overline{f^{*}}}$ and
$\operatorname{\mathrm{cov}}(\boldsymbol{f^{*}})$, which can significantly
speed up the computation.

In order to complete the regression, one needs to assess the goodness-of-fit
of the reconstruction. From Bayes’ theorem, the marginal likelihood of the
conditional distribution, eq. (2.6), is given by

$p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\mu},\mathcal{H})=\int{p(\boldsymbol{y}|\boldsymbol{f},\boldsymbol{X})p(\boldsymbol{f}|\boldsymbol{X},\boldsymbol{\mu},\mathcal{H})\mathrm{d}\boldsymbol{f}},$
(2.8)

where $\mathcal{H}$ is a set of hyperparameters of the kernel function. It is
important to note that the marginal likelihood is dependent on a choice of the
input mean function, kernel function, and its hyperparameters. In practice,
for computational reasons, the logarithm of the marginal likelihood (LML) is
often used.

####  2.1.3 Reconstruction of the derivative

GP can be used to reconstruct not only a function but also its derivative
since the derivative of a GP is also a GP. Assuming that the same kernel
function is used for the GP, the kernel function for the derivative is

$\mbox{cov}\left(f_{i},\frac{\partial{f_{j}}}{\partial{x_{j}}}\right)=\frac{\partial{k(x_{i},x_{j})}}{\partial{x_{j}}},\hskip
28.45274pt\mbox{cov}\left(\frac{\partial{f_{i}}}{\partial{x_{i}}},\frac{\partial{f_{j}}}{\partial{x_{j}}}\right)=\frac{\partial{{}^{2}k(x_{i},x_{j})}}{\partial{x_{i}}\partial{x_{j}}}.$
(2.9)

With the covariance for the derivative, the joint distribution of the
observation $\boldsymbol{y}$, the reconstruction $\boldsymbol{f^{*}}$, and its
derivative $\boldsymbol{f^{*}}^{\prime}$ are [10],

$\begin{bmatrix}\boldsymbol{y}\\\ \boldsymbol{f^{*}}\\\
\boldsymbol{f^{*}}^{\prime}\end{bmatrix}\mathtt{\sim}~{}\mathcal{N}\left(\begin{bmatrix}\boldsymbol{\mu}\\\
\boldsymbol{\mu^{*}}\\\
\boldsymbol{\mu^{*}}^{\prime}\end{bmatrix},\begin{bmatrix}K(\boldsymbol{X},\boldsymbol{X})+C&K(\boldsymbol{X},\boldsymbol{X^{*}})&K^{\prime}(\boldsymbol{X},\boldsymbol{X^{*}})\\\
K(\boldsymbol{X},\boldsymbol{X^{*}})^{T}&K(\boldsymbol{X^{*}},\boldsymbol{X^{*}})&K^{\prime}(\boldsymbol{X^{*}},\boldsymbol{X^{*}})\\\
K^{\prime}(\boldsymbol{X},\boldsymbol{X^{*}})^{T}&K^{\prime}(\boldsymbol{X^{*}},\boldsymbol{X^{*}})&K^{\prime\prime}(\boldsymbol{X^{*}},\boldsymbol{X^{*}})\end{bmatrix}\right),$
(2.10)

where

$[K^{\prime}(\boldsymbol{X},\boldsymbol{X^{*}})]_{ij}=\frac{\partial{k(x_{i},x_{j}^{*})}}{\partial{x_{j}^{*}}},\hskip
28.45274pt[K^{\prime\prime}(\boldsymbol{X^{*}},\boldsymbol{X^{*}})]_{ij}=\frac{\partial^{2}{k(x_{i}^{*},x_{j}^{*})}}{\partial{x_{i}^{*}}\partial{x_{j}^{*}}}$
(2.11)

Similarly to eq. (2.7), the conditional distribution for the derivative given
the observation also follow the Gaussian distribution:

$\displaystyle\begin{bmatrix}\boldsymbol{f}\\\
\boldsymbol{f}^{\prime}\end{bmatrix}|\boldsymbol{y}\sim\mathcal{N}\left(\begin{bmatrix}\bar{\boldsymbol{f}}\\\
\bar{\boldsymbol{f}^{\prime}}\end{bmatrix};\begin{bmatrix}\mathbf{A}-\mathbf{D}\mathbf{B}^{-1}\mathbf{D}^{T}\end{bmatrix}\right),$
(2.12)

where

$\displaystyle\begin{bmatrix}\bar{\boldsymbol{f}}\\\
\bar{\boldsymbol{f}^{\prime}}\end{bmatrix}$
$\displaystyle=\mathbf{D}\mathbf{B}^{-1}\begin{bmatrix}\boldsymbol{y}-\boldsymbol{\mu}\end{bmatrix},$
(2.13) $\displaystyle\mathbf{A}$
$\displaystyle=\begin{bmatrix}\mathbf{K}(\boldsymbol{X}^{*},\boldsymbol{X}^{*})&\mathbf{K}^{\prime}(\boldsymbol{X}^{*},\boldsymbol{X}^{*})\\\
\mathbf{K}^{\prime}(\boldsymbol{X}^{*},\boldsymbol{X}^{*})^{T}&\mathbf{K}^{\prime\prime}(\boldsymbol{X}^{*},\boldsymbol{X}^{*})\end{bmatrix}$
(2.14) $\displaystyle\mathbf{B}$
$\displaystyle=\mathbf{K}(\boldsymbol{X},\boldsymbol{X})+\mathbf{C}$ (2.15)
$\displaystyle\mathbf{D}^{T}$
$\displaystyle=\begin{bmatrix}\mathbf{K}(\boldsymbol{X},\boldsymbol{X}^{*})&\mathbf{K}^{\prime}(\boldsymbol{X},\boldsymbol{X}^{*})\end{bmatrix}.$
(2.16)

The shapes of $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{D}$ are respectively
$2n_{*}\times 2n_{*}$, $n\times n$, and $n_{*}\times n$, where $n_{*}$ is the
size of $\vec{f}*$ and $n$ the number of data points.

###  2.2 Data

We generate mock data based on the Pantheon compilation [15] for our fiducial
cosmology. This state-of-the-art SNIa compilation contains 1048 points up to
$z\simeq 2.3$. We generate a Gaussian white noise using the Pantheon
covariance matrix, and add it to the fiducial cosmology at the redshifts of
the actual Pantheon data points. This effectively creates a “Pantheon-like”
mock dataset with a Gaussian noise. Since we aim to test the performance of
GPR, we choose as a fiducial cosmology a non-$\Lambda\text{CDM}$ dark energy
with $\Omega_{\text{m0}}=0.3$, in our case, the phenomenologically emergent
dark energy (PEDE) model [16]. This model has the same number of degrees of
freedom as $\Lambda\text{CDM}$, and has a phantom equation of state in the
asymptotic past and transition to a cosmological constant ($w=-1$) in the
asymptotic future.

We compare PEDE to a commonly used parameterization of the DE equation of
state, the Chevallier-Polarski-Linder (CPL) [17, 18] parameterization. The CPL
parameterization is a Taylor expansion of $w$ at $z=0$:

$\displaystyle w(a)$ $\displaystyle=w_{0}+w_{a}(1-a).$ (2.17)

This two-parameter model has been widely used in the literature to describe
the low-redshift evolution of the equation of state.

Fig. 1 shows the equation of state and the distance modulus for the fiducial
(PEDE) model in red. In blue is the best-fit for the CPL case with free
parameters $(\Omega_{m},w_{0},w_{a})$. Although the distance moduli for both
cases are similar (difference of about 0.02 mag at $z=1.5)$, the equation of
states show different evolution. Since $\mu$ is obtained by integrating twice,
the difference between PEDE and the CPL best-fit is attenuated in $\mu$ space.
These results illustrate the fact that the fiducial cosmology is not included
in the CPL parameterization.

Figure 1:  Fiducial cosmology for mock data; PEDE model with
$(h,\Omega_{m})=(0.7,0.3)$. _Left_ : The equation of state of our fiducial
cosmology and the CPL best-fit. _Right_ : Distance moduli of fiducial
cosmology and the CPL best-fit.

##  3 Results

In this section, we study in detail different aspects of GPR. Section 3.1
shows the effect of the mean function, Section 3.2.1 is devoted to the effect
of the kernel, and Section 3.3 compare the prior and posterior approaches.

###  3.1 Role of the mean function

Although the first applications of GPR to cosmology [7, 8, 9, 10] considered
the effect of the mean function, many papers use an arbitrary mean function as
an input such as a zero mean for convenience. Refs. [7, 8, 9] used the prior
approach and modeled the equation of state $w(z)$ and its first integral
$y(z)=\int w(u)/(1+u)du$ as a GP, and fix the mean functions
$\overline{\mu}(z)=-1$ and $\overline{y}(z)=-\ln(1+z)$, corresponding to the
cosmological constant. Ref. [10] first applied an iterative smoothing method
[19, 20, 21] to prewhiten the data and obtain a mean function, then applied
the posterior approach to reconstruct $\mu(z)$ directly from the data. Because
the iterative smoothing approach converges to a unique solution, the choice of
the mean function is thus unique and defined by the data. Refs. [12, 22]
studied the LML distribution depending on the input mean function, and use the
LML as a function of $\sigma_{f}$ to reject an input mean function. Ref. [23]
used the $\Lambda\text{CDM}$ best-fit with a multiplicative parameter as a
mean function in order to give more flexibility than a constant mean function.
While a constant mean function is legitimate in the case of stationary
signals, in the case of a time-evolving quantity such as the distance moduli,
or other cosmological quantities such as the Hubble expansion history, it is
less clear. We will show that the zero mean cause a problem when
reconstructing distance moduli and illustrate how a different choice of mean
function yields very different posteriors. In this section, we will use the
posterior approach.

####  3.1.1 Zero mean function

Figure 2:  GP regression and its derivative from GaPP. The blue solid line is
the reconstructed mean of the GPR, and the blue bands show the 68% and 95%
confidence levels. _Left_ : Reconstruction of the distance moduli $\mu$ with
optimized hyperparameters assuming a zero mean function. _Right_ : GP
regression for the derivative of the distance moduli $\mu^{\prime}$ with GaPP.
This regression is unphysical (below the red line, negative $\mu^{\prime}$)
beyond redshift $z\simeq 1$.

Following Ref. [11], a large number of studies in cosmology have used the GaPP
package to reconstruct cosmological quantities. To illustrate our point, we
will also use GaPP in this section. We verified that our own GP code and GaPP
yield consistent results. Negligible differences arise from the choice of the
optimization algorithm.

The left-hand panel of Fig. 2 shows the difference between the predictive mean
function and the true value to be reconstructed when using zero as the input
mean function, as is done by GaPP. The best-fit hyperparameters for this
realization are $(l,\sigma_{f})=(0.13,37.19)$. While the reconstructed
distance moduli are consistent with the truth, oscillations of the order of
the best-fit hyperparameter $l$ are observed.

The right-hand panel of Fig. 2 shows the reconstructed derivative of the
distance moduli. Oscillations are also visible in the derivative. This is
understandable since the derivative involves a factor of order $l^{-2}$. In
addition, beyond $z\simeq 1$, the predictive mean becomes negative, which is
unphysical, since

$\displaystyle\mu^{\prime}(z)$ $\displaystyle=\frac{5}{\ln
10}\frac{d_{\mathrm{L}}^{\prime}(z)}{d_{\mathrm{L}}(z)}>0.$ (3.1)

We conclude that zero is not a suitable choice of the mean function to
reconstruct the distance moduli.

####  3.1.2 Best-fit as mean-function

Next, we illustrate the effect of the mean function by using an alternative
mean function. We choose a mean function that we expect to describe the data
relatively well, for instance the best-fit $\Lambda\text{CDM}$ case.

(a) zero mean

(b) best-fit $\Lambda\text{CDM}$ mean

Figure 3: LML as a function of the hyperparameters obtained by MCMC sampling. _Left_ : The LML for the zero mean as the input mean function. _Right_ : The LML for the best-fit $\Lambda\text{CDM}$ mean as the input mean function. | $l$ | $\sigma_{f}$ | Note
---|---|---|---
Posterior approach | Zero mean | $\mathcal{U}(0,5)$ | $\mathcal{U}(0,100)$ | Section 3.1.2
Best-fit mean |  $\mathcal{U}(0,5)$111A different prior of $\mathcal{U}(0.05,5)$ is used to reconstruct the derivative of distance modulus in Fig. 4(c) and 4(d). See Section 3.1.2. | $\mathcal{U}(0,5)$ | Section 3.1.2
Full-marginalization | $\mathcal{U}(0,5)$ | $\mathcal{U}(0,5)$ | Section 3.2
Prior approach | Full-marginalization | $\mathcal{U}(0,5)$ | $\mathcal{U}(0,5)$ | Section 3.3
Table 1:  Priors for the hyperparameters and cosmological parameters.
$\mathcal{U}$ denotes a uniform distribution.

The effect of the mean function can be immediately seen in Fig. 3, which shows
the LML as a function of the hyperparameters for both mean functions. The LMLs
are obtained by Markov Chain Monte Carlo (MCMC) algorithm as implemented in
the emcee package [24]. The priors for the hyperparameters are shown in Table
1. As noted by [11], in the zero mean-function case, the LML is sharply peaked
at $(l,\sigma_{f})=(0.13,37.19)$, and therefore it can be well approximated by
a Dirac distribution at this point. However, the situation is very different
when using the best-fit $\Lambda\text{CDM}$ as an input mean function, as can
be seen in the right-hand panel. In this case, since the fit to the data is
relatively correct, the LML is not sharply peaked and prefers a small value of
$\sigma_{f}$. We stress the difference of scales in the $\sigma_{f}$ axis. In
other words, since the input mean function is a “good enough” description of
the data, GP does not find any preference for deviations, and the LML peaks at
small values of $\sigma_{f}$. This is clearly not the case in the zero mean
function case, where the LML shows the need for a deviation from zero.

On the other hand, there is no preference found for $l$ in the LML with the
best-fit mean function as shown in Fig. 3(b). The covariance of the derivative
is roughly proportional to $l^{-2}$ when reconstructing the derivative. Thus,
if a small $l$ is chosen during MCMC sampling, its reconstruction of the
derivative will have large error bars. This will remain even after the
hyperparameters are marginalized over. To prevent too small $l$ being sampled,
a prior with a lower bound for $l$ is chosen as commented on Table 1.

(a) optimized hyperparameters

(b) marginalized hyperparameters

(c) optimized case: derivative of $\mu$

(d) marginalized case: derivative of $\mu$

Figure 4: GP regressions for different hyperparameter selection. _Top-Left_ :
Optimized hyperparameters. _Top-Right_ : Marginalized hyperparameters.
_Bottoms_ : Derivative of $\mu$ with optimized and marginalized
hyperparameters. We stress out that the y-axis in this figure is 10 times
smaller than Fig. 2. Therefore, the unphysical region $(\mu^{\prime}<0)$ is
not shown here.

Fig. 4(a) and 4(c) show the reconstructions of the distance moduli and its
derivative using the optimal hyperparameters assuming the best-fit
$\Lambda\text{CDM}$ as an input mean function. As expected from the LML
distribution in Fig. 3(b), the reconstruction is biased at less than $1\sigma$
towards the input mean function due to preference for small values of
$\sigma_{f}$.

Since no sharp peak can be seen in the LML with the best-fit
$\Lambda\text{CDM}$ mean, many choices of hyperparameters can give reasonable
fits to the data. Therefore, in a Bayesian point of view, the meaningful
approach is thus to marginalize the LML (eq. (2.8)) over the hyperparameter
space.

$p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\mu})=\iint{p(\boldsymbol{y}|\boldsymbol{f},\boldsymbol{X})p(\boldsymbol{f}|\boldsymbol{X},\boldsymbol{\mu},\mathcal{H})\mathrm{d}\boldsymbol{f}\mathrm{d}\mathcal{H}}.$
(3.2)

Fig. 4(b) and 4(d) show the reconstructions for the marginalized case.
Although the reconstructed mean functions are consistent with the optimized
case (Figs. 4(a) and 4(c)), the reconstructions with the optimized
hyperparameters have narrower confidence levels. These results show that the
_optimization_ approach significantly underestimates the errors by not
including other valid reconstructions, which are however included in the
_marginalization_ approach.

###  3.2 Full-marginalization

Figure 5: Full-marginalization with the CPL parameterization. The result is
marginalized over five parameters; two hyperparameters $(l,\sigma_{f})$ and
three cosmological parameters $(\Omega_{m},w_{0},w_{a})$.

In this section, we will introduce a rigorous application of GPR within a full
Bayesian formalism. As we saw previously, marginalizing the LML over the
hyperparameter space is a natural Bayesian concept. However, as we saw in
Section 3.1, the choice of the mean function is also not well justified.
Therefore, we suggest full-marginalization, which is done by integrating the
LML over the space of mean functions as well as hyperparameters:

$p(\boldsymbol{y}|\boldsymbol{X})=\iiint{p(\boldsymbol{y}|\boldsymbol{f},\boldsymbol{X})p(\boldsymbol{f}|\boldsymbol{X},\boldsymbol{\mu},\mathcal{H})d\boldsymbol{f}d\mathcal{H}}d\boldsymbol{\mu}.$
(3.3)

Of course, in practice it is impossible to marginalize over the infinite space
of possible mean functions. Instead, we chose a _reasonable_ family of
functions, the CPL parameterization, to marginalize over. CPL is a reasonable
model to marginalize over, since we expect it to describe the data quite well.
However, as we saw, the fiducial model is not included in the CPL family.

We used the emcee package to marginalize over the cosmological parameters
$(\Omega_{\text{m0}},w_{0},w_{a})$ and hyperparameters. This idea is similar
to the recent approach by Ref. [23], who use the best-fit $\Lambda\text{CDM}$
as a mean function and marginalize over a multiplicative parameter.

Fig. 5 shows the result with full-marginalization for one realization. Since
the reconstruction with model-independent approach relies on a particular
random realization, the result of full-marginalization has slight deviations
from the truth at low redshift. We will discuss the bias of the methods we
performed in section 3.2.2. Appendix B shows a hundred realizations.

####  3.2.1 Role of the kernel function

(a) Matérn kernel for $\nu=3/2$

(b) Matérn kernel for $\nu=9/2$

Figure 6: Full-marginalization results with different kernels in the Matérn
class.

In this section, we test the role of the kernel function. Another common
kernel is the Matérn class of kernel functions given by

$k(x,x^{\prime})=\sigma_{f}^{2}\frac{2^{1-\nu}}{\Gamma\left(\nu\right)}\left(\frac{\sqrt{2\nu}\left(x-x^{\prime}\right)}{l}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}\left(x-x^{\prime}\right)}{l}\right)$
(3.4)

where $\nu>0$ is the smoothness parameter, $l$ and $\sigma_{f}$ are
hyperparameters that play similar roles to the squared exponential case. GP
with the Matérn kernels allows to reconstruct up to the $n$th derivative of a
function for $\nu>n$. The $\nu\rightarrow\infty$ limit is the squared
exponential function. Ref. [25] argued for the use of the Matérn kernels to
reconstruct the equation of state, the comoving distance and its derivative on
mock data, as the squared exponential kernel may be deemed too smooth. Fig. 6
shows reconstructions with Matérn kernels for $\nu=3/2$ and $9/2$. The
reconstructions are essentially indistinguishable from the SE case, further
demonstrating the importance of the mean function rather than the kernel on
the reconstruction. The $\nu=3/2$ case shows more oscillations, as expected
from the less smooth kernel.

####  3.2.2 Bias of the method

(a) zero mean, optimized hyperparameters

(b) best-fit $\Lambda\text{CDM}$ mean, marginalized hyperparameters

(c) full-marginalization

Figure 7: Averaged results for 500 realizations for three different cases.
_Top-left_ : Zero mean function, optimized hyperparameters. _Top-right_ : The
best-fit $\Lambda\text{CDM}$ mean function, marginalized hyperparameters.
_Bottom_ : Full-marginalization with the CPL parameterization.

Fig. 7 shows the average of 500 realizations for three different choices of
mean functions and hyperparameter selection. The top-left panel shows the
averaged result for the zero-mean function with optimiziation of the
hyperparameters (i.e., corresponding to the GaPP approach). Clearly,
unphysical oscillations are present at $z\geq 1$, where the data become
scarce. As seen in Section 3.1.1, this leads to an unphysical derivative. The
top-right panel shows the case where the mean function is the best-fit
$\Lambda\text{CDM}$ and we marginalized over the hyperparameters. While the
behaviour is smoother than the zero-mean case, there is a small bias towards
the input mean function (in red dashedline). The bottom panel shows our
suggested improvement, by marginalizing over the hyperparameters and a family
of mean functions, in this case CPL. Even though the fiducial distance moduli
(PEDE model) are not included in the family of input mean functions, the
posteriors clearly show that our reconstructions are unbiased.

###  3.3 Prior versus Posterior approach

(a) posterior approach

(b) prior approach

Figure 8: Reconstruction of the distance moduli $\mu$ for marginalized
hyperparameters with the posterior and prior approaches. _Left_ : Posterior
approach. _Right_ : Prior approach.

In this section, we compare the prior and posterior approaches. We use the
best-fit $\Lambda\text{CDM}$ mean function and marginalization over the
hyperparameters as in the Section 3.1.2. Fig. 8 shows the reconstructions for
the prior and posterior approaches. Both cases give consistent behaviors, as
expected mathematically.

However, the prior approach fails to reconstruct the distance moduli from a
zero input mean function. Since the distance modulus increases rapidly from
$z\simeq 0$ and then slowly increases at high redshift, it is difficult for
the GP to reconstruct the distance modulus. To be specific, reconstruction
oscillates for the correlation length $(l)$. If a reconstruction increases
sharply near $z\simeq 0$, it should then decrease sharply again, which would
give a bad fit to the distance moduli data. In other words, the reconstruction
of the posterior approach with the zero input mean function is almost
impossible in practice. Therefore, we suggest giving a reasonable input mean
function for the GPR.

##  4 Discussion and conclusion

We investigated the effects of the choice of the mean function, the
hyperparameter selection in GP regression. Our findings are summarized below.

  * •

The choice of the mean function is important and should carefully be selected.
Most of the literature so far fails to consider this, with the notable
exceptions of [8, 7, 9, 10].

  * •

Using a zero mean function on a nonstationary signal results in an artificial
sharp peak in the LML, and by fixing the hyperparameters that maximize the
LML, one miss many combinations that still give a good fit to the data. We
advise against the use of zero as a mean function when the quantity is not
stationary, as is the case for most cosmological quantities.

  * •

We propose full-marginalization: using a family of mean functions, i.e., non-
flat $\Lambda$CDM or CPL, and marginalizing over the parameters associated
with this mean function and the hyperparameters. This reduces the dependence
on the input mean function. This approach is thus “more” model-independent,
albeit still not exhaustive.

  * •

Our method is unbiased, and allows to reconstruct the true mean function
without any model assumption, as shown by Fig. 7.

  * •

We verified the equivalence between the “prior” and “posterior” approaches,
which we expect from theory. The “prior” approach can be useful when fitting
different data that are not linearly dependent, and for forward-modeling.

We showed that arbitrarily setting the mean function to zero is not a neutral
choice, and impact the hyperparameter selection. Fig. 3 shows that the
observation by Ref. [11] that the marginal likelihood is sharply peaked, and
that one may just use the best-fit hyperparameters, is actually an artifact of
the choice of the mean function. If one chooses a mean function closer to the
truth, such as the $\Lambda\text{CDM}$ best-fit as shown in Fig. 3(b), the LML
peaks at low values of $\sigma_{f}$, independently of $l$ and the choice of
the hyperparameters would not be so obvious. This is because GP is very
efficient at picking deviations from the mean function. If the choice of the
mean function is a close description of the data, the data show no preference
for any deviations, and the LML peaks at $\sigma_{f}\simeq 0$. Marginalizing
over the hyperparameters becomes the obvious choice in the Bayesian point of
view.

However, choosing the best-fit $\Lambda\text{CDM}$ model as the mean function
also introduces a degree of arbitrariness, and the method is less model-
independent. The reconstructions are biased towards the mean function, as
shown in the top-right panel of Fig. 7.

In practice, we suggest as an alternative to marginalize over a family of mean
functions. We show the example of the CPL family, which yields unbiased
reconstructions, despite the fact that the data were generated using a
fiducial cosmology that is not included in CPL. As shown in Appendix A, the
results are independent of the family of mean functions, at least as long as a
“reasonable” choice of mean functions is made.

In this work, we marginalize over the CPL family of mean functions, while the
data were generated using the PEDE model. This is a reasonable choice, since
CPL contains the concordance $\Lambda\text{CDM}$ model. Should the data be too
different from a CPL, the choice of mean function should be adjusted
accordingly. Appendix A shows that the method is independent of the choice of
the family of mean functions.

Therefore, we recommend researchers in the field not to blindly use existing
packages that assume a zero mean function and fix the value of the
hyperparameters to the best-fit of the LML, but instead to use our full-
marginalization approach. Given the popularity of GP in cosmology, we hope
that the community will adopt these practices. We note that in practice, GaPP
has been applied to other quantities than the distance modulus, in particular,
direct measurements of the expansion history via cosmic chronometers.

## Acknowledgments

S.-g. H. and B. L. acknowledge the support of the National Research Foundation
of Korea (NRF-2019R1I1A1A01063740). B. L. and A. S. acknowledge the support of
the Korea Institute for Advanced Study (KIAS) grant funded by the government
of Korea. A. S. would like to acknowledge the support by National Research
Foundation of Korea NRF-2021M3F7A1082053.

## References

  * [1] A.G. Riess, A.V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P.M. Garnavich et al., _Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant_ , _AJ_ 116 (1998) 1009 [astro-ph/9805201].
  * [2] S. Perlmutter, G. Aldering, G. Goldhaber, R.A. Knop, P. Nugent, P.G. Castro et al., _Measurements of $\Omega$ and $\Lambda$ from 42 High-Redshift Supernovae_, _ApJ_ 517 (1999) 565 [astro-ph/9812133].
  * [3] A.G. Riess, S. Casertano, W. Yuan, L. Macri, B. Bucciarelli, M.G. Lattanzi et al., _Milky Way Cepheid Standards for Measuring Cosmic Distances and Application to Gaia DR2: Implications for the Hubble Constant_ , _ApJ_ 861 (2018) 126 [1804.10655].
  * [4] A.G. Riess, W. Yuan, L.M. Macri, D. Scolnic, D. Brout, S. Casertano et al., _A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km/s/Mpc Uncertainty from the Hubble Space Telescope and the SH0ES Team_ , _arXiv e-prints_ (2021) arXiv:2112.04510 [2112.04510].
  * [5] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., _Planck 2018 results. VI. Cosmological parameters_ , _A &A_ 641 (2020) A6 [1807.06209].
  * [6] C.E. Rasmussen and C.K.I. Williams, _Gaussian Processes for Machine Learning_ , MIT Press (2006).
  * [7] T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib et al., _Nonparametric reconstruction of the dark energy equation of state_ , _Phys. Rev. D_ 82 (2010) 103502 [1009.5443].
  * [8] T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib et al., _Nonparametric Dark Energy Reconstruction from Supernova Data_ , _Phys. Rev. Lett._ 105 (2010) 241302 [1011.3079].
  * [9] T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib et al., _Nonparametric reconstruction of the dark energy equation of state from diverse data sets_ , _Phys. Rev. D_ 84 (2011) 083501 [1104.2041].
  * [10] A. Shafieloo, A.G. Kim and E.V. Linder, _Gaussian process cosmography_ , _Phys. Rev. D_ 85 (2012) 123530 [1204.2272].
  * [11] M. Seikel, C. Clarkson and M. Smith, _Reconstruction of dark energy and expansion dynamics using Gaussian processes_ , _J. Cosmology Astropart. Phys_ 2012 (2012) 036 [1204.2832].
  * [12] R.E. Keeley, A. Shafieloo, B. L’Huillier and E.V. Linder, _Debiasing cosmic gravitational wave sirens_ , _MNRAS_ 491 (2020) 3983 [1905.10216].
  * [13] B. L’Huillier, A. Shafieloo, D. Polarski and A.A. Starobinsky, _Defying the laws of gravity I: model-independent reconstruction of the Universe expansion from growth data_ , _MNRAS_ 494 (2020) 819 [1906.05991].
  * [14] S. Joudaki, M. Kaplinghat, R. Keeley and D. Kirkby, _Model independent inference of the expansion history and implications for the growth of structure_ , _Phys. Rev. D_ 97 (2018) 123501 [1710.04236].
  * [15] D.M. Scolnic, D.O. Jones, A. Rest, Y.C. Pan, R. Chornock, R.J. Foley et al., _The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample_ , _ApJ_ 859 (2018) 101 [1710.00845].
  * [16] X. Li and A. Shafieloo, _A Simple Phenomenological Emergent Dark Energy Model can Resolve the Hubble Tension_ , _ApJ_ 883 (2019) L3 [1906.08275].
  * [17] M. Chevallier and D. Polarski, _Accelerating Universes with Scaling Dark Matter_ , _International Journal of Modern Physics D_ 10 (2001) 213 [gr-qc/0009008].
  * [18] E.V. Linder, _Exploring the Expansion History of the Universe_ , _Phys. Rev. Lett._ 90 (2003) 091301 [astro-ph/0208512].
  * [19] A. Shafieloo, U. Alam, V. Sahni and A.A. Starobinsky, _Smoothing Supernova Data to Reconstruct the Expansion History of the Universe and its Age_ , _Mon. Not. Roy. Astron. Soc._ 366 (2006) 1081 [astro-ph/0505329].
  * [20] A. Shafieloo, _Model Independent Reconstruction of the Expansion History of the Universe and the Properties of Dark Energy_ , _Mon. Not. Roy. Astron. Soc._ 380 (2007) 1573 [astro-ph/0703034].
  * [21] A. Shafieloo, B. L’Huillier and A.A. Starobinsky, _Falsifying $\Lambda$CDM: Model-independent tests of the concordance model with eBOSS DR14Q and Pantheon_, _Phys. Rev. D_ 98 (2018) 083526 [1804.04320].
  * [22] A. Krishak and D.K. Hazra, _Gaussian Process Reconstruction of Reionization History_ , _ApJ_ 922 (2021) 95 [2106.01728].
  * [23] J. Ruiz-Zapatero, C. García-García, D. Alonso, P.G. Ferreira and R.D.P. Grumitt, _Model-independent constraints on $\Omega$m and H(z) from the link between geometry and growth_, _MNRAS_ (2022)  [2201.07025].
  * [24] D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, _emcee: The MCMC Hammer_ , _PASP_ 125 (2013) 306 [1202.3665].
  * [25] M. Seikel and C. Clarkson, _Optimising Gaussian processes for reconstructing dark energy dynamics from supernovae_ , _arXiv e-prints_ (2013) arXiv:1311.6678 [1311.6678].

##  Appendix A Dependence on the family of mean function

(a) CPL parameterization $(\Omega_{m},w_{0},w_{a})$

(b) non-flat $\Lambda\text{CDM}$ $(\Omega_{m},\Omega_{\Lambda})$

Figure 9: Averaged full-marginalization results for two families of mean
functions.

In this section, we compare the marginalization over the CPL family of mean
functions to another family, in this case, non-flat $\Lambda\text{CDM}$. The
left and right hand panels of Fig. 9 respectively show the posteriors on the
distance moduli when applying our full-marginalization approach using the CPL
and non-flat $\Lambda\text{CDM}$ families of mean functions. Clearly, the
results do not depend on the choice of the family, as long as a reasonable
choice of mean function is made.

##  Appendix B Individual realizations

Figure 10: Full-marginalization for all realizations with the CPL
parameterization. Figure 11: Full-marginalization for all realizations with
the CPL parameterization.

The figures presented in the body of the paper show the reconstructions for
one particular realization of the mock data. Fig. 10 and 11 shows the
reconstructions for 100 of the 500 realizations used in the full
marginalization case (Fig. 7).