\section{Introduction}

Galaxy surveys have provided observations of  the growth of large-scale structure  with steadily increasing precision and range. This trend is set  to accelerate in the coming years with next-generation spectroscopic surveys such as those with {\em Euclid} \cite{Blanchard:2019oqi}, the Square Kilometre Array (SKA) \cite{Bacon:2018dui}, the Dark Energy Spectroscopic Instrument (DESI) \cite{Levi:2019ggs}, and the Nancy Grace Roman Space Telescope \cite{Dore:2019pld}. The evolution of matter perturbations gives direct insight into the underlying theory of gravity. Whether the correct description is  the standard model of cosmology $\Lambda$CDM, a dynamical dark energy model or modified gravity  (see e.g. \cite{Capozziello:2011et,Clifton:2011jh} for reviews), is still an open question.   

Studying gravity on cosmological scales has been pursued via two paths in recent years. The model-dependent approach constrains a particular model of dark energy or modified gravity with observations and compares it to $\Lambda$CDM. This approach has shown hints that alternative theories may be favoured via Bayesian evidence considerations \cite{Peirone:2019aua,Sola:2019jek,Aoki:2020oqc,DiValentino:2020kpf}. The `agnostic' approach avoids choosing a particular model and develops predictions as model independent as possible, which can be used to evaluate how they differ from those of a given model. 

Machine learning offers several model-agnostic methods, including Gaussian Process (GP) regressions, which enable data-driven reconstructions of trends \cite{Rasmussen}. The characteristic input of GP is statistical and corresponds to a kernel to build a covariance matrix and a mean prior. The use of GP is now common in cosmology and has been focused mainly on the background evolution of the universe \cite{Holsclaw:2010nb,Holsclaw:2010sk,Holsclaw:2011wi,Shafieloo:2012ht,Seikel:2012uu,Yahya:2013xma,Bester:2013fya,Busti:2014dua,Bester:2015gla,Bester:2016fbs,Joudaki:2017zhq,Gomez-Valent:2018hwc,Haridasu:2018gqm,Gerardi:2019obr,Keeley:2019hmw,Bengaly:2019oxx,Bengaly:2019ibu,Aljaf:2020eqh,Keeley:2019hmw,Colgain:2021ngq,Dhawan:2021mel,Mukherjee:2021kcu,2021arXiv210407356V,Canas-Herrera:2021qxs}. GP has also proved useful to determine the value of the Hubble constant using strong lensing data and supernova data \cite{Liao:2019qoc,Liu:2019ddm,Pandey:2019yic,Liao:2020zko,Renzi:2020fnx,Bengaly:2020neu,Zhang:2021tmg}, to improve photometric redshift estimates \cite{2016MNRAS.462..726A}, and to speed up the processing of weak lensing data \cite{Mootoovaloo:2020ott}. Less work has been done with GP on measurements of perturbations. Examples are reconstructions of the linear anisotropic stress  \cite{Pinho:2018unz} and of the growth rate of large-scale structure \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt}. 

In this paper we use GP to reconstruct the growth rate of large-scale structure. We produce forecasts for a nominal Stage IV spectroscopic survey. Our focus is not only on the growth variable $f\sigma_8$, but also on the separated growth rate $f$ and variance of matter fluctuations $\sigma_8$. We apply the concept of multi-task GP \cite{Caruana1998,Rasmussen,866abfcac5f9481d97628a546255878b,21ae58e5ad934d7494d64e5bfd4a1a52,Melkumyan,Vasudevan}, previously used in cosmology to improve reconstructions of the background acceleration \cite{Haridasu:2018gqm}.

The paper is structured as follows. We review the basics of the growth of large-scale structure and the methodology of our analysis in \autoref{sec:context}. We then analyse and compare the single- and multi-task reconstructions of mock data in \autoref{sec:perfo}, where we display the improvement from using multi-task GP. Previous single-task GP studies of growth rate in \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt} do not give much detail about the assumptions specific to their reconstruction of $f\sigma_8$. In \autoref{sec:current}, we discuss difficulties in reproducing such results and how we find current growth data to be ill-suited for robust GP reconstructions. We conclude in \autoref{sec:discussion}.

\section{Context}\label{sec:context}

We briefly review the growth of large-scale structure before presenting the GP reconstructions from the single-task to multi-task approach. Then we show how our mock data are derived and the numerical tools that we employ.

\subsection{Growth of large-scale structure}\label{sec:growthgen}

The galaxy distribution in redshift space is distorted  by galaxy peculiar velocities that are induced by the growth of cosmological structures, an effect known as redshift-space distortions (RSD) \cite{1977ApJ...212L...3S,Kaiser:1987qv,Guzzo:2008ac,Song:2008qt,Percival:2008sh}. RSD measurements determine the product of  the growth rate $f$ and variance of matter fluctuations $\sigma_8^2$, where
\begin{eqnarray}
\label{eq:defD}
    f(z)&=&-\frac{{\rm d}\ln{D(z)}}{{\rm d}\ln(1+z)}\, ,\\
\sigma^2_8(z) &=&\frac{1}{2\pi^2}\int_0^{\infty} P_m(k,z) \, \hat{W}^2(k)\, k^2 \, {\rm d} k\, .
\end{eqnarray} 
Here $D(z)=\delta_m (\bm k,z)/\delta_m(\bm k,0)$ is the linear growth factor of the matter density contrast $\delta_m$, $P_m$ is the matter power spectrum and $\hat{W}$ is the Fourier transform of the top-hat window function that  smooths over spheres of radius $8\,{\rm Mpc}/h$.

Both these quantities depend on the evolution of $\delta_m$, whose redshift dependence and sometimes scale dependence  will differ between models of gravity. The reconstruction of $f$ and $\sigma_8$ is then an ideal tool to distinguish between  different theories. In general relativity with standard dark energy, the linear growth rate does not depend on the scale $k$, but this does not necessarily hold in alternative theories of gravity.

RSD are not sensitive to the two functions $f$ and $\sigma_8$ separately, but to their product $f\sigma_8$. Disentangling $f$ and $\sigma_8$ given a measurement of $f\sigma_8$, requires their degeneracy to be broken.  It is possible to break this degeneracy through two different techniques, both requiring additional observations: (i)~by combining RSD measurements in the power spectrum and bispectrum \cite{Gil-Marin:2016wya}; (ii)~by combining RSD data with galaxy-galaxy lensing data \cite{delaTorre:2016rxm,Shi:2017qpr,Jullo:2019lgq}. 

The use of GP to reconstruct trends from current $f\sigma_8$ data has been studied previously in \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt}. Given the low precision and relative randomness of current data, we find that the reconstructions depend too strongly on the GP assumptions made. We expand on this  in \autoref{sec:current}. 

\subsection{Basics of Gaussian processes}\label{sec:stgp}

Reconstructing a latent function from data with GP essentially corresponds to assuming that the data is a random realisation of the function with Gaussian noise given by the covariance matrix of the data. The reconstruction of $\bm \mathfrak{f}^*$ as a function of $\bm X^*$, given data {$\bm y$} as a function of  $\bm X$, with a covariance matrix $C$, is characterised by the joint probability distribution \cite{Rasmussen}
\begin{equation}\label{eq:jointdistr}
\begin{bmatrix} \bm y \\ \bm \mathfrak{f}^* \end{bmatrix}\! \sim \!
\mathcal{N}\! \left(\! 
\begin{bmatrix}
\bm \mu  \\
\bm \mu^*
\end{bmatrix}
\!,
\begin{bmatrix}
K(\bm X,\bm X) +C & K(\bm X,\bm X^*)\\
K(\bm X^*,\bm X)    & K(\bm X^*,\bm X^*)
\end{bmatrix}\!
\right)\!. 
\end{equation}
Here $\bm \mu$ and $\bm \mu^*$ correspond to  priors on the means of the reconstructed functions  of $\bm X^*$ and $\bm X$ respectively. $K$ is the kernel of the GP, i.e. the covariance function that relates the values of the function to be reconstructed over the points $x$ in $\bm X$ and $x^*$ in $\bm X^*$. From the joint distribution of \autoref{eq:jointdistr}, the mean and covariance of the reconstructed function are 
\begin{align}\label{eq:gpmean}
{\rm mean}(\mathfrak{f}^*)  &=  \bm \mu^* \\\nonumber
& + K(\bm X^*, \bm X)\!\left[\!K(\bm X,\bm X) + C\right]^{-1}\!\! (\bm y -\bm \mu),  \\[2mm]
\label{eq:gpcov}
{\rm cov}(\mathfrak{f}^*)  & = K(\bm X^*,\bm X^*) \\\nonumber
& - K(\bm X^*,\bm X)\! \left[\!K(\bm X,\bm X) + C\right]^{-1} K(\bm X,\bm X^*) .
\end{align}
The marginal log likelihood of the reconstructed function can then be obtained as
\begin{align}\label{eq:gp_log_marginal}
\ln \mathcal{L} = & -\dfrac{1}{2}(\bm y - \bm \mu)^T \left[\! K(\bm X,\bm X) + C \right]^{-1} (\bm y - \bm \mu) \nonumber\\
& - \dfrac{1}{2}\ln\left| K(\bm X,\bm X) + C \right| - \frac{n}{2}\ln 2\pi  \;,
\end{align}
where $n$ is the number of data points. 

It should  now be clear that a GP allows us, given a set of training data points, to reconstruct a cosmological function without the need to assume its trend or a parametrisation, but only relying on the assumption that the data are a Gaussian realisation of the underlying function. For this reason, GP reconstructions fall in the category of non-parametric methods.  However,  this does not imply that GP are independent of  free quantities to be fixed. A kernel must be chosen, which depends on a set of `hyperparameters' that need to be evaluated during the reconstruction process. In addition, the reconstruction also requires a mean prior function.

The kernel most often used in cosmology so far is the squared-exponential
\begin{align}\label{eq:SE}
k(x,x') & = \sigma^2 \exp\Big(-\dfrac{\vert x-x' \vert ^{2}}{2 \xi^{2}}\Big).
\end{align}
Two hyperparameters are related to the characteristic scales in the data:  the correlation length  $\xi$ and signal variance $\sigma$ \cite{Seikel:2013fda}. (See Chapter 4 of \cite{Rasmussen} for more details on covariance functions.)

\subsection{Multi-tasking the growth of structures}\label{sec:mtgp}

The previous section describes how to reconstruct one function from a data set, i.e. a single-task GP. In order to generalise this for several data sets and latent functions which describe the same underlying phenomena, we use a multi-task GP \cite{Caruana1998,Rasmussen}. This extension  properly accommodates possible correlations between the data sets for each function. As described in \autoref{sec:growthgen}, it is straightforward to obtain $f\sigma_8$ measurements, while a more involved procedure is then needed to obtain the disentangled $f$ and $\sigma_8$ contributions. If these data sets originate from different surveys and different methods, the data would in principle not be correlated between components. In general there will be correlations between the data sets. The multi-task GP is able to capture such correlations.

In addition to a joint likelihood that incorporates correlations between the three data sets, multi-task GP also allows us to takes into account that the growth functions are inter-dependent, since all three are probes of the evolution of matter perturbations. Furthermore, one function is a product of the other two.

In the general case, a combined data set including $f$, $\sigma_8$ and $f\sigma_8$ data has a covariance matrix in block form,
\begin{equation}
C=\!
  \left[\! \begin{array}{lll}
   {{\rm cov}(f,f)}     & {\rm cov}(f,\sigma_8)     & {\rm cov}(f,f\sigma_8)    \\
   {\rm cov}(f,\sigma_8)  & {{\rm cov}(\sigma_8,\sigma_8)}  & {\rm cov}(\sigma_8,f\sigma_8) \\
   {\rm cov}(f,f\sigma_8) & {\rm cov}(\sigma_8,f\sigma_8) & {{\rm cov}(f\sigma_8,f\sigma_8)}    
  \end{array} \! \right],
\end{equation}
where diagonal entries are  auto-covariances  and off-diagonal are cross-covariances.  Similarly, the multi-task GP covariance function encodes not only the three auto-covariance functions, but also the cross-covariances arising from their  inter-dependence as theoretical probes of growth. It has the form
 \begin{equation}\label{eq:Ktilde}
\tilde{K} = 
  \left[ {\begin{array}{lll}
   K_{f,f}      & K_{f,\sigma_8}    & K_{f,f\sigma_8}     \\
   K_{f,\sigma_8}  & K_{\sigma_8,\sigma_8}  & K_{\sigma_8,f\sigma_8}  \\
   K_{f,f\sigma_8} & K_{\sigma_8,f\sigma_8} & K_{f\sigma_8,f\sigma_8}       \\
  \end{array} } \right] + C .
\end{equation}
Provided the data covariance matrix $C$ is block diagonal, setting all non-diagonal blocks of $\tilde{K}$ to zero leads back to the standard approach where each latent function for each data set would be reconstructed as a `single-task', i.e. independently of the other sets. In practice, this is indeed equivalent to considering separately each single-task $i$, where the relevant covariance matrix for the GP is $K_i+C_i$, and doing the respective reconstructions. By contrast, writing the full $\tilde{K}$ as in \autoref{eq:Ktilde} allows multiple functions to be reconstructed simultaneously, including the covariance between the different tasks \cite{Haridasu:2018gqm,866abfcac5f9481d97628a546255878b,21ae58e5ad934d7494d64e5bfd4a1a52,Melkumyan,Vasudevan}.

The $\tilde{K}$ kernel of the multi-task GP is constructed as follows. First, the sub-kernels for each diagonal  latent function must be chosen. Then, the corresponding cross-kernels are obtained from the convolution of the basis functions of each  sub-kernel (see for instance \cite{Melkumyan} and Section 2.2 of \cite{Haridasu:2018gqm} for details). The multi-task GP construction implies that the final reconstruction is characterised by a set of hyperparameters brought by the kernels describing the auto-covariance of each latent function, as in the single-task approach, while the cross-terms do not contribute any additional parameter to the reconstruction. For the case of the squared exponential kernel \autoref{eq:SE}, the basis function is defined as:
\begin{equation}
    g(x) = \sigma\left(\frac{2}{\pi \xi^2}\right)^{1/4}\exp{\left(-\frac{x}{\xi}\right)^2}\;,
\end{equation}
and the convolution of two squared exponential kernels is
\begin{equation}
    \label{eqn:ksese}
    k_{1 \times 2}(x,x') = \sigma_1 \sigma_2 \left(\frac{2\, \xi_1 \xi_2 }{\xi_1^2+\xi_2^2}\right)^{\!1/2}\!\! \exp{\left(-\frac{\vert x-x' \vert^2}{\xi_1^2+\xi_2^2}\right)}\!.
\end{equation}

An important issue for multi-task GP is {`negative transfer'} and {`over-fitting'}. As discussed in \cite{Haridasu:2018gqm}, these occur when one data set is much more constraining than the others and so dominates the reconstruction. It is particularly problematic in regions where data are sparse. This is something which will be of importance with future data. For instance, higher quality and quantity data on $f\sigma_8$ is likely to be released sooner than new data on $f$ or $\sigma_8$. We leave the study of these caveats for future work.

\subsection{Mock data and analysis method}\label{sec:models}

The reference model we consider is flat $\Lambda$CDM, using the parameter values from the \emph{Planck} collaboration 2018 release \cite{Aghanim:2018eyx}. We set parameters to the mean values obtained via the robust combination of constraints from\footnote{See \emph{Planck} Legacy Archive: \url{https://pla.esac.esa.int/\#cosmology}.}: cosmic microwave background (CMB) base temperature, polarisation and lensing; baryon acoustic oscillations (BAO);  supernovae type Ia (SNIa). We then use \texttt{CLASS} \cite{2011JCAP...07..034B} to compute the redshift evolution of $f$, $\sigma_8$ and $f\sigma_8$ fiducials, in order to build our mock data sets. We generate simplified mocks as follows: 
\begin{itemize}
\item
Each mock contains 20 measurements uniformly distributed between redshift 0 and 2.
\item
We assume a constant relative error of 1\%, consistent with Stage IV expectations. 
\item
We define a multivariate Gaussian distribution with mean equal to the fiducial values of each growth function, and diagonal covariance  matrix defined by the fiducial values times the relative error. We do not incorporate cross-correlations in the mock data since these are too survey-specific. Mock data are then generated as random realisations of this distribution.
\end{itemize}
\vskip1mm
We optimise the hyperparameters, i.e. find their best fit, by minimising the log marginal likelihood in \autoref{eq:gp_log_marginal} applied to the data set. We use the stochastic minimiser `differential evolution' from the python scipy package\footnote{\url{https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html}} in our GP code\footnote{ \url{https://gitlab.com/perenonlouis/GPreconstructions}}. The mean and  covariance of the final reconstruction are then obtained using this set of parameters in \autoref{eq:gpmean} and \autoref{eq:gpcov}. For the sampling of the posteriors of the hyperparameters via Markov chain Monte Carlo (MCMC) methods, we interface our GP code with \texttt{Cobaya\footnote{\url{https://cobaya.readthedocs.io/en/latest}} }\cite{Torrado:2020dgo}. This allows us to use a Metropolis Hasting sampler \cite{Lewis:2002ah,Lewis:2013hha} and the nested sampler \texttt{Polychord} \cite{Handley:2015fda,2015MNRAS.453.4384H}. We also interfaced the Affine-Invariant Ensembler Sampler \texttt{emcee\footnote{\url{https://emcee.readthedocs.io/en/stable/}} } \cite{ForemanMackey:2012ig}. All our statistical analyses and post processing of our results are done with \texttt{GetDist\footnote{\url{https://getdist.readthedocs.io/en/latest/}}} \cite{Lewis:2019xzd}.

\section{Performance and robustness}\label{sec:perfo}

Here we investigate the performance of multi-task GP compared to single-task and make some tests of the results.

\subsection{`Sharing is caring'}\label{sec:sharing}

\begin{figure}[!]
\begin{center}
\includegraphics[trim=0 77 0 0,clip,width=0.92\linewidth]{Figure_1_zplot_fsigma8}\\
\hskip3mm\includegraphics[trim=0 77 0 0,clip,width=0.88\linewidth]{Figure_1_zplot_f}
\includegraphics[width=0.92\linewidth]{Figure_1_zplot_sigma8}
\caption{Optimised single- (grey) and multi-task (red) reconstructions of $f\sigma_8$, $f$ and $\sigma_8$ with 1\% relative error mocks, using $\Lambda$CDM fiducial constrained by CMB T+P+L, BAO, SNIa (see \autoref{sec:models}), for squared-exponential kernel and null mean prior. Shading indicates 68\% confidence intervals about the mean (solid curve). The multi-task GP reconstructs the three growth functions simultaneously by minimising a common likelihood (see \autoref{sec:mtgp}) -- thus `sharing' information between the three reconstructions. As a result the single- and multi-task predictions do not match exactly: multi-task reconstruction for a given function benefits from the constraints on the other functions, thereby decreasing the error on all the reconstructions. Quantitatively, the multi-task GP provides the best fit to the mock data.}
\label{fig:stgp_vs_mtgp_plotz}
\end{center}
\end{figure}

We use a squared-exponential kernel \autoref{eq:SE} and a null mean prior as the base assumptions for the single- and multi-task GP. Note that the reconstructions are robust under the change of these assumptions for reasonable choices. This is studied in \autoref{sec:assumptions}. 

An overlay of the single- and multi-task reconstructions is shown in \autoref{fig:stgp_vs_mtgp_plotz}, for the generation of a single mock. At this level of precision, the reconstructions are very sensitive to the distribution of the data around the fiducial. We discuss  the deviations of the reconstructions in \autoref{sec:deviations}. We concentrate for now on the differences between single- and multi-task reconstructions.

\autoref{fig:stgp_vs_mtgp_plotz} further shows that multi-task GP increases the precision of reconstructions compared to single-task GP. It also changes the redshift trends relative to the single-task case. Although the error for the multi-task GP is in general reduced, this is not guaranteed at all redshifts. The changes in the redshift trends and error sizes are a clear manifestation of the sharing of information between the reconstructions of functions in the multi-task case. Quantitatively, we find that the multi-task log marginal likelihood is larger than the sum of the single-task ones.

\begin{figure*}[!]
\begin{center}
\includegraphics[scale=0.58]{Figure_2_contours_stgp_fsigma8}
\includegraphics[scale=0.58]{Figure_2_contours_stgp_f}
\includegraphics[scale=0.58]{Figure_2_contours_stgp_sigma8}\vskip4mm
\includegraphics[scale=0.45]{Figure_2_contours_mtgp}
\caption{Posterior distributions from the MCMC sampling of the hyperparameters for the single-task (grey contours) and multi-task (red contours) reconstructions shown in \autoref{fig:stgp_vs_mtgp_plotz}. Dashed lines indicate optimised values. The transfer of information between reconstructions in the multi-task case adds correlations between hyperparameters of each kernel and the hyperparameter posteriors become much wider. It also increases values of all hyperparameters relative to single-task, so that the reconstructions will be more `wiggly', as seen in \autoref{fig:stgp_vs_mtgp_plotz}.}
\label{fig:stgp_vs_mtgp_contours}
\end{center}
\end{figure*}

A complementary approach to optimising the GP is to reconstruct the full posterior distributions of the hyperparameters. This allows in particular the possibility to study the degeneracies between them. Sampling the space of hyperparameters from the log marginal likelihood in \autoref{eq:gp_log_marginal} can be achieved using Monte Carlo Markov Chains (MCMC). The multi-task case is more difficult to sample because of non-global minima. We found that a nested sampler is more robust than a Metropolis-Hasting sampler. For the same reason, optimising the hyperparameters with a stochastic minimiser is safer than any method based on gradient descents or similar. The results of the MCMC sampling corresponding to the GP of \autoref{fig:stgp_vs_mtgp_plotz} are shown in \autoref{fig:stgp_vs_mtgp_contours}. Sampling the hyperparameters also allows us to derive the marginalised reconstruction. We verify that the Gaussianity of the mocks renders each best-fit value of the reconstructions to be very close to the maximum of their marginalised posterior. This implies that the optimised reconstructions are virtually the same as the marginalised ones \cite{Seikel:2012uu}.

\autoref{fig:stgp_vs_mtgp_contours} shows that the 2-dimensional posteriors between two hyperparameters of the same kernel display correlations whatever the approach. The transfer of information between the reconstructions in the multi-task case results in correlations between hyperparameters of each kernel. The hyperparameter posteriors are therefore much wider. This means that previously ruled out values of hyperparameters for one of the functions are now allowed, as long as the others change accordingly. This transfer of information also has the effect of increasing the values of all hyperparameters relative to the single-task case. Interestingly, this implies that one effect of the convolutions in the multi-task kernel is to increase the correlation of the data as seen by the GP. The hyperparameters for each reconstruction are larger and therefore they will be more `wiggly'. This is consistent with what seen in \autoref{fig:stgp_vs_mtgp_plotz}. 

A Monte Carlo generation of $10^4$ mocks confirms that the above conclusions are robust. We find that the distribution of best-fit multi-task log marginal likelihood is higher than the sum of the three single-task ones, as shown in \autoref{fig:stgp_vs_mtgp_chi2s_x}. The multi-task thus performs a better fit to the union of the three growth data sets. This highlights an important point: the multi-task approach reduces the width of the confidence interval around the mean, but it also changes its redshift evolution. The single-task approach reconstructs more faithfully the redshift dependence of each data set -- simply because it is oblivious to the data on the other functions. Since the log marginal likelihood remains better for the multi-task case, this means that the reduction of the confidence interval width of the reconstruction dominates over the change of redshift evolution in the goodness-of-fit.

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.92\linewidth]{Figure_3_lnlike}
\caption{Distribution of best-fit values for the multi-task log marginal likelihood (red) and the sum of single-task log marginal likelihoods (grey) from MC generation of {$10^4$} mocks. The sum of the single-task log marginal likelihoods is smaller than the multi-task log marginal likelihood, indicating that the multi-task approach fits the data better.}
\label{fig:stgp_vs_mtgp_chi2s_x}
\end{center}
\end{figure}

\subsection{Assumptions}\label{sec:assumptions}

The first step to evaluate the robustness of GP reconstructions  is to compare predictions with different kernels. The more dependent the results are on this choice, the less  model-independent is the reconstruction. For the growth of structure, we reconstruct underlying smooth functions. We therefore consider only stationary kernels. For these tests, we also concentrate on the single-task reconstruction of $f\sigma_8$. 

We test all the common stationary and non-periodic kernels: exponential, squared exponential, gamma exponential, the Matern class and rational quadratic (see \cite{Rasmussen}). We find that the exponential, gamma exponential, and Matern with $\nu < 5/2$ produce jagged reconstructions, as shown in \autoref{fig:assumptions}. These are discarded  on the basis of the smoothness requirement. Reconstructions performed with the other kernels all give very similar results. The rational quadratic kernel produces virtually the same reconstruction and goodness of fit as the squared exponential. This kernel can be understood as a scale mixture, i.e. an infinite sum, of squared-exponential kernels with different correlation lengths \cite{Rasmussen}. It also has an extra  hyperparameter. The absence of difference between  reconstructions with these kernels thus highlights that adding hyperparameter freedom does not produce new features. 

The Matern class  reconstructions produce errors that increase slightly as $\nu$ decreases. In the limit $\nu \rightarrow \infty$ the Matern kernel recovers the squared exponential. Not surprisingly, the case $\nu=9/2$ is very similar to the squared exponential, as illustrated in \autoref{fig:assumptions}. On the other hand, the smaller is $\nu$ the more sensitive will this kernel be to the noise of the data. This is why a reconstruction with $\nu=3/2$ shows short wiggles on top of the evolution reconstructed by the squared exponential (see \autoref{fig:assumptions}). 
 
In conclusion, from a GP point of view the simplest kernel, the squared exponential, captures all the information about the smooth evolution underlying the data. We  do not find any need  to consider combinations of kernels.  

The next step is to look at the other input defining a GP, i.e. its mean prior. The influence of the mean prior becomes important when the latter is either far away from the data or when the data displays gaps in redshift. The former case should not occur here, given reasonable user choices. The latter could be relevant in the future, depending on survey choices. To see how sensitive the input is with our mock set-up, we compare two common choices of mean prior, i.e. equal to zero or to the mean of the data, and an unreasonable choice, equal to 20. The results are shown in \autoref{fig:assumptions}. 

The unreasonable choice displays large and noisy errors. Nonetheless, the precision of the data is such that the reconstruction does not go to the mean prior in between any neighbouring measurements,  even though the mean prior is chosen significantly far away from the data. For mean prior equal to the mean of the data, the maximum log marginal likelihood is slightly enhanced relative to the null mean prior. However, the reconstruction $\chi^2$ is slightly increased. For real future data, it could be more meaningful to use the mean of the data as mean prior. 

Alternatively, one could input the standard model prediction as mean prior. Here, this is exactly the fiducial, and it produces a flat reconstruction around the fiducial, i.e. with correlation length hyperparameter $\xi$ diverging to infinity. Overall, we avoid a mean prior that is dependent on each mock, to ensure that each mock is treated on the same ground during the MC generations. 
We conclude that the null mean prior is a sound benchmark.
\begin{figure}[!h]
\begin{center}
\includegraphics[trim=0 77 0 0,clip,width=0.9\linewidth]{Figure_7_test_kernel_fsigma8}\\
\includegraphics[trim=0 0 0 0,clip,width=0.9\linewidth]{Figure_7_test_mean_fsigma8}
\caption{Optimised single-task reconstructions of $f\sigma_8$ with 1\% relative error mock, $\Lambda$CDM fiducial constrained by CMB T+P+L, BAO, SNIa. Shading indicates 68\% confidence intervals about the mean (solid curves). Reconstructions are shown with different kernels (top) and different mean priors (bottom). The Matern kernel with $\nu =3/2$ produces jagged reconstructions, while squared exponential and Matern $\nu =9/2$ give very similar results. 
The unreasonable choice of mean prior shows large and noisy errors, while sensible mean priors  produce similar reconstructions.}
\label{fig:assumptions}
\end{center}
\end{figure}

Regarding the priors on the range of the hyperparameters, we choose them large enough that they have no effect on the reconstructions. We thus ensure that the optimised value of a hyperparameter is not  its maximum or minimum allowed. Note that these considerations do not apply when using currently available growth data (see \autoref{sec:current}), highlighting how the precision of observations is a crucial ingredient to apply machine learning tools to cosmological analyses. 

We also verified that marginalising over the hyperparameters produces virtually the same reconstructions as when optimised. This is expected, since the precision of the mock data and its Gaussian distribution around the fiducial produces posteriors that are very close to Gaussian, as shown in \autoref{fig:stgp_vs_mtgp_contours}. 

In conclusion, this analysis shows that at this level of precision of the mock data, the reconstructions are fully driven by data, indicating that our comparisons of the single- and multi-task approaches are robust. 

\subsection{Bias and deviations}\label{sec:deviations}

\begin{figure}[!]
\begin{center}
\includegraphics[trim=0 0 0 00,clip,width=0.9\linewidth]{Figure_4_average_fsigma8}
\caption{Average distribution over $10^4$ mocks of the $f\sigma_8$ reconstructions, computed by drawing 1000 realisations of each optimised GP, from the multivariate Gaussian built with their mean and covariance. We concatenate all 1000$\times${$10^4$} draws. The average fluctuates around the fiducial function, which is well recovered within the error.}
\label{fig:bias}
\end{center}
\end{figure}

Reconstructions can deviate quite significantly from their fiducials given the precision and distribution of the mock data. How likely are these departures and are any biases hidden? We run a series of tests to answer these questions.

First, we check for any overall bias in our reconstruction by performing the equivalent of an `ensemble average' of the reconstructed functions, over the $10^4$ mocks. To estimate the average distribution analytically would require the computation of a mixture of multivariate Gaussian distributions (i.e.\ the outputs from each GP optimisation) which has no closed form. Thus, we compute this average distribution numerically by drawing 1000 realisations of each optimised GP -- from the multivariate Gaussian built with their mean and covariance -- and we concatenate all 1000$\times${$10^4$} draws. We then derive the overall mean and the 68\% limits of the resulting distribution. The results for $f\sigma_8$ in are displayed \autoref{fig:bias}. We note that the average fluctuates around the fiducial function -- but nevertheless the fiducial is well recovered within the error and no significant biases can be seen. We also verified that the distributions at a given redshift are Gaussian. The same conclusions hold for the other two functions.

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.68\linewidth]{Figure_5_chi2}
\includegraphics[width=0.68\linewidth]{Figure_5_chi2_50_data}
\includegraphics[width=0.68\linewidth]{Figure_5_chi2_5_percent}
\caption{Distributions of the reduced $\chi^2/n_{\rm dof}$ (see \autoref{eq:chi2}) of the combined $f\sigma_8\;, f\;, \sigma_8$ reconstructions for the $10^4$ mocks. We vary some characteristics of the mocks: 20 points uniformly distributed between redshift 0 and 2 with 1\% relative error (top); increased number of data points to 50 (middle); decreased precision of the data to 5\% for 20 points (bottom). The distributions of the reduced $\chi^2/n_{\rm dof}$ change depending on the configurations. This highlights that the sensitivity of GP's to the probability distribution of the data depends on the method -- single-task or multi-task -- and on the characteristics of the mock data itself.}
\label{fig:stgp_vs_mtgp_chi2s}
\end{center}
\end{figure}

The second check is to quantify how significantly  the reconstructions deviate from the fiducial. In order to do so, we exploit the fact that we know the true underlying model in our case study, i.e. the fiducial. We define a hybrid $\chi^2$ between the optimised reconstructions and the fiducial as 
\begin{equation}\label{eq:chi2}
    \chi^2 = \sum_{g}\sum_{ij} \left(Y^i_{\rm fid}-Y^i_{g}\right) C^{-1}_{ij,{g}}\left(Y^j_{\rm fid}-Y^j_{g}\right)\, ,
\end{equation}
summed over the three growth reconstructions $g=\lbrace{f,\, \sigma_8,\, f\sigma_8 \rbrace}$. For each reconstruction, the $\chi^2$ is computed for a vector $\bm{Y}$ of 100 redshifts uniformly distributed between 0 and 2. $\bm{C}_{g}^{-1}$ is the corresponding inverse diagonal covariance matrix of the reconstruction. The number of degrees of freedom is therefore $n_{\rm dof}=3 \times 100$. The results are displayed in \autoref{fig:stgp_vs_mtgp_chi2s}.

We find that the reduced $\chi^2$ of the GP reconstructions is not distributed around unity (see top plot). This shows that given the characteristics of the mocks -- relative error, Gaussian probability distribution, number of data -- the reconstructions are too sensitive to the probability distribution of the data and do not recover the fiducial efficiently. In addition, the multi-task distribution is shifted to higher values and is flattened relative to the single-task. This is the result of the tighter confidence intervals produced by the former and not to any particular bias, since we found the reconstructions to be Gaussian distributed on average around the fiducial.

The third check is to test how modifying the characteristics of the mocks influences the previous conclusions. On the one hand, adding data points in the same redshift range (\autoref{fig:stgp_vs_mtgp_chi2s}, middle) does not affect significantly the distribution of reduced $\chi^2$ for the multi-task case. Interestingly, the reduced $\chi^2$ of the single-case is pulled towards reproducing that of multi-task. On the other hand, lowering the relative error of the data shifts the distributions toward a better fiducial recovery.

The fourth test we conduct is to analyse quantitatively how much  the reconstructions deviate from the fiducial. To this end, we compute for each mock the fraction of the redshift range over which the fiducial falls outside of the 95\% confidence interval of the reconstructions \cite{Seikel:2012uu}. We display the distribution of this fraction over the $10^4$ mocks for our three reconstructed functions in \autoref{fig:stgp_vs_mtgp_counts}. The results show a trend similar to our previous $\chi^2$ test -- i.e., the tighter multi-task errors exclude the fiducial function more often. Another important result highlighted by these plots is that the distributions depend significantly on the function considered, and hence on the characteristics of the mocks. For instance, the redshift evolution of the function has an importance since we consider relative errors to construct the mocks.

The tests in this section highlight some of the subtleties associated with GP regressions. Even when the reconstruction is not biased on average, the fact that it performs `too well' in some cases (i.e. has stringent error bars) can complicate the assessment of whether it excludes a given cosmological model or not. Moreover, as we showed above, such effects can depend on the specific choice of mocks, hence ultimately on the intrinsic properties of the data (noise amplitude, redshift dependence, etc.) to which we apply the GP reconstruction. A full diagnostic of such effects will be the topic of an in-depth future work.

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.7\linewidth]{Figure_6_counts_fsigma8}
\includegraphics[width=0.7\linewidth]{Figure_6_counts_f}
\includegraphics[width=0.7\linewidth]{Figure_6_counts_sigma8}
\caption{Distribution of the fraction of the redshift range over which the fiducial falls outside of the 95\% confidence interval of the reconstructions. The distribution for each  reconstruction is computed over the $10^4$ mocks for 100 points uniformly distributed between redshift 0 and 2. The distributions depend significantly the function considered.}
\label{fig:stgp_vs_mtgp_counts}
\end{center}
\end{figure}


\section{Applications to current data}\label{sec:current}

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.9\linewidth]{Figure_8_current_zplot_fsigma8}
\includegraphics[width=0.9\linewidth]{Figure_8_current_contours_fsigma8}
\caption{{\it Top:} Optimised single-task reconstructions of current $f\sigma_8$ data with squared-exponential kernel and null mean prior. Upper prior bounds on correlation length hyperparameter $\xi_{f\sigma_8}$ are either restricted the maximum redshift range of the data (grey) or unbounded (red). Shading indicates 68\% confidence intervals. Dashed curve is the $\Lambda$CDM fiducial. {\it Bottom:} The corresponding MCMC sampling of posterior distributions on the hyperparameters. Shading shows 68\% and 95\% confidence intervals. Dashed lines indicate optimised values of hyperparameters. If $\xi_{f\sigma_8}$ is unbounded from above then its best fit wants to be infinite which produces a flat reconstruction; in practice it will be very close to its maximum range prior.}
\label{fig:reco_current}
\end{center}
\end{figure}

To conclude our analysis we investigate reconstructions of current growth data. We consider first the compilation of measurements of $f\sigma_8$ in \cite{Perenon:2019dpc} (updated using \cite{Nesseris:2017vor}). By using this data from different and non-overlapping galaxy surveys, we avoid possible double counting due to overlapping areas or different releases of the same survey. (This differs from the approach of \cite{Yin:2018mvu,Benisty:2020kdt,Li:2019nux}, where an extended set of RSD measurements is considered.) 

The result (with squared-exponential kernel) is shown in \autoref{fig:reco_current}. It is apparent that current RSD data does not allow for a robust GP reconstruction, which depends on the allowed range for the correlation hyperparameter $\xi_{f\sigma_8}$ (red and grey shadings). This dependence persists for different choices of mean prior, including the standard model, and for different kernels.

Allowing the range of $\xi_{f\sigma_8}$ to be very large leads to a posterior that is unbounded from above when explored via MCMC, as shown in the bottom panel. The best-fit value for $\xi_{f\sigma_8}$ tends to infinity and in practice will always be close to the maximum range allowed. This is why the reconstruction is flat (red shading). The error on the reconstruction is then carried solely by the signal variance $\sigma_{f\sigma_8}$, whose 1-dimensional posterior distribution is on the contrary well defined.

This illustrates that current RSD data is simply too noisy for a GP reconstruction independent of assumptions. The data is effectively seen as random noise by the GP.

If we use the extended data set of \cite{Yin:2018mvu,Benisty:2020kdt,Li:2019nux}, the same problem arises. The only way we can reproduce a growth reconstruction similar to \cite{Zhang:2018gjb,Yin:2018mvu,Li:2019nux,Benisty:2020kdt} is to require the range of $\xi_{f\sigma_8}$ to be smaller than the redshift range of the data. This forces the best-fit value of $\xi_{f\sigma_8}$ to be that of the chosen maximum range, as shown in the grey shadings in \autoref{fig:reco_current}.
 
The ratio of the average size of data errors over the total variation of the function is a probe of how faithful the reconstruction is to the fine features of the data. We find that the smaller this ratio is, the more the reconstruction is faithful to data features, independent of priors  and with well-defined posteriors for the hyperparameters. For current data this ratio is much larger than for the mock Stage IV data. If we multiply the current RSD measurements by a factor, or equivalently reduce their covariance by the same factor squared, the reconstruction improves.

In conclusion, we find that the distribution and precision of current data produces GP reconstructions that are strongly dominated by the choices of priors, either the range of the hyperparameters or the mean. We tested whether a multi-task GP adding the few measurements of $f$ and $\sigma_8$ ameliorates the situation. The low number of data points and their precision does not allow us to evade the prior dominance. This scenario is also a typical example of an over-fitting configuration and should not be given much significance.

\section{Conclusion}\label{sec:discussion}

In this paper we explored Gaussian processes as a crucial tool for `agnostic' probes of gravity via the growth of large-scale structure. In particular, we focused on how GP can reconstruct the growth functions, namely the growth rate $f(z)$, the rms of matter fluctuations $\sigma_8(z)$ and their product $(f\sigma_8)(z)$.

We generated simulated measurements for the growth functions using a nominal sensitivity consistent with upcoming Stage IV surveys, and applied the GP pipeline to assess its ability to reconstruct the functions. Crucially, we reconstructed these functions not only separately, using single-task GP, but also via the multi-task GR approach. In multi-task GP, the three functions are reconstructed simultaneously, with a single likelihood. The kernel for the joint reconstruction contains sub-kernels for each function but also cross-convolution kernels that encode the links between the reconstructions. 

A multi-task approach is more faithful than single-task, since a common likelihood allows us to use a data covariance matrix which describes fully the data, i.e. it can take into account the possible correlations in the three data sets. At a more fundamental level, multi-task GP also incorporates the inter-dependence of the growth functions -- all three are probes of the evolution of matter perturbations. Furthermore, one function is the product of the other two.

We showed how the sharing of information between the reconstructions of each growth function in the multi-task approach is not only necessary but also enhances their precision. By performing a Monte Carlo generation of a large number of mocks, we show statistically that, relative to the single-task, the multi-task approach:

\begin{itemize}
    \item 
fits the data better; 
    \item 
is more precise; 
\item 
is more sensitive to the probability distribution of the data and thus can show more deviations from the fiducial. 
\end{itemize}

We verified the robustness of our findings against the change of the assumptions defining GP reconstructions.

Robustness however applies only to the simulated data sets of a nominal future survey. Indeed, we found that when applying our pipeline to currently available data, the reconstruction is very sensitive to the assumptions on the GP settings: priors on the hyperparameter ranges and the mean choices have a significant impact on the final reconstruction. We found that using our baseline settings (a zero prior mean and a very extended allowed range for the hyperparameters), the GP reconstructs a flat function. In order to reconstruct functions that behave as expected, one would instead need to fine tune the hyperparameters. This behaviour is due to the noise of current data being too random.

Next-generation surveys will reach a sensitivity that will enable the robust reconstruction of the growth functions $f(z)$, $\sigma_8(z)$ and $(f\sigma_8)(z)$, using model-agnostic methods such as GP. High-precision measurements will allow us to avoid fine tuning the parameters of the reconstruction, using only limited assumptions. The data-driven model-agnostic reconstruction of the three growth functions is further improved by a multi-task approach that incorporates the inter-dependence of these functions and possible correlations between the three data sets. A multi-task GP avoids the theoretical bias inherent in a single-task approach, which neglects inter-dependence and correlations. 

The ability of the GP to distinguish between a `true' deviation from a given model or an artefact from the noise of the data must be assessed with care. We showed that the accuracy of the reconstructions is highly dependent on the characteristics of the mock data, the function reconstructed and the GP method. Detecting a possible departure from the standard model with a GP reconstruction will require a more refined GP modelling in the future.


\section*{CRediT authorship contribution statement}

\textbf{Louis Perenon}: Conceptualization, Methodology, Software, Formal analysis, Investigation, Writing - Original Draft, Writing - Review \& Editing, Visualization, Supervision, Project administration. \textbf{Matteo Martinelli}: Conceptualization, Methodology, Software, Validation, Formal analysis, Writing - Original Draft, Writing - Review \& Editing. \textbf{St\'ephane Ili\'c}: Methodology, Validation, Writing - Review \& Editing. \textbf{Roy Maartens}: Conceptualization, Methodology, Writing - Original Draft, Writing - Review \& Editing, Supervision. \textbf{Michelle Lochner}: Methodology, Validation, Writing - Review \& Editing. \textbf{Chris Clarkson}: Conceptualization, Validation, Writing - Review \& Editing.  

\section*{Declaration of competing interest}

The authors declare that they have no known competing financial interests or personal relationships that
could have appeared to influence the work reported
in this paper.

\section*{Acknowledgements}

We thank Sambatra Andrianomena, Mario Ballardini, Stefano Camera, Ed Elson, Julien Larena, Marco Raveri and Mario Santos for useful discussions. We thank the authors of \cite{Shi:2017qpr,Jullo:2019lgq} for making their data covariance matrix available to us. We acknowledge the Sciama High Performance Compute Cluster, which is supported by the ICG at the University of Portsmouth, and the Centre for High Performance Computing, South Africa (under the project Cosmology with Radio Telescopes, ASTRO-0945), for providing computational resources for this research project. LP and RM are supported by the South African Radio Astronomy Observatory and the National Research Foundation (Grant No. 75415). MM is supported by a fellowship (code LCF/BQ/PI19/11690015) from `la Caixa' Foundation (ID 100010434) and by the Spanish Agencia Estatal de Investigacion through the grant IFT Centro de Excelencia Severo Ochoa SEV-2016-0597. RM is also supported by the United Kingdom Science \& Technology Facilities Council (STFC) Consolidated Grant ST/S000550/1. CC is supported by STFC Consolidated Grant ST/P000592/1.

\bibliographystyle{elsarticle-num-names}
