# LimberJack.jl: auto-differentiable methods for angular power spectra
analyses

J. Ruiz-Zapatero1,∗ D. Alonso1 C. García-García1 A. Nicola2 A. Mootoovaloo1 J.
M. Sullivan3,4 M. Bonici5 P. G. Ferreira1 ∗jaime.ruiz-
<EMAIL_ADDRESS>1Astrophysics, University of Oxford, DWB, Keble
Road, Oxford OX1 3RH, UK  2Argelander Institut für Astronomie, Universität
Bonn, Auf dem Hügel 71, 53121 Bonn, Germany  3Department of Astronomy,
University of California, Berkeley, CA 94720, USA  4Berkeley Center for
Cosmological Physics, University of California, Berkeley, CA 94720,USA  5INAF-
IASF, Milano, Italy

###### Abstract

We present LimberJack.jl, a fully auto-differentiable code for cosmological
analyses of 2 point auto- and cross-correlation measurements from galaxy
clustering, CMB lensing and weak lensing data written in Julia. Using Julia’s
auto-differentiation ecosystem, LimberJack.jl can obtain gradients for its
outputs an order of magnitude faster than traditional finite difference
methods. This makes LimberJack.jl greatly synergistic with gradient-based
sampling methods, such as Hamiltonian Monte Carlo, capable of efficiently
exploring parameter spaces with hundreds of dimensions. We first prove
LimberJack.jl’s reliability by reanalysing the DES Y1 3$\times$2-point data.
We then showcase its capabilities by using a O(100) parameters Gaussian
Process to reconstruct the cosmic growth from a combination of DES Y1 galaxy
clustering and weak lensing data, eBOSS QSO’s, CMB lensing and redshift-space
distortions. Our Gaussian process reconstruction of the growth factor is
statistically consistent with the $\Lambda$CDM Planck 2018 prediction at all
redshifts. Moreover, we show that the addition of RSD data is extremely
beneficial to this type of analysis, reducing the uncertainty in the
reconstructed growth factor by $20\%$ on average across redshift.
LimberJack.jl is a fully open-source project available on Julia’s general
repository of packages and GitHub .

###### keywords:

Large scale structure, cosmology, statistical inference

##  1 Introduction

Cosmology is currently sitting on a ticking time bomb that will result in an
unprecedented explosion in the quantity and quality of data. The Dark Energy
Spectroscopic Instrument (DESI) (DESI Collaboration et al., 2016) has already
recorded more galaxy spectra in its first two years of operations than
previously done in the whole history of humanity. Moreover, next-generation
surveys such as Legacy Survey of Space and Time (LSST) (Ivezić et al., 2019),
Euclid (Laureijs et al., 2011), the Nancy Grace Roman space telescope (Spergel
et al., 2015), the Simons Observatory (Ade et al., 2019), CMB-HD (Sehgal et
al., 2019) or CMB-S4 (Abazajian et al., 2016) promise to further accelerate
this trend in the next decade.

In order to match the quality of the data, physicists are starting to
incorporate into their theoretical predictions more of the physical,
observational and instrumental effects which, until now, could be overlooked.
In practice this translates into a dramatic increase in the number of
parameters that future analyses will have to consider. This combination of
large data sets with complex models will (and in many cases already does)
overwhelm the inference methods we currently use to constrain the values of
these parameters.

In a Bayesian framework, the statistical distribution of a set of parameters,
$\boldsymbol{\theta}$, given some data, $\boldsymbol{d}$, is given by Bayes
theorem:

$\mathcal{P}(\boldsymbol{\theta}|\boldsymbol{d})\propto\mathcal{L}(\boldsymbol{d}|\boldsymbol{\theta})\Pi(\boldsymbol{\theta})\,,$
(1)

where $P(\boldsymbol{\theta}|\boldsymbol{d})$ is the posterior distribution of
the parameters given the data,
$\mathcal{L}(\boldsymbol{d}|\boldsymbol{\theta})$ is the likelihood of the
data for a set of parameters and $\Pi(\boldsymbol{\theta})$ are the prior
beliefs on the distribution of the parameters. Current cosmological analyses
explore their parameter spaces by mapping out
$P(\boldsymbol{\theta}|\boldsymbol{d})$ according to some stochastic process
by which the direction of exploration is chosen. Despite the success of this
methodology, relying on random methods becomes inefficient at high dimensions
as the chances of randomly finding the parameters that describe the observed
data well decreases with the volume of the parameter space being probed. This
effect is known as the ”curse of dimensionality”.

One of the most effective ways of overcoming the curse of dimensionality is
using the gradient of the posterior;
$\nabla{P(\boldsymbol{\theta}|\boldsymbol{d})}\propto\nabla{(\mathcal{L}(\boldsymbol{d}|\boldsymbol{\theta})\,\Pi(\boldsymbol{\theta}))}$,
to guide exploration towards regions of interest in parameter space.
Algorithms that use the gradient of the likelihood to explore parameter space
are known as gradient-based samplers. Unfortunately, traditional methods of
finding numerical gradients, such as the finite-difference method, scale
poorly with the number of dimensions. Moreover, they can be prone to numerical
instabilities. This can render the use of gradient information counter-
productive.

Thankfully, in the last decades a series of algorithms known as auto-
differentiation (AD) (Griewank & Walther, 2008; Margossian, 2018; Harrison,
2021) have grown in popularity. Given a generic computer program that maps a
series of inputs to a series of outputs (i.e. a function inside a computer) AD
is a family of algorithms designed to systematically draw another computer
program that returns the gradient of the former. Unlike finite differences, AD
is subject to truncation errors and its computational cost scales much more
favourably with the dimensionality of the original function (Griewank &
Walther, 2008). Thus, the goal of this paper is to provide an AD cosmological
inference framework.

We present LimberJack.jl, an auto-differentiable angular power spectra
analysis code fully written in Julia. LimberJack.jl is designed after the Core
Cosmology Library (CCL Chisari et al., 2019) developed by the LSST Dark Energy
Science Collaboration, aiming to fulfil the similar scientific goals.
LimberJack.jl allows the user to easily compute $\Lambda$CDM model predictions
for the angular power spectra of weak lensing, CMB lensing and clustering
surveys using the Limber approximation (Limber, 1953; LoVerde & Afshordi,
2008). Most importantly, LimberJack.jl can also provide the user with accurate
gradients of these predictions in a computationally efficient way, due to its
compatibility with Julia’s AD libraries ForwardDiff.jl and ReverseDiff.jl.
While LimberJack.jl is currently more limited than CCL, its modular design
means that it can easily be extended by the community.

A number of auto-differentiable cosmological codes have been presented in the
literature. The most notable precedents are JAX-COSMO (Campagne et al., 2023)
and CosmoPower-JAX (Piras & Spurio Mancini, 2023). JAX-COSMO’s functionalities
and goals are very similar to those of LimberJack.jl both basing their design
on CCL. CosmoPower-JAX is a neural-network emulating framework for the matter
power spectrum originally developed in TensorFlow (Spurio Mancini et al.,
2022) and later ported to JAX to be paired with JAX-COSMO. The main difference
between JAX-COSMO, CosmoPower-JAX and LimberJack.jl is the AD ecosystem they
make use of. Both JAX-COSMO and CosmoPower-JAX are written in JAX, a scripting
programming language which interfaces with a lower-level language, XLA, a
just-in-time (JIT) compiled language with AD capabilities. JAX’s main
strengths are its powerful parallelisation schemes on both CPU’s and GPU’s,
its performant AD methods and the fact that it shares API with the ubiquitous
Python library NumPy. On the other hand, LimberJack.jl is fully written in
Julia, a general-purpose, JIT-compiled, programming language with native AD
capabilities. Julia’s main advantage is the lack of a lower-level programming
language as is in the case of other popular AD environments such as TensorFlow
interfacing with C++ and JAX with XLA. This transparency makes Julia an
excellent language to develop complex libraries with customised methods. A
recent practical example where this feature has played a key role is the
development of the first AD Boltzmann solver, Bolt.jl Li & Sullivan (2023), a
challenging task for JAX and TensorFlow but which, nonetheless, was possible
in Julia.

The structure of this paper is as follows: in Sect. 2 we describe the
theoretical predictions computed by LimberJack.jl, how LimberJack.jl computes
these predictions in an auto-differentiable way and how AD can be used to
speed up statistical inference. In Sect. 3, we use LimberJack.jl to reproduce
the DES Y1 3x2 analysis and to perform a Gaussian process reconstruction of
the growth factor across redshift. Finally, in Sect. 4, we summarise our work
and interpret our results. In addition to this we also present a tutorial on
how to use LimberJack.jl in App. A.

##  2 Methods

###  2.1 Angular Power Spectra

The main goal of LimberJack.jl is to compute the $\Lambda$CDM model
predictions for the angular power spectra of any two tracers using the Limber
approximation as well as their gradients. In the Limber approximation (Limber,
1953; LoVerde & Afshordi, 2008), the angular power spectrum, $C^{UV}_{\ell}$,
between two projected fields $U$ and $V$ is related to the three-dimensional
power spectrum, $P_{UV}(k,z)$ of their associated quantities $U$ and $V$ via:

$C^{UV}_{\ell}=\int\frac{d\chi}{\chi^{2}}q_{U}(\chi)q_{V}(\chi)\,P_{UV}\left(k=\frac{\ell+1/2}{\chi},z(\chi)\right)\,,$
(2)

where $q_{U}$ and $q_{V}$ are the radial kernels of the $U$ and $V$ fields
respectively.

Currently LimberJack can compute the auto- and cross-correlation of three
different fields: cosmic shear (i.e. the distortion in the shape of galaxies
caused by weak gravitational lensing), CMB lensing convergence (i.e. the
apparent magnification or reduction of CMB anisotropies due to gravitational
lensing by the matter distribution along the line of sight) and galaxy
clustering.

The radial kernels cosmic shear ($q_{\gamma}$), CMB lensing convergence
($q_{\kappa}$) and galaxy clustering ($q_{g}$) are given by:

$\displaystyle
q_{\gamma}(\chi)\equiv\frac{3}{2}H_{0}^{2}\Omega_{\rm{m}}\frac{\chi}{a(\chi)}\int_{z(\chi)}^{\infty}dz^{\prime}p(z^{\prime})\frac{\chi(z^{\prime})-\chi}{\chi(z^{\prime})},$
(3) $\displaystyle
q_{\kappa}(\chi)\equiv\frac{3}{2}H_{0}^{2}\Omega_{\rm{m}}\frac{\chi}{a(\chi)}\frac{\chi^{*}-\chi}{\chi^{*}},$
(4) $\displaystyle q_{g}(\chi)\equiv\frac{H(z)}{c}p(z)$ (5)

where $c$ is the speed of light, $H(z)$ is the Hubble expansion rate,
$H_{0}\equiv H(z=0)$, $\Omega_{\rm{m}}$ is the matter density parameter today,
$b_{g}$ is the galaxy-matter linear bias parameter,$\chi^{*}$ is the comoving
distance to the surface of last scattering and $p(z)$ is the redshift
distribution of the galaxies.

The estimation of the galaxy clustering, shear and convergence fields is
subject to a variety of systematic biases that must be accounted for in their
modelling. In galaxy clustering analyses the level of detail to which the
relation between the observed galaxy overdensity and the predicted matter
overdensity is modelled limits the range of scales that can be analysed. We
consider a single parameter linear bias model, $b_{\rm{g}}$, that multiplies
the clustering kernel. This restricts our analysis of galaxy clustering data
to relatively large scales. Cosmic shear analysis must account for possible
biases in galaxy shape estimation. Similarly, we consider a single
multiplicative bias parameter, $m$. In addition, galaxy shapes are not only
correlated by cosmic shear but also by intrinsic alignments (IAs) in the
orientation of galaxies due to local interactions. Within the so-called Non-
Linear Alignment (NLA, Hirata & Seljak, 2004) model this can be accounted for
by adding an extra contribution to the final shear kernel given by:

$q_{I}(\chi)=A_{\rm{IA}}(z)\frac{H(z)}{c}p_{\alpha}(z),$ (7)

where $A_{\rm{IA}}(z)$ is:

$A_{\rm IA}(z)=A_{{\rm
IA},0}\left(\frac{1+z}{1+z_{0}}\right)^{\eta_{\text{IA}}}\frac{0.0139\Omega_{\rm{m}}}{D(z)}\,,$
(8)

where $A_{\text{IA},0}$ and $\eta_{\text{IA}}$ are two free parameters,
$z_{0}$ is a redshift pivot (which we fix to $z_{0}=0.62$ as in Abbott et al.
(2018a); Troxel et al. (2018)), and $D(z)$ the linear growth factor.

When the shear and convergence kernels enter Eq. 2 they must be multiplied by
an $\ell$-dependent correction to account for the differences between angular
and three-dimensional derivatives of the matter density field. These
corrections are given by:

$\displaystyle
W_{\rm{\ell}}\equiv\sqrt{\frac{(\ell+2)!}{(\ell-2)!}}\frac{1}{(\ell+1/2)^{2}}\,$
(9) $\displaystyle K_{\rm{\ell}}\equiv\frac{\ell(\ell+1)}{(\ell+1/2)^{2}}\,$
(10)

for the shear and convergence kernels respectively.

Since all the fields considered are ultimately projections of the
3-dimensional matter density field we only need to consider the matter auto-
correlation power spectrum, $P(k,z)\equiv P_{mm}(k,z)$, when computing Eq. 2.
Assuming that the evolution of the matter power spectrum is scale independent,
$P(k,z)$ can be obtained by first computing the value of the linear matter
power spectrum at $z=0$ at an array of scales. As we will discuss in Sect. 2.2
this can be done either using fitting formulas or emulating the solutions of
Botlzmann codes. Then, the linear matter power spectrum can be evolved
backwards in time as $P^{\rm{L}}(k,z)=D^{2}(z)P_{0}(k)$ where $D(z)$ is the
linear growth factor found by solving the Jeans equation:

$\frac{d}{da}\left(a^{2}H\frac{dD}{d\log
a}\right)=\frac{3}{2}\Omega_{\rm{m}}(a)\,a\,H\,D\,.$ (11)

However, in order to fit the small scales it is necessary to compute the non-
linear evolution of the matter power spectrum, $P^{\rm{NL}}(k,z)$. We do so by
applying the Halofit fitting function Smith et al. (2003) with revisions from
Takahashi et al. (2012) to $P^{\rm{L}}(k,z)$. However, analyses of data from
the latest surveys (Aricò et al., 2023; Energy Survey et al., 2023) require
the use of the more accurate methods such as HMCode (Mead et al., 2021) or
baccoemu (Angulo et al., 2021) as smaller scales are included.

###  2.2 LimberJack.jl

The key feature of LimberJack.jl which sets it apart from similar, more
extensively tested codes (such as CCL Chisari et al. (2019)) are its AD
capabilities. AD methods can be classified into two groups, forwards and
backwards. In order to understand the difference, let us consider a
complicated function of an independent variable, $f(x)$, one can represent the
function as a composition of simpler functions,
$f(x)=w_{n}(w_{n-1}(...w_{1}(w_{0}))$ where $w_{0}=x$. Thus in order to obtain
the gradient of $f(x)$, one can either find the gradient of the simpler
functions with respect the independent variable or, alternatively, find the
gradient of the original function with respect the simpler functions. The
first of these different strategies to accumulate terms in the chain rule
corresponds to forward AD while the second to backwards AD. Whether to use
forwards or backwards AD will greatly depend on the nature of the problem.
Generally111The performance of backwards and forwards AD is also impacted by
the size of the operation being differentiated., given a function
$f\colon\mathbb{R^{N}}\to\mathbb{R^{M}}$, forwards AD will be more efficient
at computing $\nabla f$ if $N<M$. If, on the other hand, $M<N$, backwards AD
is preferred.

LimberJack.jl is capable of both forwards and backwards AD through the Julia
libraries ForwardDiff.jl and ReverseDiff.jl respectively. Both of these
libraries achieve AD by pushing special types of numbers through the original
computer program. On the one hand, ForwardDiff.jl uses dual numbers. Dual
numbers are expressions of the form $a+\epsilon b$ such that $\epsilon^{2}=0$
but $\epsilon\neq 0$. It can then be shown that given a generic function
$f(a+b\epsilon)=f(a)+b\epsilon f^{\prime}(a)$ (see Sect. 4 of Hoffmann (2014)
for proof) such that the gradient of the original function can be easily
obtained. This means that ForwardDiff.jl imposes little to no constraints on
the program it differentiates through222Since Julia is a typed language, some
considerations have to be made about potential type instabilities introduced
by the use of dual numbers.. On the other hand, ReverseDiff.jl uses taping
numbers, a special type of numbers that records all the operations the number
undergoes, to generate a trace of the basic operations that compose a generic
computer program (van Merriënboer et al., 2018). Given this record,
ReverseDiff.jl can then generate an expression for the gradient of the
program. Thus ReverseDiff.jl requires two passes through the original program.
First, a forward pass generates the program’s operation trace. Then, a
backwards pass computes the partial derivatives and accumulates their values
as the input is back-propagated. Because of the greater complexity of the
algorithm, ReverseDiff.jl imposes strong demands on the computer programs it
acts upon. Commonplace computations such as control flow or variable mutation
are examples of operations that must be handled carefully.

The implementation of the expressions discussed in Sect. 2.1 within
LimberJack.jl will are thus constrained by the demands of the AD methods used
to obtain their gradients. The computation of Eq. 2 fundamentally involves two
different types of quantities, the radial kernels and the matter power
spectrum. On the one hand computing the radial kernels boils down to
evaluating the expansion history of the Universe, $H(z)$, and its integral
over time, i.e. the comoving radial distance, $\chi(z)$. Calculating the
expansion history, as given by the $\Lambda$CDM model, amounts to evaluating
$H(z)=H_{0}\sqrt{\Omega_{\rm{m}}(1+z)^{3}+\Omega_{r}(1+z)^{4}+(1-\Omega_{\rm{m}}-\Omega_{r})}$
on a grid of redshifts where $\Omega_{\rm{m}}$ and $\Omega_{r}$ are the
cosmological matter and radiation density respectively. The radial comoving
distance is then obtained integrating over the grid using Simpson numerical
integration which is generically compatible with AD methods. In the left and
centre columns of Fig. 1 we show a comparison between the LimberJack.jl and
CCL predictions for the expansion history and comoving distance between
$0<z<1100$ (top panels). Each panel contains a subpanel where the relative
difference between the CCL and LimberJack.jl predictions is shown. We find
that the relative difference between the predictions of CCL and LimberJack.jl
are smaller than $10^{-}4$ at all redshift. Moreover we also show a comparison
for the derivative of said quantities with respect to $\Omega_{\rm{m}}$ when
computed by LimberJack.jl using AD and finite differences for the same
redshift range (bottom panels). Similarly, each panel contains a subpanel
where the relative difference between the AD and numerical gradient is shown.
We find that the relative difference between the two methods is smaller than
$10^{-}4$ at all redshift.

Figure 1: Left column: Top panel shows a comparison between the LimberJack.jl
(solid green) and the CCL (dashed red) computation of the expansion history
between $0<z<1100$. Bottom panel shows a comparison between the derivative of
the expansion history with respect to $\Omega_{\rm{m}}$ computed using AD
(solid blue) and finite differences (dashed orange) for the same redshift
range. Middle column: Top panel shows a comparison between the LimberJack.jl
(solid green) and the CCL (dashed red) computation of the comoving distance
between $0<z<1100$. Bottom panel shows a comparison between the derivative of
the comoving distance with respect to $\Omega_{\rm{m}}$ computed using AD
(solid blue) and finite differences (dashed orange) for the same redshift
range. Right column: Top panel shows a comparison between the LimberJack.jl
(solid green) and the CCL (dashed red) computation of the linear growth factor
$0<z<3$. Bottom panel shows a comparison between the derivative of the linear
growth factor with respect to $\Omega_{\rm{m}}$ computed using AD (solid blue)
and finite differences (dashed orange) for the same redshift range.

On the other hand computing the matter power spectrum in an AD-compatible way
is a much more difficult problem. The first and most challenging obstacle is
obtaining the linear matter power spectrum. As described in Sect. 2.1,
computing the matter power spectrum involves solving a coupled system of
linear differential equations with time- and scale-dependent coefficients.
This task has proven a major challenge for most AD environments (Campagne et
al., 2023) as it requires tweaking the lower-level programming languages. This
is however not a problem for the Julia AD ecosystem which has a (currently
under development) full Boltzmann-Einstein solver, Bolt.jl Li & Sullivan
(2023). While still in its early days, Bolt.jl can provide users with full
numerical solutions for the linear power spectrum and its derivatives with
respect to the $\Lambda$CDM parameters. LimberJack.jl can then interface with
Bolt.jl to obtain said predictions.

However, solving the full Boltzmann-Einstein equations is very computationally
expensive (being the bottleneck of most cosmological analyses) even in fast
programming languages such as C++ (CLASS), Fortran (CAMB) or Julia. Therefore,
it is common to look for ways to bypass solving the Boltzmann-Einstein
equations in cosmological analyses where speed is important. When the
evolution of the matter power spectrum can be assumed to be scale independent,
$P(z,k)$ can be computed by constructing a fitting formula that approximates
the true value of the matter power spectrum at $z=0$, $P_{0}(k)$, at an array
of scales and then evolving $P_{0}(k)$ into the past using the linear growth
factor. The first commonly applied fitting formula for $P_{0}(k)$ was the
Bardeen-Bond-Kaiser-Szalay (BBKS Bardeen et al., 1986) formula which modifies
the Harrison-Zeldovich power spectrum (Dodelson, 2003) by a transfer function:

$P_{0}(k)=\frac{8\pi^{2}}{25}\frac{A_{s}}{\Omega_{rm{m}}^{2}}T^{2}(k)\frac{k^{n_{s}}}{H_{0}^{4}k_{p}^{n_{s}-1}}\,,$
(12)

where is $A_{s}$ is the amplitude of the matter power spectrum,
$k_{p}=0.05Mpc^{-1}$ and $T(k)$ is the transfer function for a given set of
fixed cosmological parameters. The most popular fitting formula is that of
(Eisenstein & Hu, 1999, 1998, E&H from now on), which follows the same
strategy as the BBKS formula of approximating the transfer function but with a
more complex expression that includes baryonic effects such as baryonic
acoustic oscillations and small-scale power suppression. Thanks to these
inclusions the E&H formula is accurate enough to return unbiased cosmological
constraints for the 3x2 DES Y1 analysis as shown in Campagne et al. (2023).
However, it is important to keep in mind that the E&H formula will not be
accurate enough to analyse future data set as well as some current ones.

In recent years, a new family of fitting formulae, known as emulators, have
grown in popularity. Emulators are computer models (such as neural networks or
Gaussian processes) whose weights are optimized to reproduce a target function
over a certain domain. The main advantage of emulators over traditional
fitting formulae is that they automate the majority of the trial and error
process of building an accurate fitting formula. Moreover, they tend to
require less knowledge of the physical problem to be constructed. However,
emulators require vasts amounts of training data and the resulting models tend
to be far larger than traditional fitting formulae. Nonetheless, emulators
have found great success in astrophysics in recent years, offering extremely
accurate yet affordable approximations to different functions (Donald-McCann
et al., 2022; Bonici et al., 2023), including the matter power spectrum
(Spurio Mancini et al., 2022; Mootoovaloo et al., 2022).

The great advantage of these fitting formulae is that they effectively amount
to algebraic expressions through which AD algorithms can easily differentiate
as shown in Piras & Spurio Mancini (2023) and Bonici et al. (2023). Combined
with their speed, fitting formulae can often be the preferred way to obtain
estimates for the primordial power spectrum in gradient-based analyses. For
these reasons, LimberJack.jl is equipped with both a native Julia
implementation of the E&H formula and the recently developed Gaussian process
based emulator EmuPk (Mootoovaloo et al., 2022) for the linear matter power
spectrum.

The next challenge is to compute the growth factor, $D(z)$, to find the
evolution of the linear power spectrum. As discussed in Sect. 2.1, the growth
factor is obtained by solving the Jeans equation (see Eq. 11) which is an
inhomogeneous ordinary differential equation. Differentiating through the
solutions of differential equations can be done by writing a recursive
numerical scheme to solve the differential equation. These schemes amount to a
series of linear operations that update mutating variables. While mutating
variables can pose challenges for certain backwards modes of AD, the numerical
schemes can otherwise be differentiated through to yield gradients for the
solution of a differential equation. Thus, LimberJack.jl solves the Jeans
equation using a second order Runge-Kutta solver. Returning to Fig. 1, in
right column of this figure we show a comparison between the LimberJack.jl and
CCL predictions for the linear growth factor (top panels). We also shown a
comparison between the derivatives of LimberJack.jl’s $D(z)$ with respect to
$\Omega_{\rm{m}}$ when computed using AD and finite difference (bottom panels)
for the redshift range $0<z<3$. Again, each panel contains a subpanel where
the relative difference between the two compared quantities is shown.
Concerning the growth factor itself, we find that the relative difference
between the CCL and LimberJack.jl predictions are smaller than $10^{-}4$ for
all the redshift window. Concerning the derivatives of the growth factor, the
relative differences between AD and finite differences are again smaller
$10^{-}4$ for all the redshift window.

To compute the necessary non-linear corrections at small scales, LimberJack.jl
is equipped with a Juia Halofit implementation fully compatible with AD. The
biggest obstacle to accomplish this is differentiating through the root
finding process that occurs within Halofit. Similarly to solving differential
equations, root finding can be differentiated through by writing a recursive
numerical scheme of linear operations with mutable variables (Margossian &
Betancourt, 2021). The Halofit implementation within LimberJack.jl uses the
secant method for its root finding through which most AD methods can
differentiate. Putting all of these together LimberJack.jl can perform a fully
auto-differentiable computation of the non-linear matter power spectrum at any
redshfit or scale. In Fig. 2 we show a comparison between the predictions for
the non-linear matter power spectrum from LimberJack.jl and CCL for
$z=[0.0,0.5,1.0,2.0]$ and $-2<\log_{10}(k)<4$ (top panels). Both the CCL and
LimberJack.jl predictions are obtained using the E&H linear power spectrum and
Halofit prescription to obtain the non-linear corrections. Moreover we also
show a comparison between the derivative of the LimberJack.jl power spectra
with respect to $\Omega_{\rm{m}}$ when computed using AD and finite
differences. Finally, each panel contains a subpanel where the relative
difference between the two compared quantities is shown. Concerning the power
spectra, we find that the relative difference between the CCL and
LimberJack.jl predictions are smaller than $1*10^{-}3$ for all the redshift
window. However, we observe that the Halofit implementation of LimberJack.jl
systematically over-predicts the non-linear matter power spectrum at higher
redshits. As we will see this bias will manifest when considering angular
power spectra involving the CMB lensing tracer later on but it nonetheless
doesn’t prevent us from achieving our accuracy goals. Concerning the
derivatives of the same power spectra, the relative differences between AD and
finite differences are again smaller $1*10^{-}3$ for all the redshift window.

Figure 2: Top panels show a comparison between the LimberJack.jl (solid green)
and the CCL (dashed red) computation of the non-linear matter power spectrum
using Halofit for $z=[0.0,0.5,1.0,2.0]$ and $-4<log_{1}0(k)<2$. Bottom panels
shows a comparison between the derivative of the LimberJack.jl matter power
spectra with respect to $\Omega_{\rm{m}}$ computed using AD (solid blue) and
finite differences (dashed orange). Each panel contains a subpanel where the
relative difference between the two compared quantities is shown.

The only step left to compute angular power spectra is to bring the radial
kernels and the matter power spectrum together and perform the Limber integral
(see Eq. 2). Similarly to the computation of comoving distances, this is done
within LimberJack.jl by evaluating all the quantities within the integrand at
a regular array of logarithmic scales, $log_{10}(k)$, and then performing the
numerical integral for the desired multi-poles, $\ell$, using Simpson
numerical integration. One small challenge in doing so is finding the
corresponding redshift associated with the comoving distance given by the
scale and the multipole to evaluate the radial kernels. This inconvenience is
a result of LimberJack.jl defining the radial kernels as functions of redshift
instead of comoving distance. Normally, finding the redshift at a given
comoving distance would involve a costly root finding process. However,
LimberJack.jl handles this by building an interpolator between redshift and
comoving distance and then inverting the interpolator to establish a straight
forward comoving distance to redshift map for a given set of cosmological
values.

In Fig. 3 we show a comparison between the predictions of LimberJack.jl and
CCL for the angular power spectra of different types of tracers for
$10<\ell<1000$. We consider the auto- and cross-correlations of galaxy
clustering, cosmic shear and CMB convergence. In all cases the E&H formula was
paired with Halofit to obtain the matter power spectrum. We also assume a
Gaussian redshift distribution centred at $z=0.5$ with a standard deviation
$\sigma_{z}=0.05$, sampled at 1000 evenly-spaced intervals in the range
$0<z<2$. We observe that the discrepancy between LimberJack.jl and CCL is
smaller than $10^{-}3$ for all angular power spectra except for the auto-
correlation of CMB lensing where the discrpancy is larger but nontheless stays
below $10^{-}2$. This is due to the discrepancy in the evolution of the
Halofit non-linear matter power spectrum between LimberJack.jl and CCL shown
in Fig. 2. Similarly, in Fig. 4 we show a comparison for the derivatives of
the same quantities with respect to $\Omega_{\rm{m}}$ when computed using AD
and finite differences. Despite the clustering tracer introducing highly
oscillatory behaviour in the derivatives of the angular power spectra we
observe a sub-percentage-level agreement between the two approaches in all
cases.

Figure 3: Shows a comparison between the LimberJack.jl (solid green) and the
CCL (dashed red) computation of auto- and cross-correlation angular power
spectra of galaxy clustering, weak lensing and CMB lensing tracers for the
range of multipoles $1<log(\ell)<3$ Figure 4: Shows a comparison between the
derivatives of the auto- and cross-correlation angular power spectra of galaxy
clustering, weak lensing and CMB lensing tracers with respect to
$\Omega_{\rm{m}}$ computed using AD (solid blue) and finite differences
(dashed orange) for the range of multipoles $1<log(\ell)<3$.

###  2.3 Gradient-based Samplers

In the previous sections we described how the computation of angular power
spectra (Sect. 2.1) is made compatible with AD methods within LimberJack.jl
(Sect. 2.2). In this section we will describe how these newly found gradients
can be used to significantly speed up the statistical inference of
cosmological parameters, the goal of this paper.

As described in Sect. 1 stochastic inference methods become prohibitively
expensive as the number of parameters (dimensions) increases. Gradient-based
samplers provide a solution to this problem by using the gradient of the
posterior to guide the exploration of the parameter space towards the regions
where the probability density of the likelihood is the highest.

Hamiltonian Monte Carlo (HMC, MacKay, 2002) is a sampling algorithm that
simulates Hamiltonian trajectories on the parameter space to guide the
exploration. In order to do so HMC duplicates the parameter space by
introducing a set of auxiliary momentum variables ($\boldsymbol{p}$)
independent of the likelihood, and sampled from a d-dimensional unit variance
Gaussian distribution. The joint distribution of the position (i.e the
original parameters of interest $\boldsymbol{\theta}$) and momentum variables
is then defined by a Hamiltonian function which governs the dynamics of the
system:

$H(\boldsymbol{\theta},\boldsymbol{p})=T(\boldsymbol{p})+V(\boldsymbol{\theta})=\frac{\boldsymbol{p}^{2}}{2M}-\log(\mathcal{L}(\boldsymbol{\theta}))\,,$
(13)

where $M$ is the mass matrix, a positive-definite matrix that controls the
mixing of momenta variables.

In each iteration of HMC, the Hamiltonian dynamics are simulated using a
numerical integration scheme to generate the next sample. This is done using a
leapfrog scheme given its symplectic properties that preserve the volume of
the phase space along the trajectory. The sample is then accepted or rejected
using a Metropolis adjustment based on its Hamiltonian energy. This ensures
that the target posterior distribution is obtained upon marginalizing over the
momenta variables, i.e:

$\mathcal{P}(\boldsymbol{\theta})\propto\int\exp(-H(\boldsymbol{\theta},\boldsymbol{p}))\,d\boldsymbol{p}\,.$
(14)

While HMC can dramatically speed up the exploration of large parameter spaces,
it also introduces a series of challenges. Namely, the user must find the
ideal step size for the numerical solver of the Hamiltonian trajectories, the
ideal length of the trajectories and an expression for the mass matrix. The
No-U-turns sampler (NUTS) (Hoffman & Gelman, 2011) is an HMC algorithm that
dynamically tunes these quantities based only on two user provided quantities:
the Target Acceptance Probability (TAP) of the Metropolis adjustment and the
number of adaptation steps. The fundamental idea behind NUTS is to prevent the
sampler from returning to previously explored regions by avoiding U-turns in
phase space. This is done by evolving the Hamiltonian trajectories both
backwards and forwards in time (known as branching) and stopping when the
momenta start to turn on themselves on either branch. This establishes a
principled criterion for the length of the Hamiltonian trajectories that
maximises the statistical independence of the samples. Moreover, during the
adaptation phase the step size of the integrator is tuned such that the TAP is
achieved. Finally, NUTS tunes the mass matrix using the variance of the
adaptation phase samples such that the covariance matrix of the momenta
variables is as close as possible to unit-variance.

##  3 Results

###  3.1 DES Y1 3x2

Figure 5: Shows the marginalised posteriors distributions from the analysis of
DES Y1 3x2 data for the cosmological parameters $\Omega_{\rm{m}}$,
$\sigma_{\rm{8}}$ and $S_{\rm{8}}$. Lower panel compares the LimberJack.jl
(solid blue) and CCL (dashed red) analyses when using the E&H formula. Upper
triangle compares the LimberJack.jl (solid purple) and CCL (dashed green)
analyses when using EmuPk and CLASS respectively. Cobaya contours were
obtained using the MH sampler while LimberJack.jl contours were obtained using
the NUTS sampler.

In order to validate the constraints obtained using LimberJack.jl we first
replicate two different DES Y1 3x2 (Abbott et al., 2018b) analyses done using
the well-established Bayesian inference framework Cobaya which employs a
Metropolis Hastings (MH) sampler. In the first analysis, both LimberJack.jl
and Cobaya used the E&H formula to obtain the matter power spectrum. In the
second analysis, Cobaya used the CLASS code to obtain the matter power
spectrum while LimberJack.jl used the EmuPk emulator Mootoovaloo et al.
(2022), which was trained on CLASS.

With this aim, we make use of the Julia library for statistical inference
Turing.jl (Ge et al., 2018). Turing.jl is a probabilistic programming language
(PPL) that allows the user to create statistical models by explicitly writing
the relationship between the sampled parameters and the theoretical
predictions for the observed data. Turing.jl uses this information to draw a
likelihood density function compatible with the Julia AD infrastructure as
long as the internal computations of the model are also compatible with AD.
The user can then condition this function on the observations and sample it
using a series of sampling back-ends. In this work we make use of the NUTS
sampler as implemented in the Julia library AbstractHMC.jl within Turing.
Moreover, we choose the AD library ForwardDiff.jl to provide NUTS with the
gradient of the likelihood. Note that using forwards AD for statistical
inference is sub-optimal since the likelihood is a function that maps a high-
dimensional space to a single scalar. However, we found better performance
than when using the backwards-AD library ReverseDiff.jl in our particular
application.

DES Y1 3x2 Analyses Priors
---
Parameter | Prior | Parameter | Prior
Cosmology | Shear calibration bias
$\Omega_{\rm{m}}$ | $U(0.2,0.6)$ | $m^{i}$ | $N(0.012,0.023)$
$\Omega_{\rm{b}}$ | $U(0.028,0.065)$ | Galaxy bias
$h$ | $U(N(70,5),0.55,0.91)$ | $b_{g}^{i}$ | $N(0.8,3.0)$
$n_{\rm{s}}$ | $U(0.87,1.07)$ | Intrinsic Alignments
$\sigma_{\rm{8}}$ | $U(0.6,0.9)$ | $A_{\mathrm{IA,0}}$ | $U(-5.0,5.0)$
|  | $\eta_{\mathrm{IA}}$ | $U(-5.0,5.0)$
Lens redshift calibration | Sources redshift calibration
$\Delta z_{1}$ | $N(0.0,0.007)$ | $\Delta z_{1}$ | $N(-0.001,0.016)$
$\Delta z_{2}$ | $N(0.0,0.007)$ | $\Delta z_{2}$ | $N(-0.019,0.013)$
$\Delta z_{3}$ | $N(0.0,0.006)$ | $\Delta z_{3}$ | $N(0.009,0.011)$
$\Delta z_{4}$ | $N(0.0,0.01)$ | $\Delta z_{4}$ | $N(-0.018,0.022)$
$\Delta z_{5}$ | $N(0.0,0.01)$ |  |
Table 1: Prior distributions for the DES Y1 3x2 analysis parameters based on
the prior range of EmuPk (Mootoovaloo et al., 2022)

The same priors were used in the two cases for both the Cobaya and LimberJack
analyses. A summary of the priors can be found in Tab. 1. We show the
resulting posteriors in Fig. 5. We observe an excellent agreement between the
LimberJack.jl and Cobaya pipelines regardless of whether the E&H formula or
CLASS/EmuPk is used to obtain the primordial matter power spectrum.

We compare the performance of the two samplers by looking at the number of
effective samples (i.e. number of statistically independent samples in the
Markov chain Gelman et al., 2013) per number of likelihood calls. This metric
is independent of the hardware used to run the analysis and the time taken by
LimberJack.jl and Cobaya to evaluate the likelihood. Thus this metric allows
us to look at the improvement purely brought about by using gradient-based
inference methods. Analysing the Markov chains of the different samplers we
find that NUTS is approximately 1.5 times more efficient than Cobaya at
generating effective samples. However, in order to fairly compare the
performance of the two samplers we must take into account that NUTS computes
the gradient of the likelihood at every step. Given the 25 parameters of the
DES Y1 3x2 analysis, evaluating the gradient of the likelihood using
ForwardDiff.jl is ten times more expensive than evaluating the likelihood
itself. In order to compensate for this extra cost the efficiency improvement
of NUTS over MH should be equal or larger than the relative cost of computing
the gradient of the likelihood. Therefore the efficiency improvement of NUTS
is not enough to compensate for the cost of computing the gradient in this
particular application. There exist two avenues to tilt the balance in favour
of NUTS. On the one hand, increasing the efficiency of the sampler. On the
other hand, reducing the cost of the likelihood gradient. However, the current
bottleneck is the efficiency of the sampler as even more specialised AD
methods (such as Zygote.jl) would take at least twice the time to compute the
likelihood to obtain its gradient. Future works could explore initialising the
NUTS mass matrix using the estimate of the variational inference code
PathFinder.jl (Zhang et al., 2021). Therefore, in order to showcase the
strength of the methods developed on this work we need an application for
which the efficiency of gradient-based samplers truly outpaces that of
traditional inference methods.

###  3.2 Growth Factor Reconstruction

In the previous section we showed the reliability of LimberJack.jl by
reproducing the official DES Y1 3x2 analysis. In this section we will showcase
how LimberJack.jl can be used to perform statistical inference on models
outside of the reach of traditional inference methods. In order to do so we
perform a model-independent reconstruction of the growth factor using a
Gaussian process (GP) with more than a hundred parameters.

A GP is a collection of random variables (nodes), each of them sampled from a
multivariate Gaussian distribution with a non-diagonal covariance (Rasmussen &
Williams, 2006; Seikel et al., 2012). Thus a GP – $g(\textbf{x})$ where x is a
arbitrary vector representing the position of the nodes – is fully specified
by a mean function, $m(\textbf{x})\equiv\mathbb{E}[g(\textbf{x})]$ – where
$\mathbb{E}[\cdots]$ is the expectation value over the ensemble – and a
covariance function,
$k(\textbf{x},\textbf{x'})\equiv\mathbb{E}[(g(\textbf{x})-m(\textbf{x}))(g(\textbf{x'})-m(\textbf{x'}))]$.
In combination, the mean and covariance functions determine the statistical
properties of the random variables that define the family of shapes the GP can
take.

In this work we use a GP composed of 101 nodes equally spaced through the
redshift window $0\leq z\leq 3$. The mean of the GP is given by the best-fit
$\Lambda$CDM Planck 2018 (Planck Collaboration et al., 2020, P18 from now on)
prediction for $D(z)$. While this choice of mean will bias the reconstruction
towards the P18 prediction outside of the data domain, its impact in the data
dominated regions is small (see e.g. Ruiz-Zapatero et al., 2022b; Ruiz-
Zapatero et al., 2022a, for similar examples of this GP behaviour in the
literature.). Regarding the covariance matrix, we choose a square exponential
covariance function, defined as:

$\displaystyle
K(\boldsymbol{z},\boldsymbol{z}^{\prime})=\eta^{2}\exp\frac{-\mathfrak{D}(\boldsymbol{z},\boldsymbol{z}^{\prime})}{2l^{2}}\,,$
(15)

where $\eta$ is the amplitude of the oscillations around the mean and $l$ is
the correlation length between the GP realizations and
$\mathfrak{D}(\boldsymbol{z},\boldsymbol{z}^{\prime})$ is the distance matrix
of $\boldsymbol{z}$ and $\boldsymbol{z}^{\prime}$ given by
$\mathfrak{D}_{ij}=|z_{i}-z^{\prime}_{j}|^{2}$. This decision was made based
on the fact that the square exponential is a computationally inexpensive and
infinitely differentiable kernel, appropriate to model smooth fluctuations
around the mean of the GP.

Therefore, the growth factor is given by:

$D(\boldsymbol{z})=D_{\rm{P18}}(\boldsymbol{z})+L_{C}(K(\boldsymbol{z},\boldsymbol{z}))\,\boldsymbol{v}\,,$
(16)

where $\boldsymbol{v}$ is a vector of the GP nodes sampled from a unit
variance, diagonal, multivariate normal distribution, $\boldsymbol{z}$ denotes
the redshift array at which said nodes sit, $D_{\rm{P18}}(\boldsymbol{z})$ is
the P18 growth factor evaluated at $\boldsymbol{z}$ and
$L_{C}(K(\boldsymbol{z},\boldsymbol{z}))$ is lower-triangular matrix obtained
from the Cholesky decomposition of the covariance matrix. A way of
interpreting Eq. (16) is as a rotation given by $L_{C}$ on a vector of white
noise $\boldsymbol{v}$ that imposes the correlations of the GP kernel.

This model is similar to the one considered in García-García et al. (2021)
(CGG21 from now on), but using a GP to reconstruct the growth factor instead
of splines. The reasoning behind the choice of GPs over splines for this work
is threefold. First, GPs offer a well-defined measure of the uncertainty in
their predictions which makes assessing the statistical significance of their
results straightforward. Second, GPs are not subject to the strict assumptions
that can bias spline reconstructions such as the choice of linear or cubic
interpolations. However, it is important to bear in mind that GPs are far from
assumption free as their structure is constrained by the properties of the
chosen covariance kernel. Nonetheless, as we will show, these assumptions can
indeed be neglected when the data is constraining enough. Third, GPs are as
differentiable as their covariance matrix kernel (Rasmussen & Williams, 2006).
The derivative of a GP
$g(\boldsymbol{x})\sim\mathcal{N}(m(\boldsymbol{x}),K(\boldsymbol{x},\boldsymbol{x}^{\prime}))$
is another GP given by
$\dot{g}(\boldsymbol{x})\sim\mathcal{N}(\partial_{\boldsymbol{x}}m(\boldsymbol{x}),\partial_{\boldsymbol{x}}\partial_{\boldsymbol{x^{\prime}}}K(\boldsymbol{x},\boldsymbol{x}^{\prime}))$.
This is a very desirable feature that will allow us to use growth rate
measurements (i.e. measurements of the gradient of the growth factor) to
further constrain the reconstructed growth. As it will be shown, the growth
rate measurements will highly restrict the evolution of the growth.

Figure 6: Shows the $f\\!\sigma_{\rm{8}}$ data points from the different surveys used in this work across redshift. Numerical values can be found in Tab. 2 Table 2: Lists the RSD data used in this analysis. Data set |  Redshifts |  Data Points
---|---|---
BOSS DR12 Alam et al. (2017) |  0.38 - 0.61 |  3
eBOSS DR16 Alam et al. (2021) |  1.48 |  1
Wigglez Blake et al. (2012) |  0.44 - 0.73 |  3
VIPERS Pezzotta et al. (2017) |  0.60 - 0.86 |  2
6dF Beutler et al. (2012) |  0.067 |  1
FastSound Okumura et al. (2016) |  1.4 |  1

In order to constrain this model we use the south data collection described in
CGG21 composed of the 3x2 DES Y1 data, the auto-correlation of eBOSS DR16
quasars and the cross-correlation of CMB lensing data with eBOSS DR16 quasars,
DES Y1 clustering and DES Y1 weak lensing data. Therefore, our analysis
combines a total of 7 two-point statistics. Thus we will refer to it as
“7$\times$2-pt” hereafter. For a detailed description of these data we refer
the reader to Sect. 3 of CGG21. Besides the aforementioned data, we also
consider a collection of redshift-space distortion (RSD) measurements by the
BOSS DR12, eBOSS DR16, Wigglez, 6dF and VIPERS surveys. We summarize these
data in Tab. 2 and plot it in Fig. 6. For a full description of these data we
refer the reader to Ruiz-Zapatero et al. (2022b).

These data allow us to constraint the growth factor and its derivative across
redshift. In order to do so we relate $f\sigma_{\rm{8}}(z)$ with the
reconstructed $D(z)$ by

$f\sigma_{\rm{8}}(z)=-\sigma_{8\,P18}(1+z)\frac{\partial D(z)}{\partial z}\,,$
(17)

where $\sigma_{8\,P18}$ is the $\sigma_{8}$ of the $\Lambda$CDM P18 cosmology
and $D(z)$ is given by Eq. 16. The main difference in Eq. 17 with respect to
the $\Lambda$CDM model is that $\sigma_{8}$ is kept fixed. This is because in
our $\Lambda$CDM analysis we defined $D(z=0)=1.0$ such that the amplitude of
the growth factor was varies by $\sigma_{\rm{8}}$. Unlike in $\Lambda$CDM, in
our GP model the amplitude of the growth factor is not fixed, but varies with
each GP realisation. Therefore, the role of $\sigma_{\rm{8}}$ would be
completely degenerate with the amplitude of the $z=0$ GP node and hence is
kept fixed.

Once we have drawn an expression for $f\sigma_{\rm{8}}(z)$ in our GP model,
the next question is how to evaluate it. While this might seem straight
forward, the gradient of the growth factor is extremely sensitive to numerical
errors which can lead to non-physical predictions for $f\sigma_{\rm{8}}(z)$.
One solution to this problem would be to compute the gradient analytically
using the properties of GPs as discussed above. However, in practice
evaluating the derivative of the GP analytically involves treating said
derivative as a second GP with a cross-covariance matrix with the original GP.
Computing and inverting the cross-correlation matrix between the original GP
and its derivative at every step of the Markov chain for the number of nodes
used in this analysis proved too computationally expensive.

The solution found was to use the noise properties of the GP to interpolate
the GP to a finer grid such that a finite difference can be taken to obtain
the gradient to enough precision. Effectively, this is done by applying a
Wiener filter to the original GP where the kernel used for the filtering is
the covariance matrix of the GP. In this scenario the Wiener filter is just
given by
$W(\boldsymbol{z}^{*},\boldsymbol{z})=K(\boldsymbol{z}^{*},\boldsymbol{z}^{*})K^{-1}(\boldsymbol{z},\boldsymbol{z})$
where $\boldsymbol{z}^{*}$ is the finer nodes array of redshifts. The
resampled growth factor is then given by
$D(\boldsymbol{z}^{*})=W(\boldsymbol{z}^{*},\boldsymbol{z})D(\boldsymbol{z})$.
While this approach also involves inverting a matrix, this is just the
covariance matrix of the GP which is a much cheaper operation.

Having discussed our model it is time to discuss how theory and observations
are brought together to constrain the model. We built a Gaussian likelihood
assuming the RSD and angular power spectra measurements are completely
uncorrelated such that the final likelihood is the product of the likelihoods
of the individual data types. This is a fair assumption given that the RSD
data is mostly located in the northern hemisphere whereas the angular power
power spectra are located further South, resulting in a small overlap between
the surveys. We then derived constraints for the GP and cosmological
parameters by applying Bayes theorem using the priors displayed in Tab. 3.

Note that the GP hyperparameters $l$ and $\eta$ were kept fixed. On the one
hand, fixing the hyperparameters avoids the volume effects on the posteriors
expected from including these parameters. On the other hand, fixing the
hyperparameters to a particular set of values will introduce certain biases in
the final reconstruction of the growth factor specially in the regions where
data is sparse. For example, fixing the correlation length is expected to
induce spurious oscillations in the reconstructed function outside of the data
range. Nonetheless the decision to keep both parameters fixed was made as
freeing them introduces severe non-Gaussian degeneracies in the posteriors
which pose a challenge even to gradient-based methods such as HMC or NUTS.
This is due to their hierarchical relationship with the GP nodes. Future
studies that wish to undertake a completely model independent reconstruction
will need to explore inference methods specifically tailored to deal with such
non-Gaussianity such as Riemannian HMC (Betancourt, 2012).

The number of parameters listed in Tab. 3 adds up to 128, a dimensionality
which requires gradient-based samplers. For this reason, we employed a Gibbs
sampling set up where GP parameters were sampled by one NUTS sampler, keeping
all other parameters fixed, and then the rest of the parameters (cosmology and
nuisances) were sampled by their own NUTS sampler keeping the GP parameters
fixed to their last sample. Moreover, the chains were started at the DES Y1
$\Lambda$CDM best-fit cosmology with the GP nodes starting from zero. The
Gibbs scheme combined with the starting point proved to greatly increase the
efficiency of the NUTS adapting phase in this high dimensional space.

In total four analyses were run. First, we ran two $\Lambda$CDM analyses, one
of the 7$\times$2-pt data and another of 7$\times$2-pt + RSD data. In these
analyses, we used the $\Lambda$CDM model to predict the growth factor as
opposed to doing a GP reconstruction. Then the same two analyses were ran but
performing the GP reconstruction. In addition to these, we also rerun the
CGG21 analysis of 7$\times$2-pt data using the priors listed in Tab. 3 in
order to be able to compare our results against the Cobaya pipeline. To get a
rough idea of the additional computational complexity involved in sampling
high dimensional parameter spaces with NUTS, obtaining converged constraints
for the $\Lambda$CDM analysis of 7$\times$2-pt both with and withour RSD data
took approximately 24 CPU hours compared to the 240 CPU hours needed when
performing the GP analysis with an extra 101 parameters.

Growth Factor Reconstruction Priors
---
Parameter | Prior | Parameter | Prior
Cosmology | DES Y1 - Sources redshift cal.
$\Omega_{\rm{m}}$ | $U(0.2,0.6)$ | $\Delta z_{\rm{wl}}^{1}$ | $N(-0.001,0.016)$
$\Omega_{\rm{b}}$ | $U(0.028,0.065)$ | $\Delta z_{\rm{wl}}^{2}$ | $N(-0.0019,0.013)$
$h$ | $U(N(70,5),0.55,0.91)$ | $\Delta z_{\rm{wl}}^{3}$ | $N(0.009,0.011)$
$n_{\rm{s}}$ | $U(0.87,1.07)$ | $\Delta z_{\rm{wl}}^{4}$ | $N(-0.018,0.022)$
$\sigma_{\rm{8}}$ | 0.81 |  |
DES Y1 - Lens redshift cal. | Intrinsic Alignments
$\Delta z_{\rm{gc}}^{1}$ | $N(0.0,0.007)$ | $A_{\mathrm{IA,0}}$ | $U(-5,5)$
$\Delta z_{\rm{gc}}^{2}$ | $N(0.0,0.007)$ | $\eta_{\mathrm{IA}}$ | $U(-5,5)$
$\Delta z_{\rm{gc}}^{3}$ | $N(0.0,0.006)$ | eBOSS - Galaxy bias
$\Delta z_{\rm{gc}}^{4}$ | $N(0.0,0.01)$ | $b_{\rm{QSO}}^{i}$ | $U(0.8,5.0)$
DES Y1 - Shear calibration bias | Gaussian Process
$m^{i}$ | $N(0.012,0.023)$ | $\eta$ | 0.2
DES Y1 - Galaxy bias | $l$ | 0.3
$b_{\rm{gc}}^{i}$ | $U(0.8,3.0)$ | $v^{i}$ | $U(N(0,1),-2,2)$
Table 3: Prior distributions for parameters considered for the Growth Factor
Reconstruction. We sample the cosmological parameters keeping
$\sigma_{\rm{8}}$ fixed to avoid degeneracies with the Gaussian process. For
the same reason the Gaussian process hyperparameters are also kept fixed.

####  3.2.1 $\Lambda$CDM Results

We start by analysing our data using a traditional $\Lambda$CDM model where
the growth factor is obtained from the $\Lambda$CDM parameters by solving the
Jeans equation for the matter anisotropies. This allows us to establish a
frame of reference against which to compare the results of the GP
reconstruction. Moreover, it also allows us to ground our analysis by
comparing it with previous analyses of the same combination of data in the
literature, namely CGG21.

In Fig 7, we show the obtained $\Omega_{\rm{m}}$, $\sigma_{\rm{8}}$ and
$S_{\rm{8}}$ posteriors when performing a $\Lambda$CDM analysis of the
7$\times$2-pt data set with and without including additional RSD data (green
and red contours respectively). We also show the contours obtained by the
CGG21 pipeline when analysing the 7$\times$2-pt data using black dashed lines.
Our $\Lambda$CDM analysis of the 7$\times$2-pt data found
$\Omega_{\rm{m}}=0.287\pm 0.027$, $\sigma_{\rm{8}}=0.809\pm 0.045$ and
$S_{\rm{8}}=0.789\pm 0.016$. The reanalysis of CGG21 using the priors shown in
Tab. 3 found $\Omega_{\rm{m}}=0.296\pm 0.028$, $\sigma_{\rm{8}}=0.794\ 0.043$
and $S_{\rm{8}}=0.786\pm 0.015$. Thus we can observe that the constraints
produced by the LimberJack.jl pipeline presented in this work are completely
consistent with the results of the Cobaya pipeline of CGG21. Combining the
7$\times$2-pt with the RSD data we found $\Omega_{\rm{m}}=0.277\pm 0.021$,
$\sigma_{\rm{8}}=0.827\pm 0.034$ and $S_{\rm{8}}=0.793\pm 0.015$. The addition
of RSD data improved the constraints on $\Omega_{\rm{m}}$ and
$\sigma_{\rm{8}}$ by 20%. However, the constraints on $S_{\rm{8}}$ were
largely unaffected. Regardless of whether or not RSD data are included, our
results are in 2 $\sigma$ disagreement with the P18 results which found
$S_{\rm{8}}=0.832\pm 0.013$ (See Tab. 1 of Planck Collaboration et al., 2020).
The full numerical posteriors of the LimberJack.jl $\Lambda$CDM analyses of
7$\times$2-pt and 7$\times$2-pt plus RSD data can be found in the first two
columns of Tab. 4.

Figure 7: Shows the 2D and 1D marginal distributions for the parameters
$\Omega_{\rm{m}}$, $\sigma_{\rm{8}}$ and $S_{\rm{8}}$ of the $\Lambda$CDM
analysis of 7$\times$2-pt data (red) and 7$\times$2-pt plus RSD data (green).
Black dashed contours show the reanalysis of the CGG21 using Tab. 3’s priors.

####  3.2.2 GP Results

We are now in a position to start discussing the GP reconstruction of the
growth factor. We start by considering how introducing the GP affects the
parameter constraints previously discussed in Sect. 3.2.1. In order to
establish a comparison we derive constraints for
$\sigma_{\rm{8}}\equiv\sigma_{\rm{8}}(z=0)$ and $S_{\rm{8}}\equiv
S_{\rm{8}}(z=0)$ from the GP reconstruction of the growth factor. From Eq. 17,
we can see that within the GP model:

$\displaystyle\sigma_{\rm{8}}=\sigma_{\rm{8}}^{\rm{P18}}D(z=0)\,,$ (18)
$\displaystyle
S_{\rm{8}}=\sigma_{\rm{8}}^{\rm{P18}}D(z=0)\sqrt{\Omega_{\rm{m}}/0.3}\,,$ (19)

where $D(z)$ is given by 16.

Figure 8: Lower triangle: shows the 2D marginal distributions for the
parameters $\Omega_{\rm{m}}$, $\sigma_{\rm{8}}$ and $S_{\rm{8}}$ of the
$\Lambda$CDM (red) and GP reconstruction (blue) analyses of 7$\times$2-pt
data. Upper triangle: shows the 2D distributions for the parameters
$\Omega_{\rm{m}}$, $\sigma_{\rm{8}}$ and $S_{\rm{8}}$ of the $\Lambda$CDM
(green) and GP reconstruction (purple) analyses of 7$\times$2-pt data combined
with RSD data. 1D Marginals for all the previous analyses are shown along the
diagonal.

In Fig. 8 we show four different set of posteriors for the parameters
$\Omega_{\rm{m}}$, $\sigma_{\rm{8}}$ and $S_{\rm{8}}$. In the lower triangle
we show the contours obtained when using the $\Lambda$CDM model (red) and the
GP model (blue) to analyse the 7$\times$2-pt data. Similarly, in the upper
triangle we show the contours obtained when using the $\Lambda$CDM model (red)
and the GP model (blue) to analyse the 7$\times$2-pt data in addition to RSD
data. The associated numerical constraints for the GP reconstructions are
$\Omega_{\rm{m}}=0.289\pm 0.026$, $\sigma_{\rm{8}}=0.839\pm 0.067$ and
$S_{\rm{8}}=0.839\pm 0.067$ when 7$\times$2-pt alone is consider and
$\Omega_{\rm{m}}=0.286\pm 0.023$, $\sigma_{\rm{8}}=0.815\pm 0.039$ and
$S_{\rm{8}}=0.793\pm 0.017$ when RSD data is included. The full posteriors can
be found in the last two columns of Tab. 4. Performing the GP reconstruction
of the growth factor greatly reduces the constraining power of the
7$\times$2-pt data. We observe $\sim$30% wider constraints when using the GP
model to analyse the data on average across the three parameters. It is
important to note that the overwhelming majority of the impact occurs in the
parameters explicitly related to the growth factor such as $\sigma_{\rm{8}}$
and $S_{\rm{8}}$ and the linear bias parameters of the different probes. Other
parameters such as $h$ or $n_{s}$ remain virtually unchanged. The $S_{\rm{8}}$
constraint is also centred at a significantly higher value than when
performing a $\Lambda$CDM analysis. This due to the lack of data at $z=0$ to
constrain the GP, which is reflected in the large error bar in the constraint.
When RSD’s are included the degradation in constraining power when performing
the GP reconstruction is much smaller, with constraints only being $\sim$10%
wider. Moreover, the $S_{\rm{8}}$ constraint becomes centred at exactly the
same value as when performing the $\Lambda$CDM analysis of the same data. Due
the larger error bar, we observe that nonetheless the tension with the P18
$\Lambda$CDM value drops to 1.7 sigma.

Figure 9: Shows the evolution across redshift of the reconstructed $D(z)$
(first row), $S_{\rm{8}}(z)$ (second row) and $f\sigma_{\rm{8}}(z)$ (third
row). Results based on 7$\times$2-pt data are shown on the left column while
results combining 7$\times$2-pt plus RSD data are shown in the right column.
In each panel we over plot the prediction of the $\Lambda$CDM model for the
given data as well as the result of the GP reconstruction. Moreover, we also
overplot the P18 $\Lambda$CDM prediction using a black dashed line for
reference. Finally, each panel counts with a subpanel where the we show the
evolution of the $1\sigma$ confidence intervals of the plotted functions
across redshift.

In order to understand how the inclusion of RSD data leads to these changes,
we have to look at the reconstructed growth factor. In each of the rows of
Fig. 9 we show the constraints obtained for $D(z)$, $S_{\rm{8}}(z)$ and
$f\sigma_{\rm{8}}(z)$ (in this order) as functions of redshift. In each panel
we compare the constraints obtained when using the $\Lambda$CDM model against
the GP reconstruction. The left-column panels show analyses where only the
7$\times$2-pt data was used while the right column shows the respective
analyses when RSD data were included. Moreover, each panel has a subpanel
showing the evolution of the $1\sigma$ confidence interval of each function
over redshift.

Let us begin the discussion by focusing on the top panels of Fig. 9 showing
the evolution of the growth factor. In these panels we can see how introducing
a GP to reconstruct the growth factor from the data increases the error bars
in the predictions of $D(z)$ by one to two orders of magnitude when compared
to the $\Lambda$CDM prediction based on the same data. Including RSD data
significantly contributes to constraining $D(z)$ regardless of the model
considered. In the case of $\Lambda$CDM analyses, we observe a 10% reduction
in the standard deviation of $D(z)$ across redshift. The impact of RSD data is
even bigger when we instead use a GP to reconstruct $D(z)$. In this case we
observe a 20% to 30% percent improvement in the constraints depending on the
redshift window. Nonetheless, we observe that in all four cases the obtained
growth factor is compatible with the P18 prediction at all redshifts.

It is important to note that the uncertainty of the $D(z)$ $\Lambda$CDM
prediction actually falls to zero at $z=0$ since in this model the growth
factor is re-scaled to be precisely one at this redshift. Therefore, by
looking at the $\Lambda$CDM prediction of $D(z)$ we are omitting a large
contribution to the uncertainty of the growth factor in this model, its
amplitude. Therefore, it is useful to consider quantities such as
$S_{\rm{8}}(z)$ that do incorporate the uncertainty in the amplitude of the
growth factor in the $\Lambda$CDM model, encapsulated in the parameter
$\sigma_{\rm{8}}$. In the second row of panels of 9 we show the associated
$S_{\rm{8}}(z)$ constraints for the four scenarios considered. In these panels
we can see that the uncertainty in $S_{\rm{8}}(z)$ of the GP reconstruction at
low redshift is actually comparable to that of the $\Lambda$CDM model the
uncertainty in the amplitude is taken into account. However, the uncertainty
of the GP reconstruction quickly increases once the data becomes sparse as
seen in the $D(z)$ panels.

It is also interesting to note that, as expected, fixing the correlation
length of the GP induces oscillations in the reconstructed growth factor and
$S_{\rm{8}}(z)$. This is visible in the GP reconstruction of $D(z)$ based on
7$\times$2-pt data. However, including RSD data nullifies these oscillations.
In order to understand these behaviours we need to look at the third rows of
panels of Fig. 9. In these panels we show the associated $f\sigma_{\rm{8}}(z)$
predictions as a function of redshift. In the left panel we can see how the
oscillations on the growth factor induced by our assumptions on the GP
hyperparameters result in non-physical predictions for $f\sigma_{\rm{8}}(z)$
despite the growth factor itself being perfectly compatible with the P18
$\Lambda$CDM prediction. Including RSD data solves this problem by directly
constraining the possible values of $f\sigma_{\rm{8}}(z)$ given the current
data. These constraints on $f\sigma_{\rm{8}}(z)$ translate into strict demands
for the evolution of the gradient of the growth factor. Therefore, when RSD
data are included, we can see how the oscillatory behaviour present in the
reconstruction from only 7$\times$2-pt data disappears. This means that when
the data is constraining enough the GP reconstruction remains unbiased by the
assumptions introduced by fixing the hyperparameters.

Figure 10: Shows a comparison between the reconstructed growth factor in this
work using a GP (blue) and the CGG21 reconstruction based on splines (black
dots) of the same 7$\times$2-pt data. We also overplot the P18 $\Lambda$CDM
prediction for the growth factor (black dashed line) for reference.

Finally it is relevant to compare the growth factor reconstruction performed
in this work using GPs against the reconstruction performed by CGG21 from the
same data using splines. In Fig. 10 we can observe that our GP reconstruction
stands in good statistical agreement with the four constraints on
$S_{\rm{8}}(z)$ at $z=[0.24,0.53,0.83,1.5]$ obtained by CGG21. This reinforces
the notion that these different model-agnostic approaches are data-driven, and
not biased by their own individual assumptions.

Growth Factor Reconstruction Posteriors
---
Parameter | 7$\times$2-pt ΛCDM | 7$\times$2-pt +RSDΛCDM | 7$\times$2-pt GP | 7$\times$2-pt +RSDGP
$\chi^{2}$ | 658 | 645 | 652 | 635
$\chi^{2}/(d-p)$ | 1.03 | 0.99 | 1.21 | 1.16
$\Omega_{\rm{m}}$ | 0.287 $\pm$ 0.027 | 0.277 $\pm$ 0.021 | 0.289 $\pm$ 0.026 | 0.286 $\pm$ 0.023
$\Omega_{b}$ | 0.045 $\pm$ 0.009 | 0.043 $\pm$ 0.009 | 0.045 $\pm$ 0.009 | 0.044 $\pm$ 0.009
$h$ | 0.725 $\pm$ 0.042 | 0.733 $\pm$ 0.039 | 0.726 $\pm$ 0.04 | 0.731 $\pm$ 0.04
$n_{s}$ | 0.96 $\pm$ 0.056 | 0.965 $\pm$ 0.059 | 0.961 $\pm$ 0.058 | 0.948 $\pm$ 0.06
$\sigma_{\rm{8}}$ | 0.809 $\pm$ 0.045 | 0.827 $\pm$ 0.034 | 0.839 $\pm$ 0.067 | 0.815 $\pm$ 0.039
$S_{\rm{8}}$ | 0.789 $\pm$ 0.016 | 0.793 $\pm$ 0.015 | 0.82 $\pm$ 0.045 | 0.793 $\pm$ 0.017
$b_{\rm{gc}}^{1}$ | 1.456 $\pm$ 0.101 | 1.418 $\pm$ 0.082 | 1.458 $\pm$ 0.1 | 1.444 $\pm$ 0.09
$b_{\rm{gc}}^{2}$ | 1.774 $\pm$ 0.111 | 1.729 $\pm$ 0.084 | 1.84 $\pm$ 0.118 | 1.771 $\pm$ 0.098
$b_{\rm{gc}}^{3}$ | 1.749 $\pm$ 0.102 | 1.705 $\pm$ 0.078 | 1.869 $\pm$ 0.136 | 1.752 $\pm$ 0.095
$b_{\rm{gc}}^{4}$ | 2.13 $\pm$ 0.123 | 2.081 $\pm$ 0.095 | 2.292 $\pm$ 0.212 | 2.128 $\pm$ 0.121
$b_{\rm{gc}}^{5}$ | 2.19 $\pm$ 0.133 | 2.138 $\pm$ 0.106 | 2.307 $\pm$ 0.268 | 2.179 $\pm$ 0.136
$\Delta z_{\rm{gc}}^{1}$ | 0.003 $\pm$ 0.007 | 0.004 $\pm$ 0.007 | 0.003 $\pm$ 0.007 | 0.004 $\pm$ 0.007
$\Delta z_{\rm{gc}}^{2}$ | 0.003 $\pm$ 0.007 | 0.003 $\pm$ 0.007 | 0.002 $\pm$ 0.007 | 0.003 $\pm$ 0.007
$\Delta z_{\rm{gc}}^{3}$ | 0.001 $\pm$ 0.006 | 0.001 $\pm$ 0.006 | 0.0 $\pm$ 0.006 | 0.001 $\pm$ 0.006
$\Delta z_{\rm{gc}}^{4}$ | 0.002 $\pm$ 0.009 | 0.003 $\pm$ 0.01 | 0.002 $\pm$ 0.01 | 0.002 $\pm$ 0.009
$\Delta z_{\rm{gc}}^{5}$ | -0.001 $\pm$ 0.009 | -0.001 $\pm$ 0.01 | -0.001 $\pm$ 0.01 | -0.001 $\pm$ 0.01
$m_{\rm{wl}}^{1}$ | 0.02 $\pm$ 0.023 | 0.022 $\pm$ 0.023 | 0.016 $\pm$ 0.022 | 0.022 $\pm$ 0.022
$m_{\rm{wl}}^{2}$ | 0.014 $\pm$ 0.021 | 0.009 $\pm$ 0.022 | 0.012 $\pm$ 0.022 | 0.009 $\pm$ 0.022
$m_{\rm{wl}}^{3}$ | 0.011 $\pm$ 0.021 | 0.009 $\pm$ 0.021 | 0.012 $\pm$ 0.02 | 0.009 $\pm$ 0.02
$m_{\rm{wl}}^{4}$ | 0.003 $\pm$ 0.022 | 0.006 $\pm$ 0.02 | 0.007 $\pm$ 0.021 | 0.006 $\pm$ 0.021
$\Delta z_{\rm{wl}}^{1}$ | 0.006 $\pm$ 0.012 | 0.009 $\pm$ 0.013 | 0.003 $\pm$ 0.013 | 0.009 $\pm$ 0.013
$\Delta z_{\rm{wl}}^{2}$ | -0.03 $\pm$ 0.01 | -0.035 $\pm$ 0.011 | -0.031 $\pm$ 0.011 | -0.035 $\pm$ 0.011
$\Delta z_{\rm{wl}}^{3}$ | 0.006 $\pm$ 0.01 | 0.007 $\pm$ 0.01 | 0.009 $\pm$ 0.01 | 0.007 $\pm$ 0.01
$\Delta z_{\rm{wl}}^{4}$ | -0.025 $\pm$ 0.018 | -0.025 $\pm$ 0.019 | -0.02 $\pm$ 0.02 | -0.024 $\pm$ 0.018
$b_{\rm{QSO}}^{1}$ | 2.238 $\pm$ 0.19 | 2.18 $\pm$ 0.169 | 2.179 $\pm$ 0.365 | 2.259 $\pm$ 0.237
$b_{\rm{QSO}}^{2}$ | 2.422 $\pm$ 0.216 | 2.358 $\pm$ 0.187 | 2.283 $\pm$ 0.378 | 2.481 $\pm$ 0.296
$A_{IA}$ | 0.279 $\pm$ 0.17 | 0.31 $\pm$ 0.187 | 0.279 $\pm$ 0.177 | 0.312 $\pm$ 0.183
$\alpha_{IA}$ | 0.168 $\pm$ 2.334 | 0.115 $\pm$ 2.47 | -0.223 $\pm$ 2.538 | 0.013 $\pm$ 2.456
Table 4: Posterior distributions of the Growth Factor Reconstruction analyses.
Error bars denote the $1\sigma$ confidence interval.

##  4 Conclusions

In this work we have presented LimberJack.jl, an auto-differentiable
cosmological code to compute angular power spectra fully written in Julia. The
goal of LimberJack.jl is to enable the use of gradient-based inference methods
in cosmological analyses. LimberJack.jl core strength’s are:

  * •

Auto-Differentiablility: Every step in between the input cosmological
parameters and the output theoretical prediction in LimberJack.jl is
compatible with Julia’s auto-differentiation (AD) libraries ForwardDiff.jl and
ReverseDiff.jl. These methods result in gradients up to an order of magnitude
faster than when using finite differences. This is the key feature that makes
the use of gradient-based inference methods possible with LimberJack.jl

  * •

Speed: LimberJack.jl is equipped with a native implementation of the matter
power spectrum emulator EmuPk (Mootoovaloo et al., 2022) which makes it orders
of magnitude faster than CLASS.

  * •

Accuracy: The LimberJack.jl theoretical predictions have been thoroughly
tested against those of the well-established cosmological code CCL, as well as
the quality of the AD gradients against the more costly finite difference
results.

  * •

Interoperability: Thanks to its modular structure, LimberJack.jl can be easily
interfaced with other Julia libraries to increase its capabilities. For
example LimberJack.jl can be interfaced with the Bolt.jl library (Li &
Sullivan, 2023) to gain access to the first auto-differentiable Boltzmann
code.

Furthermore, we presented two examples of how LimberJack.jl can be employed to
perform present and future cosmological analyses using the gradient-based
inference algorithm NUTS (Hoffman & Gelman, 2011), a self-tuning formulation
of the Hamiltonian Monte Carlo algorithm. In the first example we reproduced
two analyses of DES Y1 3x2-point data performed using the well established
Cobaya library to ensure the reliability of LimberJack.jl. LimberJack.jl’s
constraints were found virtually identical to those obtained with Cobaya
regardless of whether the matter power spectrum was computed using the
Eisenstein and Hu formula or EmuPk (Mootoovaloo et al., 2022) to emulate CLASS
(Lesgourgues, 2011). Moreover, the NUTS sampler proved to be 1.5 times more
efficient (measured as effective samples per likelihood call) than the MH
sampler used by Cobaya. However, this improvement in efficiency proved to be
not enough to compensate for the cost of computing the gradient of the
likelihood despite using AD methods. Further work is necessary to determine
the point at which gradient-based inference methods out-weight the cost of
computing the likelihood gradient in angular power spectra analyses.

In the second example we showcased the unique capabilities of LimberJack.jl by
performing a Gaussian process (GP) reconstruction of the growth factor across
redshift adding to a total of 128 parameters. In order to constrain this model
we employed a combination of DES Y1 galaxy clustering and weak lensing data,
eBOSS QSO’s and CMB lensing (referred to as 7$\times$2-pt in the text) as well
as a collection of the latest RSD measurements of $f\sigma_{\rm{8}}$. We
started by considering a $\Lambda$CDM analysis of the aforementioned data to
establish a baseline for the GP reconstruction. Our $\Lambda$CDM analysis of
the 7$\times$2-pt data found $S_{\rm{8}}=0.789\pm 0.016$. Adding the RSD data
yielded $S_{\rm{8}}=0.793\pm 0.015$. Regardless of whether or not RSD data are
included, our results are in $2\sigma$ disagreement with the Planck 2018
results which found $S_{\rm{8}}=0.832\pm 0.013$ (See Tab. 1 of Planck
Collaboration et al., 2020). Performing the GP reconstruction instead yielded
$S_{\rm{8}}=0.839\pm 0.067$ when 7$\times$2-pt alone is considered and
$S_{\rm{8}}=0.793\pm 0.017$ when RSD data is included. Looking at the
reconstructed growth factor we observed a reasonable agreement between the GP
and the Planck 2018 $\Lambda$CDM prediction regardless of the data combination
used. However, including RSD data significantly smoothed the reconstruction of
the growth factor, disfavouring large oscillations. Moreover, it improved the
constraints on the reconstructed growth factor by 20% on average across
redshift. This stresses the importance of including RSD data in future
cosmological analyses, specially given the up coming DESI survey (DESI
Collaboration et al., 2016).

Auto-differentiable and gradient-based inference methods will play a crucial
role speeding up future cosmological analyses as well as enabling entirely new
science. For instance, analyses of multiple auto- and cross-correlations
between stage-IV surveys may contain up to a hundred free parameters. This
will be the case even when performing traditional $\Lambda$CDM analyses with
minimalistic modelling of systematics simply due to the large number of
tomographic bins that will be involved. Future surveys will, however, provide
unprecedented measurements of small scales. In order to fit these scales and
further our understanding of non-linear cosmology, more complex modelling of
baryonic effects will have to be included. Similarly, the constraining power
provided by the new data will enable analyses with beyond $\Lambda$CDM physics
as well as model-agnostic reconstructions such as the one presented in CGG21
and this work which will have a similar effect in the number of free
parameters. In addition to this, gradient-based inference methods are already
indispensable to undertake field level inference cosmology and they will
become more so in the future.

Finally, while LimberJack.jl is already a fully functional tool, there are
several avenues for future improvement:

  1. 1.

Improved predictions: the methods currently implemented in LimberJack.jl
provide enough accuracy to analyse DES Y1 data. However, these methods will
need to be improved in order to analyse DES Y3 data or future data sets such
as LSST. Here we present a non-exhaustive list of possible extensions:

     * •

Non-linear corrections to the matter power spectrum beyond the Halofit
formula.

     * •

Small scale baryonic effects on the galaxy and matter power spectra.

     * •

Angular power spectra beyond the Limber approximation.

     * •

Scale-dependent growth of structure.

     * •

Perturbatory expansion models for the matter-galaxy bias (Lepori et al.,
2023).

  2. 2.

Parallelization: currently the threading parallelisation of LimberJack.jl is
suboptimal. This is because the default Julia threading parallelisation scheme
does not handle the shared memory between the threads efficiently enough. At
the moment computing resources are best spent running different instances of
LimberJack.jl in parallel as opposed to parallelising one instance. Future
works could study alternative parallelisation schemes or manually managing the
memory between the threads to improve the multi-core performance of
LimberJack.jl.

  3. 3.

GPU’s: LimberJack.jl currently cannot run on GPU’s which are known to
significantly speed-up cosmological inference. This is due to LimberJack.jl
performing scalar indexing operations at several points of the computation of
angular power spectra. Future work could study how to bypass these operations
and to make LimberJack.jl comptible with Julia’s GPU libraries such as
CUDA.jl.

  4. 4.

Backwards-AD: currently LimberJack.jl’s preferred AD mode is forward-AD.
However, statistical inference preferred AD mode is backwards-AD, specially as
the number of parameters increases. Future works could look into making
LimberJack.jl compatible with the latest Julia AD libraries such as Zygote.jl
or Enzyme.jl to implement efficient backwards-AD.

## Author contributions

  * •

Jaime Ruiz-Zapatero: Science lead, lead designer and programmer.

  * •

David Alonso: Tracers and Limber integral design. Main advisor.

  * •

Carlos García-García: Science advisor and Cobaya analyses.

  * •

Andrina Nicola: Halofit code.

  * •

Arrykrishna Mootoovaloo: EmuPk integration.

  * •

James M. Sullivan: Bolt.jl integration.

  * •

Marco Bonici: Design advisor and benchmarks.

  * •

Pedro G. Ferreira: Science advisor.

## Software Availability

LimberJack.jl is a fully open-source project available on Julia’s general
repository of packages and GitHub .

## Acknowledgements

We would like to thank Francois Lanusse and Marius Millea for their helpful
comments. JRZ is supported by an STFC doctoral studentship. DA is supported by
the Beecroft Trust. CGG is supported by the Beecroft Trust. AN is supported
through NSF grants AST-1814971 and AST-2108126 AM is funded through Grant
62192 from the John Templeton Foundation to LSST Corporation. PGF is supported
by STFC and the Beecroft Trust. JMS is partially supported by a US Department
of Energy Office of Science Graduate Student Researcher award.

The opinions expressed in this publication are those of the authors and do not
necessarily reflect the views of LSSTC or the John Templeton Foundation.

For the purpose of Open Access, the authors have applied a CC BY public
copyright licence to any Author Accepted Manuscript version arising from this
submission.

We made extensive use of the numpy (Van Der Walt et al., 2011), scipy
(Virtanen et al., 2020), astropy (Astropy Collaboration et al., 2013, 2018),
healpy (Zonca et al., 2019), GetDist (Lewis, 2019), and matplotlib (Hunter,
2007) python packages. We also make use of the Julia packages ForwardDiff.jl
(Revels et al., 2016), AdvancedHMC.jl (Xu et al., 2020) and Turing.jl (Ge et
al., 2018).

##  Appendix A Tutorial

Figure 11: Shows an schematic of LimberJack.jl’s code.

In this section we will show a small demonstration of how to use
LimberJack.jl. We shall begin by having a quick overview of the structure of
LimberJack.jl. This will help us understand how to use the code later on. As
described in Sect. 2.2, LimberJack.jl is a modular code. A summary of the
relationship between the different modules can be found in Fig. 11. The core
of LimberJack.jl is the Cosmology structure and its homonymous constructor.
Cosmology’s role is to contain all the information on how to compute the
theoretical predictions. When the Cosmology structure is initiated it computes
theoretical predictions for the expansion history, comoving distance, the
growth factor and the matter power spectrum according to the provided
prescriptions. The computation of background quantities occurs within core.jl
itself. However, the computation of the matter power spectrum takes place
across three different modules, boltzmann.jl, growth.jl and halofit.jl, which
compute the primordial matter power spectrum, the linear growth factor and the
non-linear corrections respectively. The different predictions are evaluated
at a grid of values which are then used to build interpolators for each of the
respective quantities. The interpolators are then stored inside the Cosmology
structure. This allows the user to quickly compute the theoretical predictions
at arbitrary values by using a series of public functions. These functions
take the Cosmology structure and the necessary inputs. Inside these functions
the corresponding interpolator is evaluated at the provided inputs to return
the prediction to the user.

Let us briefly showcase how this is done in code A by computing some
relatively straight forward predictions:

⬇

1# Install LimberJack

2using Pkg

3Pkg.add("LimberJack")

4

5# Import it

6using LimberJack

7

8# create Cosmology instance

9cosmo = Cosmology()

10# same as

11# cosmology = Cosmology($\Omega_{\rm{m}}=0.3$, $\Omega_{b}=0.05$,

12# h=0.67, ns=0.96, $\sigma_{\rm{8}}=0.81$)

13

14# Compute predictions

15zs = [0.1, 0.5, 1.0, 3.0]

16H = cosmo.cpar.h*100*Ez(cosmo, zs)

17chi = comoving_radial_distance(cosmo, zs)

18Dz = growth_factor(cosmo, zs)

19fz = growth_rate(cosmo, zs)

20fs8z = fs8(cosmo, zs)

Computing the matter power spectrum is just as easy (see code A). However, the
computation is subject to a series of options that the user can alter. By
default Cosmology will use the E&H formula to find $P_{0}$ and it will not
apply non-linear corrections. These settings can be changed by specifying the
keyword arguments tk_mode and Pk_mode. In terms of transfer functions
LimberJack.jl offers two possibilities tk_mode = :EisHu (default) / :EmuPk
which correspond to using the Eisenstein and Hu formula or EmuPk. Similarly,
Pk_mode = :linear / :Halofit which determines whether or not non-linear
corrections are applied using Halofit. LimberJack.jl offers two distinct
public functions to evaluate either the linear or non-linear matter power
spectrum regardless of the choice in Pk_mode. However, if Pk_mode = :linear
the two functions will return the linear matter power spectrum.

⬇

1using NPZ

2

3cosmo_lin_EisHu = Cosmology($\Omega_{\rm{m}}=0.25$, $\Omega_{b}=0.03$,

4 h=0.70, ns=1.0,

5 $\sigma_{\rm{8}}=0.78$,

6 tk_mode=:EisHu,

7 Pk_mode=:linear)

8

9cosmo_nonlin_EisHu = Cosmology($\Omega_{\rm{m}}=0.25$, $\Omega_{b}=0.03$,

10 h=0.70, ns=1.0,

11 $\sigma_{\rm{8}}=0.78$,

12 tk_mode=:EisHu,

13 Pk_mode=:Halofit)

14

15cosmo_nonlin_emupk = Cosmology($\Omega_{\rm{m}}=0.25$, $\Omega_{b}=0.03$,

16 h=0.70, ns=1.0,

17 $\sigma_{\rm{8}}=0.78$,

18 tk_mode=:EmuPk,

19 Pk_mode=:Halofit)

20

21zs = [0.1, 0.3, 0.5]

22ks = [100, 300, 1000]

23lin_eh_Pks = lin_Pk(cosmo_lin_EisHu, ks, zs) = nonlin_Pk(cosmo_lin_EisHu,
ks, zs)

24nonlin_eh_Pks = nonlin_Pk(cosmo_nonlin_EisHu, ks, zs)

25nonlin_emupk_Pks = nonlin_Pk(cosmo_nonlin_emupk, ks, zs)

Computing angular power spectra is a slightly more involved process. An
example of how this is done in LimberJack.jl can be found in code A. In this
code we can see that first a Cosmology structure must be initiated. The
Cosmology structure automatically computes the matter power spectrum given the
user specifications. Then the user must compute the radial kernels of the
relevant tracers by proving the corresponding LimberJack.jl public functions
with the Cosmology structure and the distribution of sources. Moreover, as
described in 2.1 different tracers can be impacted by different systematic
effects. These systematics are accounted by incorporating a series of nuisance
parameters that can also be provided to the tracers public functions of
LimberJack.jl. The output of the tracer functions (see lines 10, 15 and 18 of
code A) is a Tracer structure that hosts an interpolator for the corresponding
radial kernel and the corrections described in Eqns. 9 and 10 (set to one in
the case of clustering). The angular power spectra can computed by providing
the AngularCls public function with the aforementioned Cosmology and Tracer
objects as well as the desired multipoles.

⬇

1# Initiate cosmology

2cosmo = Cosmology()

3

4# Define a distribution of sources

5z = Vector(range(0., stop=2., length=256))

6nz = @. exp(-0.5*((z-0.5)/0.05)^2)

7

8# Create tracer objects

9bias = 1.0 # Galaxy-mass bias

10tg = NumberCountsTracer(cosmo, z, nz; b=bias)

11

12mbias = 0.0 # shape multiplicative bias

13A_IA = 0.0 # Amplitude of intrinsic alignments power spectrum

14alpha_IA = 0.0 # Slope of intrinsic alignments power spectrum

15ts = WeakLensingTracer(cosmo, z, nz;

16 m=mbias, IA_params=[A_IA, alpha_IA])

17

18tk = CMBLensingTracer(cosmo)

19

20# Compute power spectra

21ls = [10.0, 30.0, 100.0, 300.0, 1000.0]

22Cl_gg = angularCls(cosmo, tg, tg, ls)

23Cl_gs = angularCls(cosmo, tg, ts, ls)

24Cl_ss = angularCls(cosmo, ts, ts, ls)

25Cl_gk = angularCls(cosmo, tg, tk, ls)

26Cl_sk = angularCls(cosmo, ts, tk, ls)

Finally, in code A we show how to use Turing.jl in unison with LimberJack to
build and sample an statistical model for a the DES Y1 3x2 analysis. The first
step is to load the data. For this purpose we will use the libraries YAML and
sacc, ubiquitous in astrophysics. Julia counts with a native implementation of
YAML but not of sacc. However, calling Python libraries from Julia is
extremely simple thanks to the PythonCall.jl library. PythonCall.jl allows us
to import sacc.py as a Julia module and read files entirely within Julia. Note
that in order to this we must first install sacc in the Python environment of
Julia or point PythonCall to our local Pytthon installation. Instructions on
how to do this can be found in the LimberJack.jl GitHub. Once the data are
loaded they must be passed to the LimberJack.jl public function make_data
which turns the files into Julia structures that LimberJack.jl can easily
manage.

After that, the user can use Turing.jl’s model macro to define an statistical
model. Inside the model, the user must define the priors for the parameters of
the model using the Distributions.jl API. Note that while the cosmology
parameters can be directly passed to the Cosmology structure constructor the
nuisance parameters must be stored inside a Julia dictionary. The name of
these parameters inside the dictionary must follow a strict convention
”tracer_name \+ __ + bin_number\+ _ +nuisance_parameter_name”.

In order to obtain the theoretical prediction for the DES Y1 3x2 analysis data
vector the user must provide the just initiated Cosmology structure and the
meta and files structures generated by make_data to the public function
Theory. Moreover, the user must also provide the nuisance parameter dictionary
using the Nuisances keyword argument. The Theory function orchestrates the
computation of the theory vector, computing the necessary radial kernels and
evaluating the angular power spectra in the correct order.

Once a theoretical prediction has been obtained, the user must define how the
data is distributed with respect to the theory prediction. In the case of a
Gaussian likelihood, this corresponds to a multivariate Gaussian distribution
with mean the theory prediction and covariance matrix the data covariance
matrix. All is left is to do is to use Turing.jl to condition the model on the
observations (see line 80 of code A), define the sampler we wish to use and
sample the model. For a more thorough explanation of how to use Turing.jl and
its different options please see Turing.jl’s documentation. The results can be
found plotted in Fig. 5

⬇

1# Imports

2using LinearAlgebra

3using Turing

4using LimberJack

5using YAML

6using PythonCall

7sacc = pyimport("sacc");

8

9# Load data

10sacc_path = "cls_FD_covG.fits"

11yaml_path = "DESY1.yml"

12sacc_file = sacc.Sacc().load_fits(sacc_path)

13yaml_file = YAML.load_file(yaml_path)

14meta, files = make_data(sacc_file, yaml_file)

15data = meta.data

16cov = meta.cov

17

18# Define model

19@model function model(data;

20 meta=meta,

21 files=files)

22 $\Omega_{\rm{m}}$ ~ Uniform(0.2, 0.6)

23 $\Omega_{b}$ ~ Uniform(0.028, 0.065)

24 h ~ TruncatedNormal(0.72, 0.05, 0.64, 0.82)

25 $\sigma_{\rm{8}}$ ~ Uniform(0.4, 1.2)

26 ns ~ Uniform(0.84, 1.1)

27

28 DESgc__0_b ~ Uniform(0.8, 3.0)

29 DESgc__1_b ~ Uniform(0.8, 3.0)

30 DESgc__2_b ~ Uniform(0.8, 3.0)

31 DESgc__3_b ~ Uniform(0.8, 3.0)

32 DESgc__4_b ~ Uniform(0.8, 3.0)

33 DESgc__0_dz ~ TruncatedNormal(0.0, 0.007, -0.2, 0.2)

34 DESgc__1_dz ~ TruncatedNormal(0.0, 0.007, -0.2, 0.2)

35 DESgc__2_dz ~ TruncatedNormal(0.0, 0.006, -0.2, 0.2)

36 DESgc__3_dz ~ TruncatedNormal(0.0, 0.01, -0.2, 0.2)

37 DESgc__4_dz ~ TruncatedNormal(0.0, 0.01, -0.2, 0.2)

38 DESwl__0_dz ~ TruncatedNormal(-0.001, 0.016, -0.2, 0.2)

39 DESwl__1_dz ~ TruncatedNormal(-0.019, 0.013, -0.2, 0.2)

40 DESwl__2_dz ~ TruncatedNormal(0.009, 0.011, -0.2, 0.2)

41 DESwl__3_dz ~ TruncatedNormal(-0.018, 0.022, -0.2, 0.2)

42 DESwl__0_m ~ Normal(0.012, 0.023)

43 DESwl__1_m ~ Normal(0.012, 0.023)

44 DESwl__2_m ~ Normal(0.012, 0.023)

45 DESwl__3_m ~ Normal(0.012, 0.023)

46 A_IA ~ Uniform(-5, 5)

47 alpha_IA ~ Uniform(-5, 5)

48

49 nuisances = Dict("DESgc__0_b" => DESgc__0_b,

50 "DESgc__1_b" => DESgc__1_b,

51 "DESgc__2_b" => DESgc__2_b,

52 "DESgc__3_b" => DESgc__3_b,

53 "DESgc__4_b" => DESgc__4_b,

54 "DESgc__0_dz" => DESgc__0_dz,

55 "DESgc__1_dz" => DESgc__1_dz,

56 "DESgc__2_dz" => DESgc__2_dz,

57 "DESgc__3_dz" => DESgc__3_dz,

58 "DESgc__4_dz" => DESgc__4_dz,

59 "DESwl__0_dz" => DESwl__0_dz,

60 "DESwl__1_dz" => DESwl__1_dz,

61 "DESwl__2_dz" => DESwl__2_dz,

62 "DESwl__3_dz" => DESwl__3_dz,

63 "DESwl__0_m" => DESwl__0_m,

64 "DESwl__1_m" => DESwl__1_m,

65 "DESwl__2_m" => DESwl__2_m,

66 "DESwl__3_m" => DESwl__3_m,

67 "A_IA" => A_IA,

68 "alpha_IA" => alpha_IA,)

69

70 cosmology = Cosmology($\Omega_{\rm{m}}$=$\Omega_{\rm{m}}$,
$\Omega_{b}$=$\Omega_{b}$,

71 h=h, ns=ns, $\sigma_{\rm{8}}$=$\sigma_{\rm{8}}$,

72 tk_mode=:EisHu,

73 Pk_mode=:Halofit)

74

75 theory = Theory(cosmology, meta, files; Nuisances=nuisances)

76 data ~ MvNormal(theory, cov)

77end

78

79# Condition model on data

80cond_model = model(data)

81

82# Define sampler

83nadapts = 500

84TAP = 0.65

85sampler = NUTS(adaptation, TAP)

86

87# Sample model

88iterations = 1000

89chain = sample(cond_model, sampler, iterations)

## References

  * Abazajian et al. (2016) Abazajian K. N., et al., 2016, arXiv e-prints, p. arXiv:1610.02743
  * Abbott et al. (2018a) Abbott T. M. C., et al., 2018a, Phys. Rev. D, 98, 043526
  * Abbott et al. (2018b) Abbott T. M. C., et al., 2018b, Phys. Rev. D, 98, 043526
  * Ade et al. (2019) Ade P., et al., 2019, J. Cosmology Astropart. Phys, 2019, 056
  * Alam et al. (2017) Alam S., et al., 2017, MNRAS, 470, 2617
  * Alam et al. (2021) Alam S., et al., 2021, Phys. Rev. D, 103, 083533
  * Angulo et al. (2021) Angulo R. E., Zennaro M., Contreras S., Aricò G., Pellejero-Ibañez M., Stücker J., 2021, MNRAS, 507, 5869
  * Aricò et al. (2023) Aricò G., Angulo R. E., Zennaro M., Contreras S., Chen A., Hernández-Monteagudo C., 2023, arXiv e-prints, p. arXiv:2303.05537
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, aj, 156, 123
  * Bardeen et al. (1986) Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15
  * Betancourt (2012) Betancourt M. J., 2012, arXiv e-prints, p. arXiv:1212.4693
  * Beutler et al. (2012) Beutler F., et al., 2012, mnras, 423, 3430
  * Blake et al. (2012) Blake C., et al., 2012, MNRAS, 425, 405
  * Bonici et al. (2023) Bonici M., Bianchini F., Ruiz-Zapatero J., 2023, arXiv e-prints, p. arXiv:2307.14339
  * Campagne et al. (2023) Campagne J.-E., et al., 2023, The Open Journal of Astrophysics, 6, 15
  * Chisari et al. (2019) Chisari N. E., et al., 2019, ApJS, 242, 2
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036
  * Dodelson (2003) Dodelson S., 2003, Modern Cosmology
  * Donald-McCann et al. (2022) Donald-McCann J., Beutler F., Koyama K., Karamanis M., 2022, Mon. Not. Roy. Astron. Soc., 511, 3768
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, ApJ, 496, 605
  * Eisenstein & Hu (1999) Eisenstein D. J., Hu W., 1999, ApJ, 511, 5
  * Energy Survey et al. (2023) Energy Survey D., et al., 2023, arXiv e-prints, p. arXiv:2305.17173
  * García-García et al. (2021) García-García C., Ruiz Zapatero J., Alonso D., Bellini E., Ferreira P. G., Mueller E.-M., Nicola A., Ruiz-Lapuente P., 2021, arXiv e-prints, p. arXiv:2105.12108
  * Ge et al. (2018) Ge H., Xu K., Ghahramani Z., 2018, in International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain. pp 1682–1690, http://proceedings.mlr.press/v84/ge18b.html
  * Gelman et al. (2013) Gelman A., Carlin J., Stern H., Dunson D., Vehtari A., Rubin D., 2013, Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science, Taylor & Francis, https://books.google.co.uk/books?id=ZXL6AQAAQBAJ
  * Griewank & Walther (2008) Griewank A., Walther A., 2008, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, second edn. Society for Industrial and Applied Mathematics, USA
  * Harrison (2021) Harrison D., 2021, arXiv e-prints, p. arXiv:2110.06209
  * Hirata & Seljak (2004) Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 063526
  * Hoffman & Gelman (2011) Hoffman M. D., Gelman A., 2011, arXiv e-prints, p. arXiv:1111.4246
  * Hoffmann (2014) Hoffmann P. H. W., 2014, arXiv e-prints, p. arXiv:1411.0583
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lepori et al. (2023) Lepori F., Schulz S., Adamek J., Durrer R., 2023, J. Cosmology Astropart. Phys, 2023, 036
  * Lesgourgues (2011) Lesgourgues J., 2011, arXiv e-prints, p. arXiv:1104.2932
  * Lewis (2019) Lewis A., 2019, arXiv e-prints, p. arXiv:1910.13970
  * Li & Sullivan (2023) Li Z., Sullivan J., 2023, In prep,
  * Limber (1953) Limber D. N., 1953, ApJ, 117, 134
  * LoVerde & Afshordi (2008) LoVerde M., Afshordi N., 2008, Phys. Rev. D, 78, 123506
  * MacKay (2002) MacKay D. J. C., 2002, Information Theory, Inference & Learning Algorithms. Cambridge University Press, Cambridge University Press, Shaftesbury Road Cambridge, CB2 8BS, United Kingdom
  * Margossian (2018) Margossian C. C., 2018, arXiv e-prints, p. arXiv:1811.05031
  * Margossian & Betancourt (2021) Margossian C. C., Betancourt M., 2021, arXiv e-prints, p. arXiv:2112.14217
  * Mead et al. (2021) Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401
  * Mootoovaloo et al. (2022) Mootoovaloo A., Jaffe A. H., Heavens A. F., Leclercq F., 2022, Astronomy and Computing, 38, 100508
  * Okumura et al. (2016) Okumura T., et al., 2016, PASJ, 68, 38
  * Pezzotta et al. (2017) Pezzotta A., et al., 2017, aap, 604, A33
  * Piras & Spurio Mancini (2023) Piras D., Spurio Mancini A., 2023, arXiv e-prints, p. arXiv:2305.06347
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Rasmussen & Williams (2006) Rasmussen C. E., Williams C. K. I., 2006, Gaussian processes for machine learning.. Adaptive computation and machine learning, MIT Press
  * Revels et al. (2016) Revels J., Lubin M., Papamarkou T., 2016, arXiv:1607.07892 [cs.MS]
  * Ruiz-Zapatero et al. (2022a) Ruiz-Zapatero J., Alonso D., Ferreira P. G., Garcia-Garcia C., 2022a, Phys. Rev. D, 106, 083523
  * Ruiz-Zapatero et al. (2022b) Ruiz-Zapatero J., García-García C., Alonso D., Ferreira P. G., Grumitt R. D. P., 2022b, mnras, 512, 1967
  * Sehgal et al. (2019) Sehgal N., et al., 2019, in Bulletin of the American Astronomical Society. p. 6 (arXiv:1906.10134), doi:10.48550/arXiv.1906.10134
  * Seikel et al. (2012) Seikel M., Clarkson C., Smith M., 2012, jcap, 2012, 036
  * Smith et al. (2003) Smith R. E., et al., 2003, MNRAS, 341, 1311
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757
  * Spurio Mancini et al. (2022) Spurio Mancini A., Piras D., Alsing J., Joachimi B., Hobson M. P., 2022, MNRAS, 511, 1771
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
  * Troxel et al. (2018) Troxel M. A., et al., 2018, Phys. Rev. D, 98, 043528
  * Van Der Walt et al. (2011) Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science & Engineering, 13, 22
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Xu et al. (2020) Xu K., Ge H., Tebbutt W., Tarek M., Trapp M., Ghahramani Z., 2020, in Symposium on Advances in Approximate Bayesian Inference. pp 1–10
  * Zhang et al. (2021) Zhang L., Carpenter B., Gelman A., Vehtari A., 2021, arXiv e-prints, p. arXiv:2108.03782
  * Zonca et al. (2019) Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, Journal of Open Source Software, 4, 1298
  * van Merriënboer et al. (2018) van Merriënboer B., Breuleux O., Bergeron A., Lamblin P., 2018, arXiv e-prints, p. arXiv:1810.11530