The Open Journal of Astrophysics submitted XXX; accepted YYY

††thanks: ⋆ E-mail<EMAIL_ADDRESS>
$\dagger$ E-mail<EMAIL_ADDRESS>
$\ddagger$ E-mail<EMAIL_ADDRESS>
#  Capse.jl: Efficient and Auto-Differentiable CMB power spectra emulation

Marco Bonici⋆1,2,3 Federico Bianchini†4,5,6 Jaime Ruiz-Zapatero‡7 1 INAF-IASF
Milano IT  2 Waterloo Centre for Astrophysics, University of Waterloo,
Waterloo, ON N2L 3G1, Canada  $3$ Department of Physics and Astronomy,
University of Waterloo, Waterloo, ON N2L 3G1, Canada  4 Kavli Institute for
Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall,
Stanford, CA, 94305, USA  5SLAC National Accelerator Laboratory, 2575 Sand
Hill Road, Menlo Park, CA 94025  6Department of Physics, Stanford University,
382 Via Pueblo Mall, Stanford, CA 94305  7Astrophysics, University of Oxford,
DWB, Keble Road, Oxford OX1 3RH, UK

###### Abstract

We present Capse.jl, a novel neural network-based emulator designed for rapid
and accurate prediction of Cosmic Microwave Background (CMB) temperature,
polarization, and lensing angular power spectra. The emulator computes
predictions in just a few microseconds with emulation errors below $0.1\sigma$
for all the scales relevant for the upcoming CMB-S4 survey. Capse.jl can also
be trained in an hour’s time on a 8-cores CPU. We test Capse.jl on Planck
2018, ACT DR4, and 2018 SPT-3G data and demonstrate its capability to derive
cosmological constraints comparable to those obtained by traditional methods,
but with a computational efficiency that is three to six orders of magnitude
higher. We take advantage of the differentiability of our emulators to use
gradient-based methods, such as Pathfinder and Hamiltonian Monte Carlo (HMC),
which speed up the convergence and increase sampling efficiency. Together,
these features make Capse.jl a powerful tool for studying the CMB and its
implications for cosmology. When using the fastest combination of our
likelihoods, emulators, and analysis algorithm, we are able to perform a
Planck $TT+TE+EE$ analysis in less than a second. To ensure full
reproducibility, we provide open access to the codes and data required to
reproduce all the results of this work

###### keywords:

Cosmology: Cosmic Microwave Background – Methods: statistical, data analysis

##  1 Introduction

The experimental cosmological landscape is extremely rich and promising.
Upcoming millimeter-wave telescopes, including the Simons Observatory (SO, Ade
et al., 2019) and CMB-S4 (Abazajian et al., 2019), are building on the legacy
of the Planck satellite (Aghanim et al., 2020a), the South Pole Telescope
(SPT, Sobrin et al., 2022) and the Atacama Cosmology Telescope (ACT, Aiola et
al., 2020) to deliver precise, sample variance dominated measurements of
cosmic microwave background (CMB) anisotropies down to arcminute scales in
intensity and polarization. Meanwhile, photometric and spectroscopic galaxy
surveys such as the Dark Energy Survey (Abbott et al., 2018), KiDS (de Jong et
al., ), BOSS (Dawson et al., 2013), DESI (Levi et al., 2013), Euclid (Laureijs
et al., 2011), the Vera Rubin Observatory (Abate et al., 2012), the Nancy
Grace Roman Space Telescope (Spergel et al., 2015), and SPHEREx (Crill et al.,
2020) are poised to accurately measure the positions and shapes of billions of
galaxies, reconstructing the large-scale structure (LSS) of the Universe out
to high redshifts. Fueled by a dramatic increase in the data throughput, these
experiments promise to map out a large portion of the Hubble volume to put
under stress the $\Lambda$CDM model, delivering groundbreaking insights into
the fundamental properties of the cosmos.

However, analyzing the large amounts of data generated by these experiments
poses significant challenges. In most cases, the analysis is carried out in
terms of summary statistics, such as the two-point function in either real or
harmonic space. At the heart of Bayesian inference pipelines, Einstein-
Boltzmann codes such as CAMB (Lewis et al., 2000) or CLASS (Blas et al., 2011)
solve for the evolution of cosmological perturbations and calculate CMB/LSS
observables, which are in turn compared to the data in the likelihood
evaluation. Depending on the dimensionality and the topology of the parameter
space, this operation is typically repeated $10^{4}$-$10^{6}$ times before
convergence thresholds are reached; considering that a call to a Boltzmann
solver might require from a second to a minute in case of high-precision
settings, these numbers show the high computational cost of Monte Carlo Markov
Chains (MCMC) analysis.

Emulation has emerged as a highly promising technique for accelerating
cosmological inference. It consists in developing surrogate models that
approximate the outputs of computationally expensive models at significantly
lower computational costs. This approach offers a twofold advantage: it
expedites calculations and enables the utilization of alternative optimization
and sampling methods, such as those based on gradients. In recent years,
emulators have proven successful in various domains with computationally
intensive forward models, including cosmology. Early contributions in this
area can be traced back to the work of Jimenez et al. (2004) and Fendt &
Wandelt (2006), which presented the initial applications of emulators based on
polynomial interpolation within the field of cosmology. The first emulator
based on neural networks (NN), called CosmoNet, was introduced in Auld et al.
(2007). Recent advancements include the work of CosmicNet (Albers et al.,
2019; Günther et al., 2022), with emulators targeting the replacement of the
most computationally demanding components within Boltzmann solvers.
Additionally, Nygaard et al. (2023), Günther (2023), and To et al. (2023)
adopted an iterative approach for emulator training to reduce the overall
computational resources required. Currently, the most precise CMB emulator
(i.e. matching the CAMB high-precision settings) that accurately coveres a
wide parameter space is CosmoPower (Spurio Mancini et al., 2022; Bolliet et
al., 2023), which directly emulates the CMB angular power spectra.

The cosmological community also developed emulators tailored for LSS analysis,
with emulators for the linear matter power spectrum (Mootoovaloo et al., 2022;
Spurio Mancini et al., 2022; Donald-McCann et al., 2022a; Aricò et al., 2021),
nonlinear corrections (Lawrence et al., 2017; DeRose et al., 2019; Knabenhans
et al., 2019, 2021; Angulo et al., 2021), power spectrum multipoles (Aricò et
al., 2021; Eggemeier et al., 2022; Donald-McCann et al., 2022b), galaxy survey
angular power spectrum (Manrique-Yus & Sellentin, 2020; Mootoovaloo et al.,
2020; Bonici et al., 2022; To et al., 2023) and even for the log-posterior
(Gammal et al., 2022). The amount of papers employing emulators is rapidly
growing, demonstrating an increasing interest of the cosmological community
toward this technique.

In this paper, we present a novel emulator of the CMB temperature,
polarization, and lensing potential anisotropies angular power spectra,
Capse.jl111https://github.com/CosmologicalEmulators/Capse.jl, which differs
from previously developed emulators in two fundamental ways. On the one hand,
we propose and implement a more efficient preprocessing, i.e. the operations
performed on the training data before being fed to the NN for the training,
and test a different technique for data compression. On the other hand, we
employ a specialized, performant NN framework in Julia, which results in our
emulator returning predictions in around 50 $\mu s$. Moreover, Capse.jl can be
trained in an hour’s time on a standard 8-cores CPU to cover a wide range of
cosmological parameters, and down to scales that will be mapped by upcoming
experiments. These characteristics make Capse.jl the CMB power spectra
emulator with the shortest evaluation time available as well as extremely
cheap to train. In addition, we also implement fully differentiable CMB
likelihoods. Coupled to our emulators, these methods enable the use of
gradient-based inference methods. The combination of these state-of-the-art
inference methods with Capse.jl allows us to analyze Planck data in less than
one second and ACT data in $\sim 10$ seconds.

This paper is structured as follows. In Sec. 2 we describe the architecture,
the preprocessing, and the training procedure of our emulators. In Sec. 3 we
describe the way we construct our likelihoods. In Sec. 4 we describe the
sampling algorithms employed in this work. In Sec. 5 we show the main results
of our work, namely the precision and the speed of our emulators coupled with
differentiable likelihoods. In Sec. 6 we compare our emulators with similar
works present in the literature. Finally, we conclude drawing the conclusions
of this study in Sec. 7.

##  2 Capse emulator

The core idea behind surrogate models is quite simple: replacing a
computationally expensive function with a sufficiently accurate approximation
that is cheaper. Although creating an emulator is conceptually
straightforward, there are several design choices involved in the process
which can have a drastic impact on the final performance. These choices
include the type and architecture of the surrogate model, the preprocessing
procedure, the application of dimensional reduction techniques, and
potentially incorporating domain knowledge.

To achieve optimal performance and train a fast and accurate emulator while
conserving computational resources, all of these choices become significant.
In this section, we provide a detailed description of our Capse.jl emulator,
focusing on the aspects that differentiate it from other emulators present in
literature. We conclude by discussing the accuracy of the trained emulators
using a validation data set that covers the entire range of emulation.

###  2.1 Architecture, preprocessing and training

The NN employed for our emulators is a multi-layer perceptron (MLP) with a
number of input features that corresponds to the number of cosmological
parameters considered, 5 hidden layers with 64 neurons each and a hyperbolic
tangent activation function. Our emulators are written in the Julia language
(Bezanson et al., 2015) and employ the NN implemented in SimpleChains.jl.
While frameworks like PyTorch, TensorFlow, and JAX offer general-purpose
capabilities to cover a wide range of use cases, from simple MLPs to complex
Large Language Models and from models running on a single CPU to those using
grids of TPUs, SimpleChains.jl has a way more limited focus: it excels at
running small NNs on CPUs. By narrowing down the scope, the library
prioritizes performance over generality, enabling several choices that improve
computational efficiency. This results in a significant boost compared to more
traditional frameworks, as demonstrated in previous comparisons with
SimpleChains.jl.222A more detailed description for the interested reader can
be found here.

A key aspect of our methodology is the preprocessing of training data, which
significantly impacts emulation performance. In this work, we perform min-max
normalization on each input and output feature. This process scales all
features to fall within a $[0,1]$ range, a technique known to enhance
emulation performance (Nygaard et al., 2023). Furthermore, we leverage domain
knowledge, exploiting the dominant linear dependence of the $C_{\ell}$’s on
the scalar amplitude $A_{s}$. To account for this, we divide each output
feature for the corresponding value of $A_{s}$, while still treating $A_{s}$
as an independent parameter within the NN. This approach analytically handles
the linear scaling, allowing the NN to focus on learning the more nonlinear
feature.333A similar rescaling for $\tau$ could potentially be employed.
However, as its effect varies across different power spectra and scales
(Alvarez et al., 2021) and since the desired accuracy was achieved with our
current procedure, further optimization for $\tau$ is left for future work. We
found that this enhances the emulation accuracy given the same amount of
training resources. This improvement can be easily explained by the fact that,
given the $\ln 10^{10}A_{s}$ emulation range, these preprocessing steps reduce
the range of the output features by almost a factor of 3, hence improving the
emulation performance. We checked it explicitly: the $A_{s}$ rescaling
enhances the emulator accuracy by a factor of $\approx 3$, almost the same as
the $C_{\ell}$’s dynamical range reduction.

The emulation space is sampled with a Latin-Hypercube Sampling (LHS) algorithm
(McKay et al., 1979) using 10,000 points. The ground-truth predictions have
been computed using the CAMB Boltzmann solver (Lewis et al., 2000) with the
same high-accuracy settings of Bolliet et al. (2023). In particular, we set
lens_potential_accuracy=8, lens_margin=2050, AccuracyBoost=2.0,
lSampleBoost=2.0, lAccuracyBoost=2.0 and DoLateRadTruncation=False. The
nonlinear corrections were computed using HMCode (Mead et al., 2020). The
emulation space is also chosen to exactly match that of Bolliet et al. (2023),
which allows us to make a fair comparison with their results; the details can
be found in Table 1. We trained a NN for each angular spectrum considered in
this work (namely $TT,\,TE,\,EE,\,\phi\phi$), using a standard laptop with an
8-core CPU using the ADAM optimizer (Kingma & Ba, 2017) with a mean square
error loss function for 50,000 epochs.

In this paper, we explore two different emulation approaches. In the first
approach, the output of the NN is a set of CMB power spectra $C_{\ell}$’s
evaluated on the multipole grid from $\ell_{\mathrm{min}}=2$ to
$\ell_{\mathrm{max}}=5000$ with a step of $\Delta\ell=1$. This emulator
delivers predictions for a single power spectrum in a meantime of $40\,\mu$s.
This improvement is truly remarkable when compared to the computation time of
predictions using a Boltzmann solver with high-accuracy settings, which
typically takes around $\mathcal{O}(60\,\text{s})$. Consequently, our
emulators achieve speeds that are approximately 1,000,000 times faster.

The second approach we consider involves a dimensionality reduction technique,
the Chebyshev polynomial decomposition. The key idea behind this approach is
to project the $\ell$-dependence of the $C_{\ell}$’s onto the Chebyshev
polynomial basis:

$C_{\ell}(\theta)\approx\sum_{n=0}^{N_{\mathrm{max}}}a_{n}(\theta)T_{n}(\ell),$
(1)

where $a_{n}$ represents the $n$-th coefficient that carries the dependence on
the cosmological parameters $\theta$, $T_{n}$ denotes the $n$-th Chebyshev
polynomial, and $N_{\mathrm{max}}$ is the highest degree of the polynomials
employed in the decomposition.444Note that we use a rescaled version of the
Chebyshev polynomials since they are defined on the interval
$\left[-1,1\right]$. The rescaling is automatically handled by the Chebyshev
projection library employed in our analysis, FastChebInterp.jl The emulator is
then trained to reproduce the cosmological dependence of the $a_{n}(\theta)$
coefficients. We choose to explore this approach for the nice mathematical
properties of Chebyshev interpolation; we refer the interested reader to
Trefethen (2012) and Wang (2023) for a detailed discussion.

The advantages of this polynomial approach are twofold. First, it
significantly speeds up the training process, since the number of output
features of the NN is reduced. We have found that using a polynomial order of
$N_{\mathrm{max}}=48$ provides accurate emulation while simplifying the
complexity of the neural network. This reduction in complexity not only lowers
the training cost but also improves the emulator speed, enabling predictions
to be generated in approximately $15\,\mu$s. This represents a 2.5-fold
increase in speed compared to the previous scenario.

The second, perhaps less obvious, advantage stems from this simple
observation. By performing a linear decomposition of the $C_{\ell}$’s and
noting that the basis functions do not depend on the cosmological parameters,
we can leverage this linear decomposition in the likelihood computation. If
the operations needed to perform the likelihood computation can be represented
as a linear operator $\mathds{L}$ acting on the theory vector, we can utilize
the linear decomposition as follows:

$\mathds{L}\sum_{n}a_{n}(\theta)T_{n}=\sum_{n}a_{n}(\theta)\mathds{L}T_{n}\equiv\sum_{n}a_{n}\hat{T}_{n},$
(2)

where we have defined:

$\hat{T}_{n}\equiv\mathds{L}T_{n}.$ (3)

If the linear operator $\mathds{L}$ does not depend on the cosmological
parameters, as it typically is the case, we can compute the action of
$\mathds{L}$ on the basis function once and store the results. As we will
demonstrate later, this approach leads to significant speed improvements in
the analysis.

Parameter | Min | Max
---|---|---
$\ln 10A_{\mathrm{s}}$ | 2.5 | 3.5
$n_{\mathrm{s}}$ | 0.88 | 1.06
$H_{0}[\mathrm{km}/s/\mathrm{Mpc}]$ | 40 | 100
$\omega_{\mathrm{b}}$ | 0.01933 | 0.02533
$\omega_{\mathrm{c}}$ | 0.08 | 0.20
$\tau$ | 0.02 | 0.12
Table 1: The parameter ranges used to generate the Latin hypercube of
cosmological parameters for computing the training data. It is important to
note that our emulators should only be used within these specified ranges.

###  2.2 Accuracy tests

We carry out two types of tests to assess the performance of the emulators.
The first check involves comparing the emulator’s predictions with the
corresponding ground truth values. This comparison allows us to measure the
discrepancy between the emulator and the Boltzmann solver output. In this
section, we will thoroughly discuss this validation approach. The second check
involves running a Markov Chain Monte Carlo (MCMC) analysis using specific
datasets. The MCMC chains produced by the emulator are then compared with the
chains obtained from the Boltzmann code. This comparison serves as an
additional validation step to assess the emulator’s accuracy and reliability.
The details and results of this MCMC analysis will be presented in Sec. 5.

We use CAMB to generate a validation dataset consisting of 20,000 combinations
of the cosmological parameters drawn from the same emulation range of the
training dataset (see Tab. 1). We emphasize that these specific combinations
are not present in the training dataset. Our objective is to assess the
accuracy of the emulators by comparing their predictions to the ground truth
values. To quantify the emulation error, we calculate the difference between
the emulated $C_{\ell}^{XY}$ and the true $C_{\ell}^{XY}$, normalized by the
corresponding spectrum variance $\sigma_{\ell}^{XY}$

$\frac{\left|C_{\ell,\text{ emulated }}^{XY}-C_{\ell,\text{ true
}}^{XY}\right|}{\sigma_{\ell}^{XY}},$ (4)

where $XY\in[TT,\,TE,\,EE,\,\phi\phi]$ represent the CMB temperature,
polarization, and lensing potential anisotropies power spectra.

To perform this test, we need both the ground truth spectra for computing the
emulation error and the expected spectra variance $\sigma_{\ell}$. The
variance can be computed, under the assumption that both
$\left[X,\,Y\right]\in\left[{T,\,E,\,\phi}\right]$ behave as Gaussian random
fields, as (e.g., Knox, 1995; Kamionkowski et al., 1997):

$\displaystyle\sigma_{\ell,\mathrm{CMB}}^{XY}=$
$\displaystyle\sqrt{\frac{1}{f_{\mathrm{sky}}(2\ell+1)}}$ (5)
$\displaystyle\times\sqrt{(C_{\ell}^{XY})^{2}+\left(C_{\ell}^{XX}+N_{\ell}^{XX}\right)\left(C_{\ell}^{YY}+N_{\ell}^{YY}\right)},$

where $f_{\mathrm{sky}}$ and $N_{\ell}^{XX}$ are respectively the observed sky
fraction and the noise power spectrum of field $X$ for the experimental
configuration under consideration. We consider two representative upcoming
ground-based CMB surveys, Simons Observatory555Available at
https://github.com/simonsobs/so_noise_models. For temperature and polarization
we use the following LAT_comp_sep_noise/v3.1.0 curves:
SO_LAT_Nell_T_atmv1_goal_fsky0p4_ILC_CMB.txt and
SO_LAT_Nell_P_baseline_fsky0p4_ILC_CMB_E.txt. For lensing we choose the
nlkk_v3_1_0_deproj0_SENS2_fsky0p4_
it_lT30-3000_lP30-5000.dat in the LAT_lensing_noise/lensing_v3_1_1/ folder.
(Ade et al., 2019) and CMB-S4666The primary CMB noise curves are available at
https://sns.ias.edu/~jch/S4_190604d_2LAT_Tpol_default_noisecurves.tgz,
specifically we use
S4_190604d_2LAT_T_default_noisecurves_deproj0_SENS0_mask_16000_ell_TT_yy.txt
and
S4_190604d_2LAT_pol_default_noisecurves_deproj0_SENS0_mask_16000_ell_EE_BB.txt
for $TT$ and $EE$ respectively. For the CMB lensing reconstruction noise we
use the kappa_deproj0_sens0_16000_lT30-3000_lP30-5000.dat curves from
https://github.com/toshiyan/cmblensplus/tree/master/example/data. (Abazajian
et al., 2019). We refer the reader to Bolliet et al. (2023) and references
therein for details regarding the surveys sensitivities. In both cases, we use
projected noise curves obtained after foreground cleaning and including multi-
frequency information from the Planck satellite. For both experiments we
assume an observed sky fraction of $f_{\rm sky}=0.4$; the results can be seen
in Figs. 1 and 5.

The degraded EE 2-point function performance in the $\ell<20$ region arises
because of two different factors. First and foremost, we only consider cosmic
variance over this multipole range since S4 will not be able to image the
largest angular scales at $\ell\lesssim 40$ and therefore they do not provide
noise curves on this range. The second effect comes from the scaling of the
2-point function with respect to $\tau$, which is proportional to $\tau^{2}$
for $\ell<30$. In fact, the dynamical range of the $C_{\ell}^{EE}$ in this
specific multipole range, across all the cosmologies considered, is 30 times
bigger than the one outside. We expect that, when including the $\tau$
rescaling in a future work, we will be able to improve the low-$\ell$ accuracy
as well. If we exclude the $C_{\ell}$’s with $\ell<13$, the results are even
better, with an emulation error that is below $0.08\,\sigma$ for $99\%$ of the
points in the validation dataset.

Figure 1: Comparison of the relative differences between Capse.jl and high-
accuracy predictions from CAMB for the four 2-point statistics emulated in
this study. The differences are normalized to the forecasted CMB-S4
statistical uncertainties ($f_{\mathrm{sky}}=0.4$) and measured in units of
cosmic variance for $\ell<40$. The shaded regions represent the 68%, 95%, and
99% percentiles of the respective distributions.

##  3 Cosmological inference framework

After providing a detailed description of the Capse.jl emulator, our focus
shifts to leveraging its CMB power spectra predictions for cosmological
inference. Specifically, we are interested in drawing samples from the
posterior distribution of the model parameters $\boldsymbol{\theta}$ given the
data $\boldsymbol{d}$, $P(\boldsymbol{\theta}|\boldsymbol{d})$, after
specifying a likelihood function
$\mathcal{L}(\boldsymbol{d}|\boldsymbol{\theta})$ and the priors
$\Pi(\boldsymbol{\theta})$. To this end, we use the Julia Bayesian statistical
inference library Turing.jl (Ge et al., 2018).

Planck | ACT
---|---
Parameter | Prior | Parameter | Prior
$\ln 10A_{\rm{s}}$ | $U(2.5,3.5)$ | $\ln 10A_{\rm{s}}$ | $U(2.8,3.2)$
$n_{\rm{s}}$ | $U(0.88,1.06)$ | $n_{\rm{s}}$ | $U(0.88,1.06)$
$H_{\rm{0}}$ | $U(40,100)$ | $H_{\rm{0}}$ | $U(40,100)$
$\omega_{\rm{b}}$ | $U(0.01933,0.02533)$ | $\omega_{\rm{b}}$ | $U(0.01933,0.02533)$
$\omega_{\rm{c}}$ | $U(0.08,0.20)$ | $\omega_{\rm{c}}$ | $U(0.08,0.20)$
$\tau$ | $\mathcal{N}(0.0506,0.0086)$ | $\tau$ | $\mathcal{N}(0.065,0.015)$
$A_{\mathrm{Planck}}$ | $\mathcal{N}(1.0,0.001)$ | $y_{\rm{p}}^{2}$ | $U(0.7,1.3)$
Table 2: Priors on the cosmological and calibration parameters varied in the
analyses of Planck and ACT datasets. $U(a,b)$ denotes a uniform distribution
between $[a,b]$, while $\mathcal{N}(\mu,\sigma)$ indicates a Gaussian
distribution with mean $\mu$ and variance $\sigma^{2}$.

Turing.jl is a Probabilistic Programming Language (PPL). Unlike traditional
programming languages, which are designed for deterministic computations, PPLs
are designed to explicitly incorporate probabilistic constructs. This allows
the user to explicitly declare the priors for the input parameters and write
their relationship to the observations within a statistical mode. Moreover,
the user can also declare the distribution of the observations. Then,
Turing.jl uses this information to automatically draw an expression for the
likelihood by conditioning the model on the observations.

Below is a Julia code snippet demonstrating how a likelihood is implemented
using Turing.jl:

⬇

@model function CMB_planck(data, covariance)

#define the priors

ln10As ~ Uniform(2.5, 3.5)

ns ~ Uniform(0.88, 1.06)

H0 ~ Uniform(40, 100)

omb ~ Uniform(0.01985, 0.025)

omc ~ Uniform(0.08, 0.20)

tau ~ Normal(0.0506, 0.0086)

Aplanck ~ Normal(1.0, 0.0025)

#create a parameters vector

params = [ln10As, ns, H0, omb, omc, tau]

#perform the theoretical computation

pred = theory_planck(params) ./ (Aplanck^2)

#compute the likelihood

data ~ MvNormal(pred, covariance)

return nothing

end

This is slightly different from the actual likelihood being used, which
performs a reparametrization as explained in Appendix A.

In this work, we use the distributions listed in Tab. 2 as priors for the
input parameters of the Capse.jl emulator. The only exception is $\tau$, for
which we employ a Gaussian prior, stemming from the analysis of the low-$\ell$
$EE$ angular spectrum. Moreover, since both Planck, ACT, and SPT-3G
measurements are Gaussian-distributed (see Sec. 5.1), we use Turing.jl to draw
a Gaussian likelihood of the shape:

$\boldsymbol{d}\sim\mathcal{N}(\boldsymbol{t}(\boldsymbol{\theta}),\boldsymbol{C})\,,$
(6)

where $\boldsymbol{t}(\boldsymbol{\theta})$ is a theory vector computed as a
function of the model parameters $\boldsymbol{\theta}$ and $\boldsymbol{C}$ is
the covariance matrix of the data vector $\boldsymbol{d}$.

While the bottleneck is typically the computation of the theory vector
$\boldsymbol{t}(\boldsymbol{\theta})$, even for the other emulators described
in the literature, this is not true anymore in our case, especially for the
data and theory vectors of the Planck and ACT analyses. The binning of the
$C_{\ell}$’s and the $\chi^{2}$ computation are more expensive by a factor of
7 than our emulators computations.

In order to overcome this issue, we follow two strategies. The first is the
one depicted in Eq. 2. As anticipated in Sect. 2.1, all the operations
required to compute the Planck and ACT likelihoods can be represented as a
linear operator $\mathds{L}$ on $\boldsymbol{t}(\boldsymbol{\theta})$. In this
case, we can compute its action on the basis function once
$\hat{T}_{n}=\mathds{L}T_{n}$ and store it. During the MCMC, we multiply
$\hat{T}_{n}$ with $a_{n}(\boldsymbol{\theta})$ and then compute the
$\chi^{2}$. This procedure, while speeding up the analysis, does not introduce
any theoretical error: the computation of the binned $C_{\ell}$’s using the
two approaches matches up to machine-precision.

A further optimization that we can introduce is using MOPED, a linear data
compression algorithm (Heavens et al., 2000; Heavens et al., 2017). Given a
$n$-dimensional data vector whose theoretical prediction depends on $m$
parameters, MOPED generates $n$ $m$-dimensional data vector,
$\boldsymbol{b}_{i}$. Moreover, if the noise of the data does not depend on
the parameters of the model, as in the case of Planck and ACT measurements,
the original and compressed data vectors can be shown to have identical Fisher
information matrices with respect the parameters of the model (Heavens et al.,
2000). The MOPED compression vectors $\boldsymbol{b}_{i}$ can be computed
using:

$\boldsymbol{b}_{i}=\frac{\boldsymbol{C}^{-1}\boldsymbol{t}_{,i}-\sum_{j=1}^{i-1}\left(\boldsymbol{t}_{,i}^{T}\boldsymbol{b}_{j}\right)\boldsymbol{b}_{j}}{\sqrt{F_{ii}-\sum_{j=1}^{i-1}\left(\boldsymbol{t}_{,i}^{T}\boldsymbol{b}_{j}\right)^{2}}},$
(7)

where $\boldsymbol{t}_{,i}$ is the derivative of the theory vector with
respect to the $i$-th parameter and $F_{ii}$ is the $i$-th element on the
diagonal of the Fisher matrix $\boldsymbol{F}$:

$F_{ij}=-\frac{\partial^{2}\mathcal{L}(\theta)}{\partial\theta_{i}\partial\theta_{j}}.$
(8)

As stated in Campagne et al. (2023), it is pretty straightforward to evaluate
both the Fisher matrix and the $\boldsymbol{t}$ derivatives using automatic
differentiation (AD).

Therefore, we use MOPED to perform a compression of the likelihood shown in
Eq. 6. This allows us to transform a relatively expensive $n$-dimensional
Gaussian likelihood into the product of $m$ one dimensional unit-variance
Gaussian likelihoods:

$\mathcal{L}(\boldsymbol{d}|\boldsymbol{\theta})\approx\Pi^{m}_{i}\mathcal{N}(\boldsymbol{t}(\boldsymbol{\theta})M_{i},1)\,,$
(9)

where $M_{i}$ is the MOPED compression matrix for the $i$-th parameter of the
model. However, unlike our first strategy, using the MOPED compression changes
the actual likelihood being targeted, opening the possibility for changes the
results of our analysis.

##  4 Samplers

We now outline the methods employed to explore the likelihoods introduced in
Sec. 3 and derive parameter constraints for the Planck, ACT, and SPT-3G
datasets. The algorithms employed, all based on leveraging the likelihood
gradient, are:

  * •

NUTS, a version of the Hamiltonian Monte Carlo sampler (Hoffman & Gelman,
2011). This is considered the state-of-the-art algorithm to be employed in
MCMC analysis.

  * •

MicroCanonical Hamiltonian Monte Carlo, a recently proposed extension of
Hamiltonian MonteCarlo that promises to be competitive with NUTS (Robnik et
al., 2022; Robnik & Seljak, 2023).

  * •

Pathfinder, a variational inference methods (Zhang et al., 2022). Although
less robust than the previous methods, as it lacks of asymptotic convergence
properties, Pathfinder excels at quickly delivering draws from the typical
set, that can be representative of the posteriors in the case of Gaussian or
mildly non-Gaussian posterior, or used to initialize standard MCMC.

In this section we describe all these methods, while addressing the interested
reader to the relevant works in the literature.

###  4.1 Hamiltonian MonteCarlo

Hamiltonian Monte Carlo (HMC) (Betancourt, 2018) is a MCMC algorithm that
explores a parameter space by simulating the dynamics of a Hamiltonian system.
The fundamental idea behind HMC is to introduce auxiliary momentum variables
that are independent of the target function. The joint distribution of the
position (i.e. the original parameters) and momentum variables is then defined
by a Hamiltonian function, which governs the dynamics of the system. The
Hamiltonian is typically chosen to be the sum of the potential energy, defined
as the negative logarithm of probability density of the likelihood, and the
canonical kinetic energy of the momentum variables.

In each HMC iteration, a proposal is generated by simulating the dynamics of
the system for a fixed number of steps using a numerical integration scheme.
The acceptance probability of the proposal is then computed using the
Metropolis-Hastings (MH) algorithm based on the Hamiltonian energy of the
sample. Then, it can be shown that the target distribution can be returned by
marginalizing the momenta variables:

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

One of the challenges of HMC is choosing the step size and number of
integration steps, which can have a significant impact on the performance of
the algorithm. To address this challenge, the No-U-Turn Sampler (NUTS)
(Hoffman & Gelman, 2011) algorithm was proposed as an extension of HMC. NUTS
addresses this issue by introducing a recursive algorithm that determines the
optimal number of steps and step size for each iteration. The algorithm
generates a ”tree” of proposals by evolving the Hamiltonian dynamics forwards
and backwards in time. The number of steps is then chosen such that the
trajectory does not turn on itself. This prevents the sampler from returning
to previously explored regions of the parameter space. The algorithm then
evaluates the Hamiltonian at each new proposal and computes an acceptance
probability using MH. The step size is chosen such that a given target
acceptance probability is met.

###  4.2 MicroCanonical HMC

Recently, Robnik et al. (2022) developed a new family of HMC sampling
algorithms based on the Micro-Canonical ensemble, named Micro-Canonical
Hamiltonian Monte Carlo (MCHMC). The key feature of MCHMC is that it only
makes use of a single Hamiltonian energy level to explore the whole parameter
space; i.e. the Hamiltonian energy of the system is conserved. Mathematically
this means that the target distribution is obtained as:

$\mathcal{P}(\boldsymbol{x})\propto\int\delta(H(\boldsymbol{x},\boldsymbol{p})-E)\,d\boldsymbol{p}\,.$
(11)

As shown in Robnik et al. (2022), there is a multitude of Hamiltonians that
satisfy Eq. 11. In this work, we use a Julia implementation of the variable
mass Hamiltonian as presented in Robnik et al. (2022).

MCHMC poses several advantages over traditional HMC. First, it does not
require the Metropolis adjustment to return unbiased posteriors as in
traditional HMC. Second, the bias on the estimated distributions depends on
the Hamiltonian error per dimension while in traditional HMC it depends on the
absolute Hamiltonian error. Third, the dynamics of the variable mass
Hamiltonian result in the particle moving the slowest in the regions where the
likelihood density is the highest. This is the opposite behaviour to
traditional HMC where the particle moves the fastest where the likelihood
density is the highest. In combination, these three features allow MCHMC to
explore high dimensional parameter spaces more efficiently than traditional
HMC specially for very large numbers of dimensions, as shown in Bayer et al.
(2023).

###  4.3 Pathfinder

Sometimes exploring the parameter space can become too expensive even when the
gradient of the likelihood is known. For this reason, the scientific community
has developed a series of posterior approximation methods.

The most popular of such methods is variational inference (VI) (Ranganath et
al., 2013; Kucukelbir et al., 2016). The goal of VI is to find a tractable
approximate distribution to the posterior. In order to do so, VI uses a family
of distributions (known as the variational family) in terms of the original
parameters of the problem and a series of variational parameters. VI then
finds the variational parameters that minimize the Kullback-Leibler (KL)
divergence with respect to the posterior. This results in faithful
approximations to the posterior even if the actual posterior falls out of the
variational family.

Pathfinder is a VI algorithm that locates approximations to the target density
along a minimization path (Zhang et al., 2022). Pathfinder uses L-BFGS, a
quasi-Newton method that relies on gradients to estimate the Hessian matrix.
While moving along the minimization path, samples are drawn from a
multivariate normal distribution characterized by the estimated Hessian
matrix. Pathfinder outperforms other VI methods in providing quick estimates
of the posterior, as shown in Zhang et al. (2022) for a series of benchmark
problems.

Moreover, Pathfinder can be run in parallel, greatly improving its estimates
of the posterior by combining the approximations of each independent run using
Pareto-Smoothed Importance Sampling (PSIS) at a low computational cost
(Vehtari et al., 2022). Furthermore, PSIS provides a quantitative metric which
can be used to assess the quality of the posterior found by Pathfinder.

Last but not least, Pathfinder can be used to accelerate MCMC methods, using
its draws to start chains, whose initial points are usually drawn from the
prior and hence are far from the typical set. Since Pathfinder draws are (very
likely to be) on the typical set, this reduces the length of the so called
burn-in phase, where chains approach the high mass probability region.777It is
possible to further reduce the burn-in phase using the estimated covariance
from Pathfinder as the mass matrix provided to the HMC sampler, but we did not
take advantage of this feature in this work.

##  5 Results

In this section we discuss the main results of this work.

We start by testing the accuracy of our emulators performing the MCMC analysis
of three different datasets, the Planck (Aghanim et al., 2020b), ACT DR4
(Aiola et al., 2020), and 2018 SPT-3G (Balkenhol et al., 2023) datasets. The
analyses are performed using both CAMB and Capse.jl and the obtained
posteriors are compared in order to assess the emulator accuracy.

###  5.1 Datasets

####  5.1.1 Planck (2018)

In this subsection, we show the constraints obtained using data from the
Planck satellite. Specifically, we consider the compressed temperature and
polarization Planck 2018 Plik_lite likelihood, which is constructed from the
high-$\ell$ $TT+TE+EE$ CMB bandpowers and covariances that have been
premarginalized over astrophysical contributions in the spectrum estimation
step (Aghanim et al., 2020a). The only nuisance parameter left in this
compressed likelihood is the overall Planck calibration factor
$A_{\mathrm{Planck}}$ which simply rescales the foreground-marginalized
bandpowers as $\hat{C}_{\ell}\to\hat{C}_{\ell}/A_{\mathrm{Planck}}^{2}$. This
likelihood is available in Cobaya as planck_2018_highl_plik.TTTEEE_lite.

In Fig. 2, we compare the contours obtained with CAMB and Capse.jl. The former
has been used within Cobaya and a nuisance-marginalized Planck likelihood,
with a termination criterion of $|\hat{R}-1|<0.01$, requiring 6 hours to
converge with 4 parallel chains, each of them with 4 threads, for a total of
around 100 CPU-hours. The latter employed
PlanckLite.jl,888https://github.com/JuliaCosmologicalLikelihoods/PlanckLite.jl
a pure Julia implementation of the Planck lite likelihood, and used 6 NUTS
chains, with 500 burn-in and 1,000 accepted steps, on a laptop in around 4
minutes of wall-clock time using less than $0.5$ CPU-hours, but still reaching
an $\hat{R}-1$ smaller than 0.002 for all parameters. We emphasize that the
NUTS analyses was performed with only a handful of parameters to tune: the
number of parallel chains, the burn-in and sampling steps per chains, and the
target acceptance ratio. We did not provide any input such as a covariance
matrix or a proposal step-size for each parameter. This shows the robustness
of the NUTS algorithm, which is able to perform a Bayesian analysis with
little to no input from the user. This is very different from the MH sampler
implemented in Cobaya, which greatly benefits in its performance from being
finely tuned by the user.

Regarding the accuracy of the analysis, the two contours are very similar,
with differences on the marginalized 1D posteriors smaller than $0.1\,\sigma$,
as is shown in the lower half of Fig. 2.

Figure 2: Triangle plot, showing the standard Planck lite chains and the one
obtained using Capse.jl in combination with Turing.jl. The lower triangular
part of the plot focuses on the comparison between CAMB and Capse.jl, in order
to ensure the precision of our emulators and likelihoods. The upper half of
the plot shows only Capse.jl contours, with a comparison of the results
between the samplers employed in our analysis, namely NUTS, MCHMC, and
Pathfinder.

####  5.1.2 ACT DR4

Here we show the results from the chains based on the ACT Data Release 4 (DR4)
likelihood (Aiola et al.,
2020).999https://github.com/ACTCollaboration/pyactlike This likelihood is
constructed from the temperature and polarization $TT/TE/EE$ power spectra
measured with ACT data taken between 2013 and 2016 covering more than 17,000
deg2 of the sky with spatially varying depth. Specifically, the ACT likelihood
relies on cleaned CMB bandpowers in the multipole range
$300\lesssim\ell\lesssim 4000$ that have been marginalized over foreground
emissions including the Sunyaev-Zel’dovich effects, diffuse Galactic emission,
and the cosmic infrared background. The covariance of the CMB power spectra
already includes contributions from foreground, calibration, and beam
uncertainties. The only nuisance parameter that is allowed to freely vary is
the overall polarization efficiency $y_{p}$, which linearly and quadratically
rescales the $C_{\ell}^{TE}$ and $C_{\ell}^{EE}$ spectra, respectively. Also
in this case we have coded a pure Julia version of this likelihood,
ACTPolLite.jl.101010https://github.com/JuliaCosmologicalLikelihoods/ACTPolLite.jl
For the CAMB$+$Cobaya chains we used the same hardware settings from before.
However, since we used the high-precision settings for CAMB, the computational
resources employed are higher than before, being around 480 CPU-hours required
to reach convergence. Regarding our Capse.jl$+$Turing.jl, we run as before 6
chains in parallel with 500 burn-in steps and 1,000 accepted steps, obtaining
convergent chains with $|\hat{R}-1|<0.002$ using around one CPU-hour. Also in
this case the agreement between the CAMB and the Capse.jl chains is excellent,
with differences on the posteriors smaller than $0.1\,\sigma$, as shown in the
lower half of Fig. 3.

Figure 3: Triangle plot, showing the standard ACTPolLite chains and the one
obtained using Capse.jl in combination with Turing.jl. The lower triangular
part of the plot focuses on the comparison between CAMB and Capse.jl, in order
to ensure the precision of our emulators and likelihoods. The upper half of
the plot shows only Capse.jl contours, with a comparison of the results
between the samplers employed in our analysis, namely NUTS, MCHMC, and
Pathfinder.

####  5.1.3 SPT-3G

The final CMB likelihood considered in this study pertains to the
multifrequency 2018 SPT-3G temperature and polarization dataset, as detailed
in Balkenhol et al.
(2023).111111https://pole.uchicago.edu/public/data/balkenhol22/ This dataset
comprises measurements of the CMB $TT/TE/EE$ power spectra across
$\sim\,$1,500 deg2 of the southern Sky. Specifically, the $TT$ power spectra
span the angular multipole range from $750<\ell\leq 3,000$,while the $TE$ and
$EE$ spectra extend from $300<\ell\leq 3,000$. Distinct from the
PlancklLite.jl and ACTPolLite.jl likelihoods, which are premarginalized over
astrophysical emissions, the SPT-3G one employs a parametric model for
foreground characterization in both temperature and polarization.
Additionally, the SPT-3G likelihood accounts for the effects of instrumental
calibration and beam uncertainties, aberration due to the relative motion with
respect to the CMB rest frame, and super-sample lensing (Dutcher et al.,
2021). Owing to its large number of foreground and nuisance parameters, the
SPT-3G likelihood serves as an exemplary testbed for gradient-based sampling
techniques. More precisely, in addition to the 6 $\Lambda$CDM parameters, we
sample over 6 calibration parameters, 26 foreground parameters, and a nuisance
parameter (super-sample lensing convergence), for a total of 39 parameters. A
thorough description of the data model can be found in Balkenhol et al.
(2023). A pure Julia implementation of the SPT-3G likelihood,
SPTLikelihoods.jl, is made available in this case
too.121212https://github.com/JuliaCosmologicalLikelihoods/SPTLikelihoods.jl
For the CAMB$+$Cobaya chains we used the same hardware settings from before.
However, since we used the high-precision settings for CAMB, the computational
resources employed are higher than before, being around 1,600 CPU-hours
required to reach convergence. We also run 6 chains, with 500 burn-in steps
and 1,000 accepted steps, using Capse.jl$+$Turing.jl: the obtained chains have
a convergence with $|\hat{R}-1|<0.002$ using 14 CPU hours. Also with this more
challenging posterior there is an excellent agreement between the standard and
the emulator chains, with differences on the posterior smaller than
$0.1\,\sigma$, as shown in Fig. 4.

Figure 4: Triangle plot, showing the standard SPT-3G chains and the one
obtained using Capse.jl in combination with Turing.jl. The lower triangular
part of the plot focuses on the comparison between CAMB and Capse.jl, in order
to ensure the precision of our emulators and likelihoods. The upper half of
the plot shows only Capse.jl contours, with a comparison of the results
between the samplers employed in our analysis, namely NUTS and MCHMC.

###  5.2 Samplers comparison

|  CAMB+Cobaya |  Capse.jl+NUTS |  Capse.jl+MCHMC |  Capse.jl+Pathfinder
---|---|---|---|---
Planck, Parameter and 95% limits
$\ln 10^{10}A_{s}$ | $3.051^{+0.035}_{-0.035}$ | $3.051^{+0.033}_{-0.033}$ | $3.051^{+0.033}_{-0.033}$ | $3.050^{+0.031}_{-0.031}$
$n_{s}$ | $0.9626^{+0.0089}_{-0.0087}$ | $0.9630^{+0.0085}_{-0.0087}$ | $0.9630^{+0.0085}_{-0.0084}$ | $0.9628^{+0.0092}_{-0.0090}$
$H_{0}[\mathrm{km}/s/\mathrm{Mpc}]$ | $67.0^{+1.2}_{-1.2}$ | $67.0^{+1.2}_{-1.2}$ | $67.1^{+1.2}_{-1.2}$ | $67.1^{+1.2}_{-1.3}$
$\omega_{b}$ | $0.02233^{+0.00029}_{-0.00029}$ | $0.02234^{+0.00029}_{-0.00029}$ | $0.02234^{+0.00029}_{-0.00028}$ | $0.02234^{+0.00029}_{-0.00030}$
$\omega_{c}$ | $0.1208^{+0.0028}_{-0.0028}$ | $0.1207^{+0.0027}_{-0.0027}$ | $0.1207^{+0.0027}_{-0.0026}$ | $0.1207^{+0.0029}_{-0.0028}$
$\tau$ | $0.057^{+0.017}_{-0.016}$ | $0.057^{+0.016}_{-0.016}$ | $0.057^{+0.016}_{-0.016}$ | $0.056^{+0.015}_{-0.015}$
$A_{\textrm{Planck}}$ | $1.0004^{+0.0047}_{-0.0048}$ | $1.0005^{+0.0049}_{-0.0049}$ | $1.0005^{+0.0048}_{-0.0047}$ | $1.0006^{+0.0047}_{-0.0046}$
ACT, Parameter means and 95% limits
$\ln 10^{10}A_{s}$ | $3.047^{+0.056}_{-0.057}$ | $3.048^{+0.060}_{-0.059}$ | $3.047^{+0.058}_{-0.059}$ | $3.045^{+0.053}_{-0.050}$
$n_{s}$ | $1.011^{+0.029}_{-0.028}$ | $1.010^{+0.031}_{-0.030}$ | $1.010^{+0.030}_{-0.029}$ | $1.013^{+0.031}_{-0.030}$
$H_{0}[\mathrm{km}/s/\mathrm{Mpc}]$ | $68.4^{+3.0}_{-2.9}$ | $68.3^{+3.0}_{-2.9}$ | $68.3^{+2.9}_{-2.9}$ | $68.3^{+2.7}_{-2.6}$
$\omega_{b}$ | $0.02155^{+0.00060}_{-0.00058}$ | $0.02155^{+0.00061}_{-0.00060}$ | $0.02155^{+0.00059}_{-0.00059}$ | $0.02150^{+0.00061}_{-0.00061}$
$\omega_{c}$ | $0.1167^{+0.0071}_{-0.0073}$ | $0.1168^{+0.0073}_{-0.0072}$ | $0.1168^{+0.0072}_{-0.0070}$ | $0.1167^{+0.0066}_{-0.0063}$
$\tau$ | $0.064^{+0.026}_{-0.027}$ | $0.064^{+0.029}_{-0.028}$ | $0.064^{+0.028}_{-0.028}$ | $0.064^{+0.025}_{-0.025}$
$y_{p}$ | $1.0006^{+0.0095}_{-0.0084}$ | $1.0005^{+0.0090}_{-0.0091}$ | $1.0005^{+0.0089}_{-0.0090}$ | $1.0004^{+0.0091}_{-0.0091}$
SPT, Parameter means and 95% limits
$\ln 10^{10}A_{s}$ | $3.034^{+0.041}_{-0.041}$ | $3.034^{+0.039}_{-0.038}$ | $3.034^{+0.040}_{-0.041}$ | $-$
$n_{s}$ | $0.971^{+0.032}_{-0.033}$ | $0.972^{+0.030}_{-0.031}$ | $0.972^{+0.033}_{-0.032}$ | $-$
$H_{0}[\mathrm{km}/s/\mathrm{Mpc}]$ | $68.2^{+3.1}_{-2.9}$ | $68.2^{+3.0}_{-2.9}$ | $68.2^{+3.0}_{-2.9}$ | $-$
$\omega_{b}$ | $0.02225^{+0.00061}_{-0.00062}$ | $0.02225^{+0.00061}_{-0.00059}$ | $0.02224^{+0.00062}_{-0.00062}$ | $-$
$\omega_{c}$ | $0.1168^{+0.0074}_{-0.0074}$ | $0.1168^{+0.0073}_{-0.0073}$ | $0.1168^{+0.0075}_{-0.0074}$ | $-$
$\tau$ | $0.052^{+0.015}_{-0.014}$ | $0.053^{+0.015}_{-0.015}$ | $0.052^{+0.014}_{-0.015}$ | $-$
Table 3: Mean and 95% limits of the 1D marginal posteriors of the cosmological parameters obtained from the Planck, ACT, and SPT analyses, respectively. Different columns show the constraints obtained performing the same analyses using the sampling algorithms MH, NUTS, MCHMC and Pathfinder. The mean and limits have been computed using GetDist. Given the highly unstable results obtained when applying Pathfinder to the SPT posterior, we did not quote any number for this specific application. | Planck | ACT | SPT
---|---|---|---
| ESS/s | Comp. (s) | ESS/s | Comp. (s) | ESS/s | Comp. (s)
CAMB+Cobaya | 0.004 | 345,600 | 0.0003 | 1,728,000 | 0.0003 | 5,760,000
Capse.jl+NUTS | 1.597 | 7,014 | 1.914 | 11,813 | 0.09 | 50,400
Capse.jl+MCHMC | 3.686 | 4,374 | 4.136 | 7,792 | 0.335 | 43,200
Capse.jl+Pathfinder | N/A | 15 | N/A | 16 | - | -
Table 4: Shows the ESS per second and computation time achieved by the Cobaya,
NUTS, MCHMC and Pathfinder inference algorithms when performing the Planck
(2nd and 3rd columns), ACT (4th and 5th columns), and SPT (6th and 7th column)
analyses. Note that Cobaya used the cosmological code CAMB while the other
samplers were paired with the Capse.jl emulator. Given the highly unstable
results obtained when applying Pathfinder to the SPT posterior, we did not
quote any number for this specific application.

After validating our emulators, we now focus on the different inference
algorithms and study their impact on the speed, accuracy, and robustness of
the analyses.

The compatibility of the Capse.jl emulating framework with AD methods makes it
extremely synergistic with a plethora of gradient-based samplers. Here we
compare the performance of three gradient-based samplers (see Sec. 4 for
details). First, we compare the state-of-the-art NUTS sampler used to obtain
our fiducial results in Sec. 5.1 with the recently developed MCHMC sampler.
Second, we compare NUTS against a VI method, Pathfinder131313We use the
publicly available Julia implementations of NUTS, MicroCanonicalHMC.jl and
Pathfinder.jl.. In all cases we use the same likelihoods described in Sec. 5.1
to draw the comparisons. This is possible thanks to the modularity of the
Julia ecosystem, which allows us to write a likelihood in Turing.jl and then
sample it using different libraries and algorithms.

####  5.2.1 NUTS vs MCHMC

Before comparing the performance of the two samplers, it is essential to
verify that NUTS and MCHMC return the same posterior distributions. As MCMC
methods, both NUTS and MCHMC should asymptotically converge to the same target
distribution. The comparison of the NUTS and MCHMC posteriors for Planck, ACT,
and SPT analyses can be found in Figs. 2, respectively. The corresponding
numerical values for the quoted constraints can be found in Tab. 3. In this
table, we observe that the constraints obtained from both NUTS and MCHMC are
within a 0.1$\sigma$ distance of their Metropolis-Hastings (MH) counterpart.
The excellent agreement between the two samplers indicates that differences
are well within the sampling noise.

Having established that both samplers converge to the same posterior, the
pivotal question now is how quickly they achieve convergence. To address this,
we assess the effective sample size (ESS) per second, as it provides insight
into the convergence rate of MCMC chains (Brooks et al., 2011). The auto-
correlation among samples in an MCMC chain introduces uncertainty in
estimating the mean and standard deviation of the posterior distributions. A
higher number of effective samples per second indicates quicker convergence.

We compute the ESS of all chains using the Python library GetDist. When
analyzing Planck, we find an ESS/s of 1.6 for the NUTS sampler, while MCHMC
returns an ESS/s of 3.7, as can be seen in Tab. 4. In the case of ACT, NUTS
achieves 1.9 ESS/s while MCHMC reaches 4.1. Thus in both analyses we see that
MCHMC is twice as fast as NUTS, as expected on the findings of Robnik et al.
(2022). The difference is even bigger for the SPT analysis, since NUTS and
MCHMC delivers an ESS/s of, respectively, 0.09 and 0.335141414The overall
decreased performance for the SPT analysis, compared to the Planck and ACT
ones, is mostly driven by the more expensive likelihood, since more operations
are needed to go from the input $C_{\ell}$ spectra to the log-likelihood
calculation.. Overall, our findings show that the MCHMC sampler has an higher
efficiency compared to the NUTS for the applications considered, as in (Bayer
et al., 2023).151515We would like to emphasize that MicroCanonicalHMC.jl is
still an alpha-software which is being actively developed: performance might
improve with a more stable and refined implementation of MCHMC within this
library.

####  5.2.2 NUTS vs Pathfinder

Comparing NUTS to Pathfinder is a more difficult task since Pathfinder is a VI
method while NUTS is an MCMC method. VI methods do not guarantee asymptotic
convergence to the target distribution, unlike MCMC methods. Instead, VI
approximates the posterior distribution as a combination of possible
distributions (see Sec. 4 for details). Moreover, the concept of ESS is not
applicable to VI methods since they do not inherently produce MCMC chains,
although MCMC-like samples can be generated by sampling the approximated
posterior.

Despite these caveats, we observe excellent agreement between the posteriors
obtained with Pathfinder and NUTS for all the likelihoods considered. This can
be observed in the top half of Figs. 2 and 3. In these plots it is possible to
see that Pathfinder reproduces not only the mean and standard deviation of the
marginal posteriors but also the different correlations between the
parameters. Looking at Tab. 3 we can see that the Pathfinder constraints are
once again within $0.2\,\sigma$ of the CAMB+MH constraints.161616These results
have been generated using the ensemble version of Pathfinder, multi
Pathfinder, that runs several VI processes in parallel to produce better
approximations; in particular, for each multi Pathfinder realization, we run 6
single Pathfinder realization in parallel.

To assess the reliability of Pathfinder, we computed 1,000 realizations of it
and found that the mean obtained for each parameter follows a distribution
which is close to a Gaussian with a standard deviation corresponding to less
than $0.1,\sigma$, measured on each parameter when performing a full MCMC
analysis.

Regarding the SPT analysis, we found Pathfinder to provide unstable results
when applied to this likelihood. For this reason, we decided not to show any
result related to this combination of likelihood and sampler.

Regarding performance, since the ESS is not well-defined for VI methods, we
compare the speed of the two methods by examining pure wall times. As shown in
Tab. 4, Pathfinder successfully approximates the posteriors of the Planck and
ACT analyses in around 10 seconds. This represents a 5 to 6 orders of
magnitude improvement compared to the Cobaya analyses using CAMB and a 2
orders of magnitude improvement compared to using Capse.jl with either NUTS or
MCHMC. However, we emphasize that the robustness of variational methods can
significantly decrease when the posterior is everything but a multivariate
normal distribution. However, even in these cases, methods such as Pathfinder
still can prove useful as they can be used to initialize standard MCMCs,
reducing the burn-in phase even when considering challenging posteriors with
funnel-like geometries (Zhang et al., 2022). This is actually what we did for
the SPT analysis: the Pathfinder draws were unstable and likely to be wrong,
as detected by the PSIS. Nevertheless, they were way closer to the typical set
than random draws from the prior and hence used them to initialize our NUTS
and MCHMC chains.

###  5.3 MOPED compression & Chebyshev decomposition

All the MCMC results shown in the previous section are based on the standard
emulator and likelihood. Although the favourable comparison with CAMB$+$Cobaya
approach, there is still room for some computational improvements. In this
section we show results based on PlanckLite.jl and on a MOPED compressed
likelihood and the Chebyshev-based emulators.

We begin by discussing the standard likelihood results. As discussed in Sec.
2.1, we exploit the Chebyshev decomposition to achieve peak performance. The
PlanckLite.jl likelihood is composed of a sequence of linear operators that
remain independent of cosmological parameters. Utilizing this property, we can
apply these operators to the Chebyshev basis functions and precompute and
store the results. During each step of the MCMC, we only need to multiply the
stored quantities by the expansion coefficients, which solely contain the
cosmological dependence. This approach eliminates the need to repeatedly apply
the linear operator, resulting in enhanced computational efficiency.

Traditionally, optimizing the computation of the likelihood is not a priority
since the bottleneck is the computation of the theory prediction. However,
this is no longer the case when emulators are employed. For Capse.jl, the
$C_{\ell}$’s computation is around 7 times faster than the computations
required to compute the log-likelihood such as the $C\ell$’s binning. This
highlights the performance of our emulators, since the bottleneck is now
represented by the likelihood itself. When using NUTS and MCHMC we find an
ESS/s of 13 and 32 respectively. Meanwhile, the computation time when using
Pathfinder is reduced to 1.5 seconds. This shows an almost uniform performance
improvement by a factor of 8 across the different inference methods
considered.

Finally, we employ a MOPED likelihood (Prince & Dunkley, 2019), a data
compression technique which reduces the dimensions of the data vector to one
number per parameter of interest, requiring an even lower amount of
computational resources. We find an exquisite agreement between the
PlanckLite.jl and the MOPED posteriors, with differences on marginalized
parameters smaller than $0.1\,\sigma$ as shown in Fig.6, but with an improved
computational performance: NUTS and MCHMC reach an ESS/s of, respectively, 60
and 150, while Pathfinder runs in 0.5 seconds. Compared to the performance
obtained when using the standard emulator and PlanckLite.jl we report an
improvement of a factor of 40-50. This brings down the wall time for the
Planck analysis to a record one second. The finding of this section are
summarized in Table 5.

| Planck
---|---
| ESS/s | Comp. (s)
PlanckLite.jl+NUTS | 13 | 1031
MOPED+NUTS | 60 | 228
PlanckLite.jl+MCHMC | 32 | 605
MOPED+MCHMC | 150 | 157
PlanckLite.jl+Pathfinder | N/A | 1.5
MOPED+Pathfinder | N/A | 0.5
Table 5: Shows the ESS per second and computation time achieved by the NUTS,
MCHMC and Pathfinder inference algorithms when performing the Planck Analysis.
Note that here we are using the Chebyshev based emulator.

##  6 Comparison with other emulators

In this section we compare the emulators presented in this work to two closely
related alternatives in the literature, ClassNet (Albers et al., 2019) and the
recently developed CosmoPower171717Unless otherwise stated, in this section we
consider the CosmoPower implementation of Bolliet et al. (2023). (Spurio
Mancini et al., 2022; Bolliet et al., 2023). ClassNet emulates the source
functions required to compute the CMB angular power spectra using Python-based
neural network emulators. The result of the emulation is then passed to CLASS,
which computes the line-of-sight integrals. The main advantage of this
approach is that ClassNet does not need to include $A_{s}$ and $n_{s}$ in the
emulation space, since the primordial power spectrum is analytically computed
and integrated against the emulated source functions. CosmoPower is a
TensorFlow based emulator that follows a design philosophy closer to Capse.jl,
as it directly emulates the CMB angular power spectra. Initially, it focused
on the Planck analysis but it was later extended to
$\ell_{\mathrm{max}}=10,000$ and the same high-accuracy settings employed here
(Bolliet et al., 2023).

Comparing emulators is challenging due to their distinct architectures and
optimization goals. Here we highlight the differences between ClassNet,
CosmoPower and Capse.jl, explaining the pros and cons of each framework.

We decided not to perform a comparison with other frameworks, such as those
presented in Nygaard et al. (2023); Günther (2023); To et al. (2023) as they
have a different goal: performing an analysis using the lowest possible number
of Boltzmann code evaluations using active learning schemes.

In order to compare Capse.jl, ClassNet, and CosmoPower we will use several
metrics. First, we will compare the computational performance of each
emulator. Second, we will compare the training efficiency. Third, we will look
at the flexibility of the emulator, i.e. how often would it need to be
retrained based on the changes to the analysis that one wishes to perform.
Finally, we will also study the efficiency of the sampling algorithms
compatible with each emulator.

We report the speedup obtained by the three emulation frameworks when compared
to a traditional Boltzmann solver in the first row of Tab. 6. In this table,
we can see that ClassNet, CosmoPower, and Capse.jl emulators report a speedup
of a factor of 3,181818We report that in Günther et al. (2022) the authors
performed a Planck data analysis, without focusing on high-precision settings
in which they would likely achieve a higher speed-up. 1,000, and 1,000,000
respectively. We would like to emphasize that the numbers shown in the
comparison were taken directly from the ClassNet and CosmoPower release
papers; their performance might improve with an hyperparameter optimization.

| ClassNet | CosmoPower | Capse.jl
---|---|---|---
Speedup | 3 | 1,000 | 1,000,000
Architecture | NN∗ | 4x512-NN + PCA | 5x64-NN + Chebyshev
Training points | $10^{4}$ | $10^{5}$ | $10^{4}$
Tessellation | Hyperellipsoid | LHC | LHC
Emulation domain
$\ln 10^{10}A_{\mathrm{s}}$ | N/A | [2.5, 2.5] | [2.5, 3.5]
$n_{\mathrm{s}}$ | N/A | [0.8812, 1.0492] | [0.88, 1.06]
$\tau$ | $N(0.0561,0.0071)$ | [0.02, 0.12] | [0.02, 0.12]
$H_{0}[\mathrm{km}/s/\mathrm{Mpc}]$ | $N(67.66,0.42)$ | [39.99, 100.01] | [40, 100]
$\omega_{\mathrm{b}}$ | $N(0.02242,0.00014)$ | [0.0193, 0.02533] | [0.0193, 0.02533]
$\omega_{\mathrm{c}}$ | $N(0.11933,0.00091)$ | [0.08, 0.2] | [0.08, 0.2]
$\ell_{\mathrm{max}}$ | 2,500 | 11,000 | 5,000
Table 6: A comparison between the emulating frameworks ClassNet, CosmoPower,
and Capse.jl. Note that the architecture of ClassNet involves 7 different NN
with 3 to 4 layers each adding to a total 5366 neurons. Regarding CosmoPower,
here we refer to the emulator presented in Bolliet et al. (2023).

Comparing the training efficiency of the three emulators is not
straightforward, primarily due to two main differences between them. Firstly,
the number of inputs (i.e. the cosmological parameters) and their ranges vary
from emulator to emulator. While CosmoPower and Capse.jl were trained on the
same number of input features, ClassNet requires two fewer inputs,
significantly reducing the volume of training data needed to cover its domain.
Secondly, the tessellation method to cover the parameter space can also
differ. CosmoPower and Capse.jl use a Latin hypercube to span their domains,
while ClassNet was trained on a hyperellipsoid centered around the best-fit
Planck TT,TE,EE+lowE+lensing+BAO $\Lambda$CDM cosmology (Aghanim et al.,
2020b), spanning 6 $\sigma$’s. Regarding the number of learned features,
ClassNet learns source functions on a $k$-grid with $\mathcal{O}(700)$
elements, chosen to accurately compute the Planck $C_{\ell}$’s. In contrast,
Capse.jl computes all $C_{\ell}$’s up to $\ell_{\mathrm{max}}=5,000$ while
CosmoPower reaches an $\ell_{\mathrm{max}}$ of 11,000, making them both
suitable for next-generation CMB experiments.

Considering the aforementioned caveats, training Capse.jl requires $10^{4}$
training points, which is of the same order of magnitude as ClassNet, and an
order of magnitude less than CosmoPower Bolliet et al. (2023). In terms of
network architecture, CosmoPower employs a one-size-fits-all approach with 4
fully connected hidden layers, each containing 512 neurons. ClassNet, on the
other hand, tailors its networks for each emulator, utilizing different
architectures with 2-4 hidden layers and 100-500 neurons. Capse.jl employs a
consistent architecture across all networks, consisting of 5 hidden layers
with 64 neurons. In terms of training time, Capse.jl takes slightly less than
one hour using a CPU, while the training of CosmoPower lasted one hour using a
GPU. After training, Capse.jl is slightly more precise than CosmoPower for the
$\phi\phi$ coefficients, while CosmoPower is 2-3 times more precise for the
$TT$, $TE$, $EE$ coefficients.

Regarding flexibility, ClassNet is the most versatile of the emulators. As
previously described, it only emulates the transfer functions and computes the
line-of-sight integrals using CLASS. This allows the user to change the
computation of the primordial power spectrum or the method to compute the
line-of-sight integrals without having to retrain the emulator, as it would be
the case for Capse.jl and CosmoPower.

Finally, we studied the efficiency of sampling algorithms compatible with each
emulator. The wall time users will experience depends on both the speed of the
emulator and the sampling algorithms available for Bayesian inference.
ClassNet employs traditional gradient-free sampling methods while Capse.jl and
CosmoPower are compatible with state-of-the-art gradient-based algorithms. The
original CosmoPower was written in TensorFlow and, coupled to
TensorFlowProbability191919https://www.tensorflow.org/probability, has in
principle access to a rich variety of gradient-based samplers. Moreover,
CosmoPower was recently re-written in JAX (Piras & Spurio Mancini, 2023) and
embedded in the JAX-COSMO framework, explicitly designed for gradient-based
sampling methods. On the other hand, using these methods with ClassNet would
be technically challenging due to its mixed C and Python implementation.
Capse.jl was coupled with the Turing.jl PPL, enabling the use of NUTS, MCHMC,
and Pathfinder inference algorithms, fully leveraging the differentiability of
our emulators.

When comparing the wall time required to compute converging chains, ClassNet
reports a speedup by a factor of 3. In Spurio Mancini et al. (2022), a Planck
analysis was performed in $\mathcal{O}(10\,s)$ using a GPU; however, as can be
seen on the official notebook prepared by the authors, the chains were
initialized on the typical set, using the knowledge coming from previous
analyses. Capse.jl is able to perform a Planck analysis in around a minute
with MCMC methods or $\mathcal{O}(1\,s)$ with variational methods. We
emphasize that Capse.jl’s runs were performed randomly initializing our
algorithms from the prior. While there is nothing conceptually wrong in
initializing the chains on the typical set, we do think that starting using
random draws from the prior is more representative of what happens when
analyzing a new dataset for the first time.

Finally, we want to emphasize that the results shown in this comparison are a
direct consequence of the design choices of the emulators. ClassNet is the
most flexible framework, but this comes at the expense of a lower speedup.
CosmoPower delivers a higher speedup than ClassNet and has already been
trained to analyze next-generation CMB surveys, but its training requires an
expensive NN architecture and training dataset. Capse.jl is faster than
CosmoPower, can be used to analyze actual datasets and is cheaper to train,
but its multipole range and accuracy are slightly lower than the one of
CosmoPower. According to the scenario considered, one of the three emulators
considered will better suit the requirements of the analysis. We postpone to a
further publication a detailed investigation of the origin of the higher
runtime performance of Capse.jl.

##  7 Conclusions

In this paper, we introduced Capse.jl, a novel emulator framework for CMB
angular power spectra. We discussed the preprocessing steps and the
characteristics of the NN library employed, highlighting the distinctions from
similar emulators in the literature (Spurio Mancini et al., 2022; Günther et
al., 2022). We additionally developed an emulator based on Chebyshev
polynomial decomposition, which significantly enhanced the computational
performance of the original emulator by reducing the number of output
features. Leveraging the Chebyshev linear decomposition further improved the
efficiency of the log-likelihood evaluation.

To validate the accuracy of our emulators, we conducted two tests. First, we
compared Capse.jl and CAMB computations over 20,000 combinations of input
parameters across the entire emulation range. Our emulators exhibited errors
below $0.1\,\sigma$ for all scales relevant to future CMB surveys. Second, we
analyzed the Planck and ACT DR4 nuisance-marginalized datasets as well as the
full 2018 SPT-3G dataset, comparing the posteriors with those obtained using
CAMB and Cobaya. The derived contours showed a remarkable agreement, better
than $0.1\,\sigma$.

Next, we examined the performance impact of different sampling algorithms. For
this purpose, we considered several gradient-based samplers, namely NUTS, the
state-of-the-art Hamiltonian Monte Carlo sampler, MicroCanonical Hamiltonian
Monte Carlo, a recently developed sampler by Robnik et al. (2022) and Robnik &
Seljak (2023), and Pathfinder, a state-of-the-art variational inference
algorithm. While CAMB paired with Cobaya required around 100, 480, and 1,600
CPU hours to analyze Planck, ACT and SPT-3G data respectively, Capse.jl
combined with the NUTS algorithm required only 0.5 CPU hours for the Planck
analysis, one CPU hour for the ACT analysis, and 14 hours for the SPT-3G
analysis. Thus our analyses were at least two orders of magnitude more
efficient. When comparing NUTS to MCHMC, we found the latter to be more
efficient by a factor of 2-3, while producing virtually identical posteriors.
Pathfinder was able to perform the Planck and ACT analysis in 1-10 s, with a
very good precision. However, this approach is less robust when compared to
MCMC methods as it lacks asymptotic convergence guarantees.

In this work we also tested a data reduction technique, the Chebyshev
polynomial decomposition. This technique was instrumental in improving
Capse.jl computational performance, as it improved the Planck analysis
efficiency by almost an order of magnitude.

While our emulators and methods exhibits high precision and computational
performance, there remain opportunities for future improvements. First, we aim
to enhance the preprocessing procedure, addressing the reduced precision of
Capse.jl for low-$\ell$ $EE$ 2-point function spectra by employing a proper
$\tau$ rescaling. Implementing a weighted mean square error loss function
(similar to Nygaard et al. (2023)) based on the CMB power spectra variance
could also lead to a more homogeneous emulation error at no additional
computational cost. We also plan to work on refining the Chebyshev polynomial
decomposition approach. Although it demonstrated precision and efficiency
within the Planck multipole range, we intend to explore ways to improve its
robustness when considering higher $\ell_{\mathrm{max}}$ values. In
particular, some preliminary tests show that we can reach the same
$\ell_{\mathrm{max}}$ of CosmoPower, while maintaining the computational
efficiency of our framework and even improving the precision of our emulators
with respect to what we showed in this paper.

Working on the Chebyshev polynomial decomposition approach might not only be
relevant for CMB analyses but also for spectroscopic galaxy clustering
analyses. The estimators used to extract cosmological information from galaxy
surveys (Feldman et al., 1994; Yamamoto et al., 2006) suffer from an imprint
of the survey window geometry. In order to take into account this effect, the
community has developed several methods tailored for power spectrum (Blake et
al., 2013; Beutler et al., 2017) estimation. However, only recently this issue
has been addressed for the bispectrum (Pardede et al., 2022; Alkhanishvili et
al., 2023). We plan to investigate whether the Chebyshev polynomial
decomposition can be useful for this kind of application.

## Acknowledgements

We would like thank Boris Bolliet and Alessio Spurio Mancini for
clarifications regarding CosmoPower. MB is expecially grateful to Martin White
for the initial suggestion on employing the Chebyshev polynomial
decomposition, to Marius Millea for his advices while coding the SPT-3G
likelihood in Julia and to Emanuele Castorina, Giulio Fabbian, SImone Ferraro
and Uros Seljak for useful discussions. MB thanks University of California
Berkeley for hospitality during his visit, during which this project was
started. MB acknowledges financial support from INAF MiniGrant 2022. JRZ is
supported by an STFC doctoral studentship. F.B. acknowledges support by the
Department of Energy, Contract DE-AC02-76SF00515.

## References

  * Abate et al. (2012) Abate A., et al., 2012, Large Synoptic Survey Telescope: Dark Energy Science Collaboration (arXiv:1211.0310)
  * Abazajian et al. (2019) Abazajian K., et al., 2019, CMB-S4 Decadal Survey APC White Paper (arXiv:1908.01062)
  * Abbott et al. (2018) Abbott T. M. C., et al., 2018, Phys. Rev. D, 98, 043526
  * Ade et al. (2019) Ade P., et al., 2019, JCAP, 02, 056
  * Aghanim et al. (2020a) Aghanim N., et al., 2020a, Astron. Astrophys., 641, A5
  * Aghanim et al. (2020b) Aghanim N., et al., 2020b, Astron. Astrophys., 641, A6
  * Aiola et al. (2020) Aiola S., et al., 2020, JCAP, 12, 047
  * Albers et al. (2019) Albers J., Fidler C., Lesgourgues J., Schöneberg N., Torrado J., 2019, JCAP, 09, 028
  * Alkhanishvili et al. (2023) Alkhanishvili D., Porciani C., Sefusatti E., 2023, Astron. Astrophys., 669, L2
  * Alvarez et al. (2021) Alvarez M. A., Ferraro S., Hill J. C., Hložek R., Ikape M., 2021, Phys. Rev. D, 103, 063518
  * Angulo et al. (2021) Angulo R. E., Zennaro M., Contreras S., Aricò G., Pellejero-Ibañez M., Stücker J., 2021, Mon. Not. Roy. Astron. Soc., 507, 5869
  * Aricò et al. (2021) Aricò G., Angulo R. E., Zennaro M., 2021, doi:10.12688/openreseurope.14310.2
  * Auld et al. (2007) Auld T., Bridges M., Hobson M. P., Gull S. F., 2007, Mon. Not. Roy. Astron. Soc., 376, L11
  * Balkenhol et al. (2023) Balkenhol L., et al., 2023, Phys. Rev. D, 108, 023510
  * Bayer et al. (2023) Bayer A. E., Seljak U., Modi C., 2023, in 40th International Conference on Machine Learning. (arXiv:2307.09504)
  * Betancourt (2018) Betancourt M., 2018, A Conceptual Introduction to Hamiltonian Monte Carlo (arXiv:1701.02434)
  * Beutler et al. (2017) Beutler F., et al., 2017, Mon. Not. Roy. Astron. Soc., 466, 2242
  * Bezanson et al. (2015) Bezanson J., Edelman A., Karpinski S., Shah V. B., 2015, Julia: A Fresh Approach to Numerical Computing (arXiv:1411.1607)
  * Blake et al. (2013) Blake C., et al., 2013, Mon. Not. Roy. Astron. Soc., 436, 3089
  * Blas et al. (2011) Blas D., Lesgourgues J., Tram T., 2011, JCAP, 07, 034
  * Bolliet et al. (2023) Bolliet B., Spurio Mancini A., Hill J. C., Madhavacheril M., Jense H. T., Calabrese E., Dunkley J., 2023, High-accuracy emulators for observables in $\Lambda$CDM, $N_{\mathrm{eff}}$, $\Sigma m_{\nu}$, and $w$ cosmologies (arXiv:2303.01591)
  * Bonici et al. (2022) Bonici M., Biggio L., Carbone C., Guzzo L., 2022, Fast emulation of two-point angular statistics for photometric galaxy surveys (arXiv:2206.14208)
  * Brooks et al. (2011) Brooks S., Gelman A., Jones G., Meng X.-L., 2011, Handbook of Markov Chain Monte Carlo. CRC press
  * Campagne et al. (2023) Campagne J.-E., et al., 2023, Open J. Astrophys., 6, 1
  * Crill et al. (2020) Crill B. P., et al., 2020, in Space Telescopes and Instrumentation 2020: Optical, Infrared, and Millimeter Wave. SPIE, pp 61–77, doi:10.1117/12.2567224, https://www.spiedigitallibrary.org/conference-proceedings-of-spie/11443/114430I/SPHEREx-NASAs-near-infrared-spectrophotometric-all-sky-survey/10.1117/12.2567224.full
  * Dawson et al. (2013) Dawson K. S., et al., 2013, Astron. J., 145, 10
  * DeRose et al. (2019) DeRose J., et al., 2019, Astrophys. J., 875, 69
  * Donald-McCann et al. (2022a) Donald-McCann J., Beutler F., Koyama K., Karamanis M., 2022a, Mon. Not. Roy. Astron. Soc., 511, 3768
  * Donald-McCann et al. (2022b) Donald-McCann J., Koyama K., Beutler F., 2022b, Mon. Not. Roy. Astron. Soc., 518, 3106
  * Dutcher et al. (2021) Dutcher D., et al., 2021, Phys. Rev. D, 104, 022003
  * Eggemeier et al. (2022) Eggemeier A., Camacho-Quevedo B., Pezzotta A., Crocce M., Scoccimarro R., Sánchez A. G., 2022, Mon. Not. Roy. Astron. Soc., 519, 2962
  * Feldman et al. (1994) Feldman H. A., Kaiser N., Peacock J. A., 1994, Astrophys. J., 426, 23
  * Fendt & Wandelt (2006) Fendt W. A., Wandelt B. D., 2006, Astrophys. J., 654, 2
  * Gammal et al. (2022) Gammal J. E., Schöneberg N., Torrado J., Fidler C., 2022, Fast and robust Bayesian Inference using Gaussian Processes with GPry (arXiv:2211.02045)
  * 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
  * Günther (2023) Günther S., 2023, Uncertainty-aware and Data-efficient Cosmological Emulation using Gaussian Processes and PCA (arXiv:2307.01138)
  * Günther et al. (2022) Günther S., Lesgourgues J., Samaras G., Schöneberg N., Stadtmann F., Fidler C., Torrado J., 2022, JCAP, 11, 035
  * Heavens et al. (2000) Heavens A., Jimenez R., Lahav O., 2000, Mon. Not. Roy. Astron. Soc., 317, 965
  * Heavens et al. (2017) Heavens A., Sellentin E., de Mijolla D., Vianello A., 2017, Mon. Not. Roy. Astron. Soc., 472, 4244
  * Hoffman & Gelman (2011) Hoffman M. D., Gelman A., 2011, The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo (arXiv:1111.4246)
  * Jimenez et al. (2004) Jimenez R., Verde L., Peiris H., Kosowsky A., 2004, Phys. Rev. D, 70, 023005
  * Kamionkowski et al. (1997) Kamionkowski M., Kosowsky A., Stebbins A., 1997, Phys. Rev. D, 55, 7368
  * Kingma & Ba (2017) Kingma D. P., Ba J., 2017, Adam: A Method for Stochastic Optimization (arXiv:1412.6980)
  * Knabenhans et al. (2019) Knabenhans M., et al., 2019, Mon. Not. Roy. Astron. Soc., 484, 5509
  * Knabenhans et al. (2021) Knabenhans M., et al., 2021, Mon. Not. Roy. Astron. Soc., 505, 2840
  * Knox (1995) Knox L., 1995, Phys. Rev. D, 52, 4307
  * Kucukelbir et al. (2016) Kucukelbir A., Tran D., Ranganath R., Gelman A., Blei D. M., 2016, Automatic Differentiation Variational Inference (arXiv:1603.00788)
  * Laureijs et al. (2011) Laureijs R., et al., 2011, Euclid Definition Study Report (arXiv:1110.3193)
  * Lawrence et al. (2017) Lawrence E., et al., 2017, Astrophys. J., 847, 50
  * Levi et al. (2013) Levi M., et al., 2013, The DESI Experiment, a whitepaper for Snowmass 2013 (arXiv:1308.0847)
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473
  * Manrique-Yus & Sellentin (2020) Manrique-Yus A., Sellentin E., 2020, Mon. Not. Roy. Astron. Soc., 491, 2655
  * McKay et al. (1979) McKay M. D., Beckman R. J., Conover W. J., 1979, Technometrics, 21, 239
  * Mead et al. (2020) Mead A., Brieden S., Tröster T., Heymans C., 2020, Mon. Not. Roy. Astron. Soc.
  * Mootoovaloo et al. (2020) Mootoovaloo A., Heavens A. F., Jaffe A. H., Leclercq F., 2020, Mon. Not. Roy. Astron. Soc., 497, 2213
  * Mootoovaloo et al. (2022) Mootoovaloo A., Jaffe A. H., Heavens A. F., Leclercq F., 2022, Astron. Comput., 38, 100508
  * Nygaard et al. (2023) Nygaard A., Holm E. B., Hannestad S., Tram T., 2023, JCAP, 05, 025
  * Pardede et al. (2022) Pardede K., Rizzo F., Biagetti M., Castorina E., Sefusatti E., Monaco P., 2022, JCAP, 10, 066
  * Piras & Spurio Mancini (2023) Piras D., Spurio Mancini A., 2023, Open J. Astrophys.
  * Prince & Dunkley (2019) Prince H., Dunkley J., 2019, Phys. Rev. D, 100, 083502
  * Ranganath et al. (2013) Ranganath R., Gerrish S., Blei D. M., 2013, Black Box Variational Inference (arXiv:1401.0118)
  * Robnik & Seljak (2023) Robnik J., Seljak U., 2023, Microcanonical Langevin Monte Carlo (arXiv:2303.18221)
  * Robnik et al. (2022) Robnik J., Luca G. B. D., Silverstein E., Seljak U., 2022, Microcanonical Hamiltonian Monte Carlo (arXiv:2212.08549)
  * Sobrin et al. (2022) Sobrin J. A., et al., 2022, Astrophys. J. Supp., 258, 42
  * Spergel et al. (2015) Spergel D., et al., 2015, Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report (arXiv:1503.03757)
  * Spurio Mancini et al. (2022) Spurio Mancini A., Piras D., Alsing J., Joachimi B., Hobson M. P., 2022, Mon. Not. Roy. Astron. Soc., 511, 1771
  * To et al. (2023) To C.-H., Rozo E., Krause E., Wu H.-Y., Wechsler R. H., Salcedo A. N., 2023, JCAP, 01, 016
  * Trefethen (2012) Trefethen L. N., 2012, Approximation Theory and Approximation Practice.. SIAM
  * Vehtari et al. (2022) Vehtari A., Simpson D., Gelman A., Yao Y., Gabry J., 2022, Pareto Smoothed Importance Sampling (arXiv:1507.02646)
  * Wang (2023) Wang H., 2023, Analysis of error localization of Chebyshev spectral approximations (arXiv:2106.03456)
  * Yamamoto et al. (2006) Yamamoto K., Nakamichi M., Kamino A., Bassett B. A., Nishioka H., 2006, Publ. Astron. Soc. Jap., 58, 93
  * Zhang et al. (2022) Zhang L., Carpenter B., Gelman A., Vehtari A., 2022, Pathfinder: Parallel quasi-Newton variational inference (arXiv:2108.03782)
  * (73) de Jong J. T. A., et al.,, The Kilo-Degree Survey

##  Appendix A Affine transformation of Multivariate distribution

Let us consider a random variable $X$, distributed according to:

$X\sim\mathcal{N}(\mu_{X},\Sigma_{X}).$ (12)

If we define a new variable:

$Y=LX+u,$ (13)

then $Y$ follows:

$Y\sim\mathcal{N}(\mu_{Y},\Sigma_{Y})=\mathcal{N}(u+L\mu_{X},L\Sigma_{X}L^{T}).$
(14)

This can be proved as follows:

$\mu_{Y}=\mathbb{E}[Y]=\mathbb{E}[LX+u]=\mathbb{E}[LX]+u=L\mu_{X}+u$ (15)
$\displaystyle\Sigma_{Y}$
$\displaystyle=\mathbb{E}\left[\left(Y-\mu_{Y}\right)\left(Y-\mu_{Y}\right)^{\top}\right]$
(16)
$\displaystyle=\mathbb{E}\left[\left(LX+u-L\mu_{X}-u\right)\left(LX+u-L\mu_{X}-u\right)^{\top}\right]$
$\displaystyle=\mathbb{E}\left[\left(L\left(X-\mu_{X}\right)\right)\left(L\left(X-\mu_{X}\right)\right)^{\top}\right]$
$\displaystyle=\mathbb{E}\left[L\left(X-\mu_{X}\right)\left(X-\mu_{X}\right)^{\top}L^{\top}\right]$
$\displaystyle=L\mathbb{E}\left[\left(X-\mu_{X}\right)\left(X-\mu_{X}\right)^{\top}\right]L^{\top}$
$\displaystyle=L\Sigma_{X}L^{\top}.$

In this work we consider likelihood of the form:

$d\sim\mathcal{N}\left(t(\theta),\Sigma\right)$ (17)

where $d$ represents either Planck or ACT data and $\Sigma$ their covariance.
Given the Capse.jl emulator performance, computing the inverse of the
covariance matrix gives a not negligible computational overhead. Eq. 14 is
interesting for our discussion as it gives us a way to avoid computing the
inverse of the covariance matrix at each step of the MonteCarlo. In
particular, if we define:

$\Gamma=\mathrm{sqrt}(\Sigma),$ (18) $\mathcal{D}=\mathrm{inv}(\Gamma)d,$ (19)

then:

$\mathcal{D}\sim\mathcal{N}(\mathrm{inv}(\Gamma)t(\theta),\mathds{I}).$ (20)

We use this parametrization in our likelihoods.

##  Appendix B Simons Observatory residuals distribution

Figure 5: Distribution of the relative difference between Capse.jl and the
high-accuracy predictions from CAMB across the four 2-point statistics
emulated in this work, normalized to the forecasted SO statistical
uncertainties ($f_{\mathrm{sky}}=0.4$). The differences are in units of cosmic
variance below $\ell<30$. The shaded regions denote the 68%, 95%, and 99%
percentiles of the distributions, respectively.

##  Appendix C Planck MOPED likelihood contours

Figure 6: Contours obtained analyzing the Planck premarginalized dataset. As
in Fig. 2, we compare the CAMB chains with the Capse.jl ones, the main
difference being that we used the MOPED likelihood for the emulator-based
runs. Also in this case, the comparison of the contours is very good, being
the difference on marginalized parameters smaller than $0.1\,\sigma$.