# OLÉ– Online Learning Emulation in Cosmology

Sven Günther, Lennart Balkenhol, Christian Fidler, Ali Rida Khalife, Julien Lesgourgues, Markus R Mosbech, Ravi Kumar Sharma

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sven Günther, Lennart Balkenhol, Christian Fidler, Ali Rida Khalife, Julien Lesgourgues, et al., OLÉ–Online Learning Emulation in Cosmology. Journal of Cosmology and Astroparticle Physics, 2025, 09, pp.059. 10.1088/1475-7516/2025/09/059. [hal-05029734]

HAL Id: hal-05029734

https://hal.science/hal-05029734v1

Submitted on 31 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
PAPER • OPEN ACCESS

OLÉ — Online Learning Emulation in cosmology

To cite this article: Sven Günther *et al* JCAP09(2025)059

View the [article online](#) for updates and enhancements.

You may also like

- Dark Quest. I. Fast and Accurate Emulation of Halo Clustering Statistics and Its Application to Galaxy Clustering
Takahiro Nishimichi, Masahiro Takada, Ryuichi Takahashi et al.

- The Aemulus Project. III. Emulation of the Galaxy Correlation Function
Zhongxu Zhai, Jeremy L. Tinker, Matthew R. Becker et al.

- Applications of emulation and Bayesian methods in heavy-ion physics
Jean-François Paquet
# OLÉ — Online Learning Emulation in cosmology

Sven Günther<sup>D,a</sup>, Lennart Balkenhol,<sup>b</sup> Christian Fidler,<sup>a</sup> Ali Rida Khalife<sup>D,b</sup>

Julien Lesgourgues<sup>D,a</sup>, Markus R. Mosbech<sup>D,a</sup> and Ravi Kumar Sharma<sup>a</sup>

E-mail: sven.guenther@rwth-aachen.de, lennart.balkenhol@iap.fr, fidler@physik.rwth-aachen.de, ridakhal@iap.fr, lesgourg@physik.rwth-aachen.de, mosbech@physik.rwth-aachen.de, rksharma@physik.rwth-aachen.de

**ABSTRACT:** In this work, we present OLÉ, a new online learning emulator for use in cosmological inference. The emulator relies on Gaussian Processes and Principal Component Analysis for efficient data compression and fast evaluation. Moreover, OLÉ features an automatic error estimation for optimal active sampling and online learning. All training data is computed on-the-fly, making the emulator applicable to any cosmological model or dataset. We illustrate the emulator's performance on an array of cosmological models and data sets, showing significant improvements in efficiency over similar emulators without degrading accuracy compared to standard theory codes. We find that OLÉ is able to considerably speed up the inference process, increasing the efficiency by a factor of 30 – 350, including data acquisition and training. Typically the runtime of the likelihood code becomes the computational bottleneck. Furthermore, OLÉ emulators are differentiable; we demonstrate that, together with the differentiable likelihoods available in the `candl` library, we can construct a gradient-based sampling method which yields an additional improvement factor of 4. OLÉ can be easily interfaced with the popular samplers `MontePython` and `Cobaya`, and the Einstein-Boltzmann solvers `CLASS` and `CAMB`. OLÉ is publicly available at https://github.com/svenguenther/OLE.

**KEYWORDS:** cosmological parameters from CMBR, cosmological parameters from LSS, Machine learning, Statistical sampling techniques

ARXIV ePRINT: 2503.13183

<figure><img src="image_3.png" /></figure>

<footnote><sup>a</sup>Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, Sommerfeldstr. 16, D-52056 Aachen, Germany</footnote>

<footnote><sup>b</sup>Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>OLÉ strategy</td><td>2</td></tr><tr><td>2.1</td><td>Learning strategy</td><td>3</td></tr><tr><td>2.2</td><td>Emulator design</td><td>5</td></tr><tr><td>2.3</td><td>Code design</td><td>11</td></tr><tr><td>2.4</td><td>OLÉ in MCMC analyses</td><td>11</td></tr><tr><td>3</td><td>Examples with Metropolis Hastings sampling</td><td>12</td></tr><tr><td>3.1</td><td>ΛCDM with current data using Cobaya and CLASS</td><td>15</td></tr><tr><td>3.2</td><td>Extended cosmology with current data using Cobaya and CAMB</td><td>17</td></tr><tr><td>3.3</td><td>Stage-IV LSS forecast for an extended cosmology using MontePython and CLASS</td><td>20</td></tr><tr><td>3.4</td><td>Early dark energy with current data using MontePython and TriggerCLASS</td><td>21</td></tr><tr><td>4</td><td>Staged NUTS sampling with differentiable likelihoods</td><td>24</td></tr><tr><td>5</td><td>Conclusion</td><td>26</td></tr><tr><td>A</td><td>Using OLÉ</td><td>27</td></tr><tr><td>A.1</td><td>Implemented sampler</td><td>27</td></tr><tr><td>A.2</td><td>MontePython interface</td><td>27</td></tr><tr><td>A.3</td><td>Cobaya interface</td><td>29</td></tr><tr><td>B</td><td>Precision settings</td><td>30</td></tr><tr><td>B.1</td><td>Default precision parameters</td><td>31</td></tr><tr><td>B.2</td><td>Extended cosmology example</td><td>32</td></tr><tr><td>B.3</td><td>Stage-IV LSS forecast example</td><td>32</td></tr><tr><td>B.4</td><td>Early dark energy example</td><td>32</td></tr><tr><td>B.5</td><td>Staged NUTS sampling example</td><td>32</td></tr><tr><td>C</td><td>Cosmological parameters</td><td>32</td></tr><tr><td>D</td><td>Introduction to Gaussian processes</td><td>33</td></tr><tr><td>D.1</td><td>Overview</td><td>33</td></tr><tr><td>D.2</td><td>Gaussian process regression</td><td>34</td></tr><tr><td>D.3</td><td>Kernel functions</td><td>34</td></tr><tr><td>D.4</td><td>Advantages and limitations</td><td>34</td></tr></table>

# 1 Introduction

With the ever increasing precision and complexity of data sets gathered by modern cosmological experiments [1–7], the requirements for computational accuracy rise in tandem. For the use of
these data sets to remain feasible for inference, we must also improve the efficiency with which
we compute our theoretical predictions without degrading their accuracy. Emulators are a
common solution to this problem, allowing for fast calls of a function trained to reproduce the
output of a more computationally expensive theory code. These are often pre-trained [8–12],
relying on a library of pre-computed calls to the theory code for a variety of input parameters.
This offers the advantage of fast predictions out of the box, at the cost of flexibility, as the
emulator can only be assumed to be accurate within the parameter space it is trained on.
In particular, this is a problem in cases where the relevant parameter space or the required
accuracy is a priori unknown, leading to overly conservative coverage of the parameter space.
This also limits the usefulness of pre-trained emulators when considering models outside
of the standard Cosmological Constant + Cold Dark Matter (ΛCDM) cosmology and its
most common extensions.

In this work, we present an alternative type of emulator, relying on online learning<fnref n="1" /> to
work jointly with the theory code to produce training data and train on-the-fly for whichever
model and parameter space it is applied to. This method is an extension of the one presented
in [14], which was used in [15] as a proof of concept for its functionality. Our method
allows to fit any combination of data, since the training is performed specifically for the
observables relevant to a given data set, while the emulator accuracy is matched automatically
to a fraction of the observational errors. This approach works extremely well to speed up
inference pipelines relying on many calls to already relatively fast theory codes, such as the
Cosmic Linear Anisotropy Solving System, CLASS [16, 17], or the Code for Anisotropies in
the Microwave Background, CAMB [18, 19]. Our emulator is called OLÉ,<fnref n="2" /> standing for Online
Learning Émulator. In this work, we describe the design of OLÉ, and showcase its applicability
to inference using different cosmological data sets.

The paper is structured as follows. In section 2 we describe the emulator and its design. In section 3, we present a few examples of use cases of the emulator, and we compare its results and performance to traditional methods. We develop and benchmark a gradient-based staged sampling approach in section 4. We present our conclusions in section 5. In appendices A and B, we provide more details on the code and how to use it.

## 2 OLÉ strategy

Emulators find widespread use in research that relies on computationally expensive simulations to perform parameter inference. Applications in cosmology are particularly common for N-body codes [9, 12, 20–42], hydrodynamic simulations [10, 44–48] or emulators of Einstein-Boltzmann-Solvers [8, 11, 49–58]. Beyond the architecture — such as neural networks or other machine learning techniques — one of the key design choices for an emulator is the distribution of support points used for training. A larger number of support points can enhance accuracy and enable reliable predictions across a broader parameter space. However, this comes

<footnote n="1">¹Online learning is a training strategy relying on training data becoming available sequentially, allowing the model to be trained progressively as more data becomes available. One example is by training on (a subset of) steps in a sampling algorithm (see [13, 14] for more details).</footnote>

<footnote n="2">²Available at https://github.com/svenguenther/OLE, with a detailed documentation available at https://ole.readthedocs.io.</footnote>
at a significant computational cost, affecting both training and emulation, particularly in high-dimensional parameter spaces.

One way to reduce the number of training samples is through active sampling, where the selection of additional support points is optimized based on the emulator’s current state. This approach aims to improve accuracy in regions with high prediction uncertainty and relevance to the specific task the emulator is designed for. In the context of parameter inference, new support points should be concentrated in regions where high likelihoods are possible, but the emulator’s accuracy remains insufficient. As the relevant parameter space only becomes available during the inference run itself, it comes natural to use online learning, thus training the emulator during the inference on the samples it has seen before to have a suited emulator for all future samples to come.

OLÉ makes use of both active sampling and online learning, with the aim of deploying fast and compact emulators that rely on small but important sets of samples. Similarly to CONNECT [54], OLÉ utilizes the sequential data of the inference process to train the emulator. This ensures that the collected samples are already concentrated in the most relevant regions of the parameter space. Unlike CONNECT, OLÉ includes a built-in accuracy check, which can be applied to every individual call. This allows OLÉ to dynamically switch between the emulator and the standard theory code, to ensure accurate results and improve itself on-the-fly by adding support points whenever it is deemed necessary. To the best of our knowledge, this is the first emulator for cosmology using such an approach.

This strategy allows for an extremely data efficient construction of emulators that are trained on a set of support points which is a magnitude smaller than for comparable emulators.<fnref n="3" /> OLÉ requires no pre-training, as all data acquisition and training is done during the inference process itself.

## 2.1 Learning strategy

In this section, we outline the main features of the OLÉ algorithm, which is summarized via pseudocode in Algorithm 1. We consider N ordered samples of an arbitrary sampling algorithm. For each sample i, we have a set of parameter values $\theta_i$, observables $x_i$, and a likelihood value $l_i$, such that a complete chain is:

$$ \vec{X} = (x_1, \dots, x_N), \qquad (2.1) $$

$$ \vec{\Theta} = (\theta_1, \dots, \theta_N), \qquad (2.2) $$

$$ \vec{L} = (l_1, \dots, l_N). \qquad (2.3) $$

For example, in a cosmological context, $\theta_i$ would correspond to cosmological parameters, while $x_i$ could represent the power spectra of the Cosmic Microwave Background (CMB). The parameters, observables, and likelihood values are related to each other via a theory

<footnote n="3">³For example, the OLÉ run outlined in section 3.2 on an extended set of cosmological parameters required 940 calls to CLASS, of which 240 were stored in the cache for training. Comparable analyses with CONNECT used ~ 40 times as many CLASS calls [54].</footnote>
code $\mathcal{T}$ and a likelihood code $\mathcal{L}$,

$$x_i = \mathcal{T}(\theta_i), \qquad (2.4)$$

$$l_i = \mathcal{L}(x_i). \qquad (2.5)$$

The emulator $\mathcal{E}$ can, once trained, alternatively perform the role of the theory code. We denote emulated quantities via $\tilde{x}_i$

$$\tilde{x}_i = \mathcal{E}(\theta_i). \qquad (2.6)$$

OLÉ also defines an error estimation function `IsEmulatorAccurate(\theta_i)` in charge of testing whether relying on the emulator is a valid choice for a set of parameters $\theta_i$. This function estimates a statistic $\Delta_{\log l}$ that is proportional to the relative error on the likelihood induced by the use of the emulator. It returns a boolean flag, stating whether the emulator is acceptable or not on the basis of both $\Delta_{\log l}$ and $l_i$ itself, requiring more accuracy in the region where the likelihood is high. Further details on the construction of this function and criterion can be found in section 2.2.3.

The algorithm begins by initializing the OLÉ internal cache, which will store the parameters $\theta_j$, the observables $x_j$, and the corresponding likelihoods $l_j$. The process of filling the cache occurs in two stages. In the first stage, OLÉ populates the cache with all points generated by the sampler that fall within a specified margin of the highest likelihood found so far.

During this burn-in stage, the margin$^4$ discriminates between relevant samples and outliers that are removed as they divert the focus of the emulator from the relevant parameter space. Once a sufficient number of support points have been collected, the emulator is trained for the first time.

In the second stage, for each new sample $\theta_i$, the emulator is evaluated and tested using the function `IsEmulatorAccurate(\theta_i)` defined above. If its accuracy flag returns true, the emulated result $\tilde{x}_i$ is used to compute the likelihood, and no additional training data is required. In contrast, if the emulator fails the accuracy test, this means that the emulator lacks precision in this region of parameter space. In that case, OLÉ calls the theory code $\mathcal{T}$ to recompute the observables at this point instead, potentially adding the new sample as a support point to the cache and updating the emulator. This ensures good emulator performance for subsequent samples in this part of the parameter space. This approach guarantees that the emulator selectively adds points only in regions where its accuracy is insufficient, maintaining a low-correlation training set and thereby optimizing the learning process.

Ultimately, the emulator becomes accurate across the full region where the likelihood is high. Then, the function `IsEmulatorAccurate(\theta_i)` returns a positive result at each new point drawn from the sampler. In this final stage, the theory code no longer needs to be called, and the emulator greatly speeds up the sampling rate for the remainder of the run.

<footnote>$^4$This margin is computed from quantile function of the $\chi^2$-distribution with the argument of the probability $p$ and the dimensionality of the inference task. Thus, we would expect the occurrence of random samples of the posterior within the margin to happen with a likelihood of $1-p$. The probability $p$ is inferred from user-defined precision parameters, see appendix B.</footnote>
**Algorithm 1 OLÉ sampling algorithm.**

Require: (functions) Sampler, FindLikelihoodThreshold, IsEmulatorAccurate

Require: (OLÉ emulator) $\mathcal{E}$

Require: (theory) $\mathcal{T}$

Require: (likelihood) $\mathcal{L}$

Require: (int) $min\_data\_points$ (OLÉ hyperparameter), $N_{\text{steps}}$ (number of sampler steps)

1: cache $\leftarrow$ list()

2: for $i \leftarrow 1, N_{\text{steps}}$ do

3: $\theta_i \leftarrow$ Sampler()

4: if $size(cache) \le min\_data\_points$ then

5: $x_i \leftarrow \mathcal{T}(\theta_i)$

6: $l_i \leftarrow \mathcal{L}(x_i)$

7: $l_{\text{threshold}} \leftarrow \text{FindLikelihoodThreshold}(\text{cache})$

8: if $l_i > l_{\text{threshold}}$ then

9: cache.append($x_i$)

10: end if

11: if $size(cache) = min\_data\_points$ then

12: $\mathcal{E}$.train(cache)

13: end if

14: else

15: if $\text{IsEmulatorAccurate}(\mathcal{E}, \theta_i) = \text{True}$ then

16: $\tilde{x}_i \leftarrow \mathcal{E}(\theta_i)$

17: $l_i \leftarrow \mathcal{L}(\tilde{x}_i)$

18: else

19: $x_i \leftarrow \mathcal{T}(\theta_i)$

20: $l_i \leftarrow \mathcal{L}(x_i)$

21: if $l_i > l_{\text{threshold}}$ then

22: cache.append($x_i$)

23: $l_{\text{threshold}} \leftarrow \text{FindLikelihoodThreshold}(\text{cache})$

24: $\mathcal{E}$.train(cache)

25: end if

26: end if

27: end if

28: end for

## 2.2 Emulator design

For our active sampling strategy, we require an emulator that can estimate its own accuracy. Similarly to [59], we find that *Gaussian Processes* (GPs) applied to a compressed version of the data can achieve this goal, while reaching excellent accuracy even with a moderate training set size.
### 2.2.1 Data compression

The observables we consider, such as the CMB angular power spectra and the matter power spectrum, are high-dimensional ($O(10^3) - O(10^4)$) but encode a relatively small number of features. To efficiently extract these features, we employ *Principal Component Analysis* (PCA) for data compression (see reference [60] for more details on PCA decomposition).

Our compression scheme transforms the data linearly, mapping it from the highly correlated and high-dimensional *data space* into the uncorrelated *feature space* of PCA components. We find that between 10 and 40 orthogonal PCA components are usually sufficient to represent the cosmological observables considered in this work. More complex cosmological models, which involve a larger number of parameters, typically require more PCA components. In practice, the optimal number of components is determined in OLÉ based on its precision parameters, ensuring minimal information loss during the compression of observables. Since the features are computed from the data in the cache itself, with a growing cache the PCA components will eventually change. This results in a better resolution of subdominant physical effects and a better attribution of relevance between the found features. However, this change of PCA components is negligible for individual data points. Therefore, we keep the PCA components constant for most of the added points and then update them when a certain number of new samples is accumulated.<fnref n="5" /> With OLÉ's important attribute of removing outliers (see section 2.1), and since PCAs are sensitive to outliers, this strategy ensures that the compressed PCA components are the relevant features of the data.

An ordinary PCA is agnostic to the experimental sensitivity and assigns equal weight to all features. However, since we have some prior knowledge of the experiment, we can optimize the compression by normalizing the data using the observational errors.<fnref n="6" /> We find that this normalization reduces the number of required PCA components while improving the accuracy of the most significant features in the data.

### 2.2.2 Gaussian processes

After transforming to feature space, we exploit the orthogonality property of the PCA components to construct a collection of independent one-dimensional emulators, each modeled using GPs, summarised in appendix D. This approach is highly beneficial since multidimensional GPs can be computationally expensive.<fnref n="7" />

The first few PCA components, representing the most prominent features of the data, often exhibit a relatively simple dependence on the input parameters $\theta = \{\theta_1, \dots, \theta_d\}$, typically following linear or quadratic relationships. In contrast, the high-order PCA components tend to encode a mixture of subdominant physical effects or even numerical noise induced by the theory code, leading to more complex or even random dependency on the input parameters.

<footnote n="5">⁵This strategy saves computational resources since a set of new PCA components requires an entire training of the emulator, while the addition of a single data point necessitates only a computationally cheaper update of the emulator.</footnote>

<footnote n="6">⁶When no observational errors are provided to OLÉ, we assume them to be proportional to the sampled variance of the training set. This is equivalent to not prioritizing any specific part of the observable.</footnote>

<footnote n="7">⁷Multidimensional GPs scale quadratically with the emulated dimension because all cross-correlations must be inferred. In our case, since PCA components are approximately uncorrelated, the complexity is reduced from $O(N^2)$ to $O(N)$.</footnote>
Every Gaussian Process relies on a choice of kernel and mean function, which must be optimized to accurately capture the response of a given PCA component to each parameter in the $\theta$ basis. To simplify the model, we normalize the data such that a zero (vanishing) mean function can be assumed — an approach that is commonly adopted when working with Gaussian Processes. Hence, we select the prior of the Gaussian Process to be zero. Thus, without any training data the GP will return zero with a variance of one.

To efficiently describe the dependence of the PCA components on the input parameters, we employ a composite kernel comprising the sum of three components: an *anisotropic radial basis function* (RBF) kernel, a *linear kernel*, and a *white-noise kernel*. This combination enables the emulator to better distinguish between physically meaningful variations and numerical noise, improving the robustness of the model.

The anisotropic RBF kernel allows the emulator to describe smooth variations of the PCA components in the vicinity of their distribution. Its hyperparameters consist of one length scale per input dimension (the dimensionality being $d$, the number of cosmological parameters). These length scales are fitted to the training set. They can be interpreted as the range over which a support point provides an accurate prediction, as neighboring points within this scale remain highly correlated. If no support point lies within this scale, the RBF kernel pushes the emulator to revert to the mean value of the PCA components while assigning a large emulation error. This typically prompts the addition of a new support point to the data cache, enhancing the emulation accuracy in future emulator calls.

The linear kernel instead captures the overall linear dependence of the PCA components on the input parameters and, in combination with the RBF kernel, enables the response of the PCA components to each parameter to be tilted overall. The linear kernel adds one hyperparameter per input dimension, which quantifies the amount of linear correlation that is fitted to the data. It is particularly beneficial for emulation far from the best-fit point in regions without nearby support points, where the linear kernel yields a suppression of the likelihood alongside a large prediction error. By incorporating the linear kernel, we can effectively disregard points that are too distant from the best-fit, improving the robustness of the emulator and preventing unnecessary additions to the support points far from the central emulation region, as long as the prediction error remains consistent with our precision goals.

The white noise kernel adds one single hyperparameter accounting for the level of added white noise. We do not fit it to the data, but fix it according to the information loss incurred during data compression.<fnref n="8" /> This ensures that the inclusion of the noise kernel does not substantially inflate our overall noise budget. By setting an upper bound on the achievable accuracy of the GPs, the noise kernel prevents them from overfitting minor artifacts that have negligible practical impact on the final results. This leads to a more robust and efficient hyperparameter estimation. By refining the precision settings of the data compression process, this source of uncertainty is also automatically reduced.

<footnote n="8">⁸We compute the variance of the residuals introduced by the truncation of PCA components above a given order. We use this variance to estimate the overall white noise level. We distribute this white noise evenly between all PCA components, such that its impact is negligible for the dominant PCA components and significant for the noise-dominated ones.</footnote>
### 2.2.3 OLÉ accuracy estimation

GPs offer a significant advantage in that they not only predict the expected value but also quantify their own uncertainty. This capability allows us to systematically assess whether the emulator is reliable within a given region of parameter space or if additional training points are still required. When assessing the error of an emulated observable $\tilde{x}$ at a given sample $\theta$, it is relevant to account for its location in parameter space and, consequently, its likelihood $l$. While an error of a few percent may be acceptable far away from the global maximum likelihood, higher accuracy is required close to the best-fit point. Therefore, rather than imposing a uniform accuracy requirement for the GPs, we adapt an error threshold based on the likelihood value. This approach provides a natural interpretation that directly links the emulator's accuracy to the precision demands of the parameter inference task, ensuring that higher accuracy is prioritized in regions where it has the most impact on the final results.

Furthermore, the required precision of an emulated observable depends on its inherent measurement accuracy. An observable measurable at per-mille precision imposes stricter emulation accuracy requirements compared to one that is only determined at order unity.

**Algorithm 2 OLÉ algorithm for the IsEmulatorAccurate function.**

**Input:** (sample) $\theta$

**Output:** (bool) accurate_flag

**Require:** (OLÉ emulator) $\mathcal{E}$

**Require:** (likelihood) $\mathcal{L}$

**Require:** (function) PrecisionCriterion

**Require:** (int) $N_{\text{quality\_samples}}$ (OLÉ hyperparameter)

1: $l_{\text{samples}} \leftarrow \text{list()}$

2: **for** $i \leftarrow 1, N_{\text{quality\_samples}}$ **do**

3:     $\tilde{x}_{\text{PCA,sample}} \leftarrow \mathcal{E}.\text{sample\_from\_GP}(\theta)$

4:     $\tilde{x}_{\text{sample}} \leftarrow \mathcal{E}.\text{decompress}(\tilde{x}_{\text{PCA,sample}})$

5:     $\tilde{l}_{\text{samples}}.\text{append}(\mathcal{L}(\tilde{x}_{\text{sample}}))$

6: **end for**

7:     $\tilde{x} \leftarrow \mathcal{E}(\theta)$

8:     $\tilde{l} \leftarrow \mathcal{L}(\tilde{x})$

9:     $\Delta_{\log \tilde{l}} \leftarrow \sqrt{\text{median}((\log \tilde{l}_{\text{samples}} - \log \tilde{l})^2)}$

10:    $\Delta_{\log \tilde{l}, \text{max}} \leftarrow \text{PrecisionCriterion}(\log \tilde{l})$

11: **if** $\Delta_{\log \tilde{l}} < \Delta_{\log \tilde{l}, \text{max}}$ **then**

12:     accurate_flag $\leftarrow$ True

13: **else**

14:     accurate_flag $\leftarrow$ False

15: **end if**

Eventually, we estimate an error on the likelihood by sampling from the GP emulator and computing the response of the likelihood function (see the pseudocode of the *IsEmulatorAccurate* Algorithm 2). This procedure samples a set of observables from the error estimate of the GPs that are computed for a given point $\theta_i$. By evaluating the likelihood
for each of these plausible observables, we obtain a distribution of likelihood values, which allows us to quantify the uncertainty of the likelihood estimation.

Instead of directly measuring the standard deviation from this collection, which is susceptible to large sample variance, we construct a more robust statistic by sampling the GPs at their one-sigma boundary. We then compare the median difference of these samples to the likelihood obtained from the mean prediction of the emulator. This approach provides a more stable and reliable measure of the uncertainty in the likelihood estimation, reducing sensitivity to outliers and improving the robustness of our accuracy assessment.

In order to combine the errors of our individual GP emulators into one final error for the likelihood, we assume that all emulators are entirely uncorrelated, despite relying on the same support points, and that they contribute independently to the likelihood evaluation. In practice, these assumptions are often not fulfilled, leading to an overly optimistic estimation of the uncertainty. Nevertheless, these effects primarily introduce a constant scaling factor in our error estimates. This scaling factor can be empirically determined during the deployment of the emulator by comparing the estimated errors with the direct results obtained from the theory code. This calibration step allows us to correct for the overestimation of precision and improve the reliability of our uncertainty quantification.

In the following we look at examples of the scaled error estimate. In figure 1 we compare the difference between the (log-)likelihood of samples that were computed using the emulator ($\tilde{l}$) and the (log-)likelihood of samples that used the theory code ($l$). All samples were randomly selected from the Markov chain Monte Carlo (MCMC) analysis of Example 1 in section 3.1. In figure 1(a), one observes the distribution of the samples' likelihood and the distribution of the Hubble parameter. The largest likelihood can be obtained in the center of the Hubble parameter distribution. This supports our expectation that the best-fit point is located near the center of our MCMC samples of single-peaked, (approximately) symmetric distributions. The difference between the simulation-inferred likelihood and the emulated one is encoded by the color. The largest deviations can be found at the smallest likelihoods, where their influence is smaller. Overall, we observe no systematic over- or underestimation which would directly lead to systematic biases in the estimate of the posterior. In figure 1(b) we show a histogram of the differences in the (log-)likelihoods scaled by the error estimator. As before, the histogram is centered at zero which indicates that there is no overall bias. The standard deviation is about unity, which indicates that the error estimate is accurate.<fnref n="9" />

Having validated the statistic of the error estimate of OLÉ, we still need to construct a validity criterion that dictates whether the OLÉ prediction, given its uncertainty, will or will not be used. In Algorithm 1, we denote this function `IsEmulatorAccurate`. We choose a validity criterion such that the error estimate $\Delta_{\log \tilde{l}}$ must be smaller than the maximum acceptable error $\Delta_{\log \tilde{l}, \max}$ for the accuracy to be deemed acceptable. This upper bound consists of three terms

$$ \Delta_{\log \tilde{l}, \max} = c_0 + c_1 \log \frac{\tilde{l}}{l_{bf}} + c_2 \log^2 \frac{\tilde{l}}{l_{bf}}, \qquad (2.7) $$

<footnote n="9">⁹We find that the distribution is not necessarily Gaussian but exhibits wider tails on both sides. However, this does not change our interpretation.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Evaluation of the log-likelihoods computed with predictions from CLASS (log *l*) and OLÉ (log *l̃*) in the example of section 3.1. On the left they are displayed as function of the Hubble rate *h*. The over or undershooting of the estimated likelihood from using the OLÉ prediction is indicated by color. On the right the difference between the likelihoods are scaled by the error estimate of the emulator $\Delta_{\log l̃}$. The solid and dashed black lines indicate the one standard deviation region and mean, respectively. The likelihood values calculated using the emulator prediction are unbiased and the emulator's uncertainty estimation is accurate.</figcaption></figure>

that depend on the three user-given precision parameters $c_i$ for $i \in (0, 1, 2)$, the estimated likelihood of the emulator $\tilde{l}$ and the likelihood of the current best-fit that was found by the sampler $l_{bf}$. The first parameter $c_0$ represents the acceptable ground level error around the maximum of the distribution. Note that this parameter should be chosen as a function of the dimension of the sampled distribution.<fnref n="10" /> The other two terms serve to degrade the required accuracy at the edges of the distribution.

For Gaussian posterior distributions, the linear precision parameter $c_1$ has an intuitive interpretation. In this case, the log-likelihood is quadratic in the sampled parameters $\theta$. Hence, when inferring the relation between the log-likelihood and $\theta$, the first derivative is linear in $\theta$. Thus, the credible intervals in $\theta$ (e.g. around the $1\sigma$ boundary) are related to credible intervals in $\log l/l_{bf}$, such that a change in $\log l/l_{bf}$ will be directly proportional to a change in $\theta$. For the relative errors, one finds that $\frac{\Delta \log l}{\log l/l_{bf}} = 2\frac{\Delta \theta}{\theta}$. In the case where equation (2.7) is dominated by the linear term one can directly relate $c_1$ to $2\frac{\Delta \theta}{\theta}$. Thus, in order to determine the position of the $1\sigma$ contour to a percent level, one should set $c_1 = 0.02$.

The quadratic term can be seen as a cutoff scale at which the demand on accuracy falls off. The default value of this term is set such that it dominates equation (2.7) at log-likelihoods outside of the 3–4 $\sigma$ parameter region at which little calls are expected anyways.

<footnote n="10">¹⁰An error of 0.1 log-likelihood is very relevant for a 1 dimensional distribution while it might not be for a higher dimensional one. For example, the difference between the best-fit and the $1\sigma$ bound for 1-dimensional case is at 0.5 log-likelihood for the 1-dimensional case, while in 27 dimensions it is at ~14.9.</footnote>
## 2.3 Code design

The code can be interfaced with the cosmological samplers MontePython<fnref n="11" /> [61] and Cobaya<fnref n="12" /> [62]. Therefore, OLÉ is designed to be parallelized with MPI by constructing an emulator for each MPI-process, each corresponding to an MCMC chain. In order to increase efficiency, all chains share one cache such that simulation calls can be shared among all processes. This significantly reduces the duration of the initial data collection to train the emulator. Furthermore, the training of the emulator is parallelized and MPI communication is implemented such that updates to one process are forwarded to all other processes, making all emulators based upon the same set of support points.

A major acceleration of the emulator can be achieved by using the JAX software package [63], which is specialized on high-performance numerical computing. The just-in-time (jit) compilation of JAX transforms the functions of the emulator and compiles it during the runtime of the code. This enables the emulator to be executed efficiently, accelerating it by a factor of ~ O(10²). However, this compilation also results in an increase in memory consumption compared to using just the standard theory code.<fnref n="13" />

JAX also exposes OLÉ to an automatic differentiation algorithm. This makes it easy to obtain robust derivatives of cosmological observables with respect to parameters. Using a series of examples, references [64–66] illustrate how access to such derivatives facilitates faster and more efficient analysis pipelines; we add to this with the development of a staged gradient-based sampler in section 4.

Lastly, for the implementation of the GPs, we use the python package GPJax [67] that comes with an efficient JAX-compatible implementation for the construction, training and evaluation of GPs.

## 2.4 OLÉ in MCMC analyses

Much of modern cosmological inference is done via MCMC analyses with Metropolis Hastings (MH) sampling, which requires a large number of calls to a cosmological theory code to compute the observables for a given set of input parameters. Using OLÉ greatly reduces the time spent on this aspect, speeding up the acquisition of samples, eventually reaching a level of speed-up where the runtime of the likelihood code is the dominant step.

During the initial burn-in phase, the total runtime is dominated by evaluating the theory code in order to burn-in and accumulate training data. In this phase, we use oversampling (keeping the ‘slow’ cosmological parameters fixed during several steps while sampling over ‘fast’ nuisance parameters<fnref n="14" />) to converge more quickly to the center of the posterior distribution. Once the training is done, OLÉ proceeds with the test phase in which the execution time is largely taken up with testing/training and updating the emulator. However, the testing rate

<footnote n="11">¹¹https://github.com/brinckmann/montepython_public.</footnote>

<footnote n="12">¹²https://github.com/Cobaya Sampler/cobaya.</footnote>

<footnote n="13">¹³An example MCMC with 12 parallel MPI processes requires ~ 5 GB of RAM with standard CLASS and 35-80 GB with OLÉ, depending on precision settings. The biggest memory requirement comes from the ‘jitted’ emulator, as each MPI process needs to keep this in its memory. The more points in the cache, the more memory each process will use. The total memory usage can thus be reduced by running fewer parallel processes, but in practice this should rarely be a problem on a modern computing cluster.</footnote>

<footnote n="14">¹⁴See reference [62] for details.</footnote>
is one of the limiting factors in the achievable acceleration of the MCMC. As this necessitates multiple calls to the likelihood code, too frequent testing can quickly become a computational bottleneck. Therefore, the emulator gradually lowers the testing rate until it reaches a user-defined minimum. Eventually, once the emulator is confident and the testing rate is reduced, the runtime is dominated by the evaluation of the likelihood rather than the theory code evaluation. When OLÉ reaches a level of accuracy where only a negligible part of the total MCMC runtime is spent evaluating the theory code, we can derive an upper bound on the total acceleration induced by the emulator. Assuming that the evaluation time of the emulator ($T_E$) is much smaller than the one of the theory ($T_T$) and likelihood ($T_L$) codes, we get:

$$ \textit{Speed-up} = \frac{T_T + T_L}{T_L}. \qquad (2.8) $$

This speed-up is generally achievable for a variety of likelihood codes and cosmologies, as we present here, and is particularly evident in figure 2. Another consequence of the likelihood code dominating the runtime is that oversampling is no longer useful. It even slows down the acquisition rate of points with new cosmological parameters in the chains. Oversampling is only beneficial when the emulator remains slower than the likelihood code or during the initial burn-in phase, when the emulator is rarely used. Hence, in this stage of the MCMC run, OLÉ automatically switches off the fast/slow parameter decomposition.

The progress of an example chain, and with it of the emulator, is displayed in figure 3, using a chain from our standard ΛCDM example (see section 3.1). We distinguish between the initial burn-in phase of about 300 calls to CLASS during which the chain moves towards the center of the posterior and the acquisition of the first data points used to train the emulator. While about 2/3 of those calls are regarded as outliers of little use to the emulator in the latter stage (which is dominated by calls in the most central region of the parameter space), about 1/3 of the initial calls remain in the training set. Once ~ 80 training points are gathered, the emulator is trained for the first time. In the early phase after the initial training (and in particular after new points are added to the cache), most emulator calls are tested to ensure a sufficient accuracy of the emulated quantities. During this phase, “holes” in the training set that lead to inaccurate predictions are detected and filled by additional calls to the theory code. Eventually, most gaps are closed and the emulator is confident with its predictions (achieving the user-defined desired accuracy level).

The most important OLÉ precision parameters are presented in appendix B.1, along with their default values.

## 3 Examples with Metropolis Hastings sampling

We demonstrate the accuracy and efficiency of OLÉ as a tool for parameter inference through a selection of examples covering a variety of samplers, theory codes, cosmological models and observational likelihoods. In this section, we limit ourselves to the use of the MH algorithm as implemented in the Cobaya or MontePython samplers, and to the theory codes CLASS and CAMB. The next section features additional sampling algorithms and parameter inference packages.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Fraction of elapsed computation time spent on each component of the process, spanning the first ~1.5 hours of an MCMC analysis with OLÉ. It is clear that the initial phase is dominated by running the traditional theory code until enough points have been computed to train the emulator. Then comes a phase dominated by testing the emulator, as its accuracy must be verified. Finally, when confidence in the emulator is high, the runtime becomes dominated by the likelihood. As the run progresses, it will become even more likelihood-dominated.</figcaption></figure>

**Quantifying efficiency.** When running Cobaya, we sample until reaching a desired value of the convergence criterion,<fnref n="15" /> while with MontePython we run our chains until hitting a given number of points or a maximum wall-clock time. To quantify the efficiency of OLÉ versus traditional runs, we could simply compare the time it takes on a given platform to reach a given convergence level or a given number of points. Such a test would not be robust and informative enough, because such a time could be strongly affected by arbitrary choices for the number of chains running in parallel, the number of cores per chain, the oversampling factor for nuisance parameters, the initial proposal density, etc. In order to have a uniform metric, we refer to the effective sample size (ESS) of the cosmological parameters of MCMC chains (see reference [69]). The ESS provides an estimate of the effective number of uncorrelated samples in a chain. It is obtained by first estimating the scale on which the samples of a chain are correlated, and then dividing the number of total samples by this correlation scale.<fnref n="16" /> Note that oversampling allows to acquire samples at a higher rate in traditional runs at the expense of increasing the correlation scale, such that each drawn (cosmological) sample contains less information. Thus, using the ESS instead of the number of samples allows us to have a meaningful comparison between the performance of runs performed with or without oversampling. This is useful in our context since, as mentioned in section 2.4, OLÉ switches off oversampling when no longer useful. The CPU time needed on a given platform to reach a given ESS directly estimates the efficiency of a MH sampling pipeline, since it quantifies the speed at which relevant information accumulates. For the CPU time, we should use the effective amount of time during which the CPUs were active,

<footnote n="15">¹⁵We use the Gelman-Rubin statistic $|R - 1|$ [68].</footnote>

<footnote n="16">¹⁶In this work, we compute the ESS using `GetDist` [70] and `ArviZ` [71].</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Value of the reduced Hubble parameter $h$ for the first 15000 theory evaluations in one of the chains of a Cobaya MCMC run, as well as 6000 evaluations later in the chain for comparison. The points are colored according to the type of call to the standard theory code or emulator: CLASS calls saved in the cache and used to improve the emulator (black cross), CLASS calls deemed not relevant for the emulator and therefore not cached (gray dots), emulator calls followed by a successful accuracy test (green dots), emulator calls followed by a failed accuracy test (red dots), and emulator calls not followed by an accuracy test (blue dots). Gray dots appear only at the beginning during the burn-in phase. Note that red dots (failed accuracy tests) are immediately followed by black dots (CLASS calls added to the cache). As performance checks are passed and confidence in the emulator grows, testing is gradually decreased, such that green dots give way to blue dots. Note that points very far from the high- likelihood region would always be computed by the standard theory code to avoid an overly expensive and conservative training of the emulator, but this situation did not occur within the displayed sample.</figcaption></figure>

which reflects the energy consumption of the runs. This time can be obtained by multiplying the wall-clock time of a run by the number of chains, the number of cores per chain, and the average CPU usage efficiency reported by the computer at the end of the run. We denote the effective CPU time in hours as CPUh. In summary, we will compare the efficiency of different runs performed with or without OLÉ by looking at the ratio of their ESS per effective CPU time in hours, ESS/CPUh.

**Quantifying accuracy.** Various estimators can be defined to quantify the agreement between parameter inference runs based on OLÉ and traditional runs. For simplicity and concision, we will focus on the shift $D_x$ of the mean value of the 1-dimensional posterior for each parameter $x$ normalized to the standard deviation found in the traditional run, $\sigma_x$. This dimensionless shift is given as

$$D_x = \frac{|\langle x \rangle - \langle \tilde{x} \rangle|}{\sigma_x}, \qquad (3.1)$$

where $\langle \tilde{x} \rangle$ refers to the mean found in the OLÉ-based run, and $\langle x \rangle$ to that of the traditional run. Even when comparing two well-converged runs, sampling noise usually leads to shifts
$D_x$ in the range from 0.01 to 0.1 [53], indicating that the mean value is estimated with good accuracy compared to the parameter uncertainty. In the following sections, we will quote the relative shifts $D_x$ obtained for various parameters, assuming different cosmological models and likelihoods. Values of $D_x$ smaller than 0.1 will indicate that the **OLÉ** emulator is not significantly biasing the posteriors.

## 3.1 ΛCDM with current data using Cobaya and CLASS

**Model.** As a first example, we investigate the performance of the emulator by running Cobaya for the standard ΛCDM model, both with and without the use of OLÉ. We parameterize the ΛCDM model via: the expansion rate today $h$ in units of 100 km/s/Mpc, the fractional density of baryons $\Omega_b$ and cold dark matter $\Omega_{cdm}$, the amplitude $\ln(10^{10}A_s)$ and spectral index $n_s$ of the primordial power spectrum of scalar fluctuations, as well as the optical depth to reionization $\tau_{reio}$. One of the three neutrino species is assumed to have a fixed mass of 0.06 eV. A definition of all cosmological parameters referenced in this work can be found in table 4 of appendix C.

**Data.** We include Planck 2018 data for temperature, polarization [4] and lensing [72] anisotropies, as well as Supernova luminosity data from the Pantheon sample [73] and Baryon Acoustic Oscillation (BAO) data from the SDSS DR7 [74] and DR12 [75] samples. In both types of runs, we use CLASS as the theory code to compute observables, and the MH sampler of Cobaya to run 8 parallel chains up to a Gelman-Rubin statistic of $|R-1| < 0.01$. The likelihood codes used here are not as highly parallelized as the theory codes. Since the OLÉ runtime is dominated by likelihood evaluations, it is not extremely efficient to run chains on many cores in the OLÉ case. A higher number of cores would speed up the initial phase during which the theory code is often called, but then, once the run time gets dominated by non-parallelized-likelihood evaluations, it would result in an inefficient use of resources with no benefit. We thus run each of the 8 chains on 4 cores each in the traditional run, and on a single core each in the OLÉ run. We also mentioned in section 2.4 that oversampling is no longer useful with OLÉ. In the traditional run, we use the default Cobaya settings for oversampling the 21 nuisance parameters, but in the OLÉ run, the fast/slow parameter decomposition gets switched off automatically after the initial training phase. We perform the OLÉ run using the default values of the OLÉ precision parameters described in appendix B.1.

**Accuracy.** We find excellent agreement between the posteriors derived from the two runs, as can be seen in figure 4. The normalized shifts of the means $D_x$ of equation (3.1) are found to be below 0.03. Such tiny deviations are consistent with the stochastic variance of the posteriors estimated from MCMCs [53].

**Efficiency.** In order to estimate the maximum achievable speed-up according to equation (2.8), we measure $T_\tau \sim 5$ sec and $T_\ell \sim 0.05$ sec on 4 cores in this particular case. Consequently, we may expect an acceleration by at most about a factor of 100. Note that this is only a rough estimate, which does not take into account several effects such as the initial burn-in phase, the computation efficiency, the cost of calling CLASS to build the cache, that of testing the emulator accuracy, or the impact of oversampling over nuisance
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Posteriors of the ΛCDM model parameters inferred from CLASS and Cobaya with Planck 2018, Pantheon+ and BAO data, as outlined in section 3.1. We see an excellent match between the posteriors computed using either a standard MCMC or the OLÉ emulator, with a maximum deviation of posterior means by 0.03σ.</figcaption></figure>

parameters.<fnref n="17" /> In practice, by comparing the sampling rate ESS/CPUh in the two runs, we find that OLÉ accelerates the whole parameter inference pipeline by a factor ~98, which is very close to the expected speed-up.

We conclude that OLÉ is able to boost the efficiency of MCMC runs by a considerable amount without degrading their accuracy. However, it is worth checking that this conclusion applies also to more complex models and other likelihoods.

<footnote n="17">¹⁷Even if, as explained in section 3, using the ESS instead of the number of points in the chains mitigates the impact of using oversampling.</footnote>
## 3.2 Extended cosmology with current data using Cobaya and CAMB

**Model.** To investigate the performance of OLÉ in the case of computationally more expensive models, we extend the ΛCDM cosmology to include a non-zero spatial curvature parametrised with the effective density Ω$_k$, a (degenerate) neutrino mass leading to the fractional neutrino density Ω$_ν$, and a dark energy (DE) fluid with an equation of state $w(a)$ obeying the Chevallier-Polarski-Linder (CPL) parametrization [76, 77]

$$w(a) = w_0 + w_a(1 - a), \qquad (3.2)$$

with free parameters $w_0$ and $w_a$. For the minimal ΛCDM parameters, to illustrate the versatility of the code, we use the same parameters as in the previous subsection except for the total non-relativistic matter density Ω$_m$ (replacing Ω$_{cdm}$) and the current matter power spectrum amplitude parameter $\sigma_8$ (replacing $A_s$). Note that for the OLÉ emulator, the effect of the neutrino mass (or equivalently of Ω$_ν$) on the CMB spectra is potentially more challenging to capture than that of other parameters. As a matter of fact, for very small values of Ω$_ν$, the response of the CMB spectrum coefficients $C_ℓ$ to variations of Ω$_ν$ is far from linear.

**Data.** We base this analysis on Planck 2018 temperature, polarization and lensing measurements, in conjunction with DESI 2024 BAO data [3]. The observables are computed with CAMB using the Parametrized Post-Friedmann (PPF) DE parametrization [78]. We use again the MH sampler of Cobaya. As in the previous example, we carry out two runs, a traditional one with 8 chains using 4 cores each, and an OLÉ run with 8 chains using a single core each. Both runs proceed until reaching a Gelman-Rubin statistic of $|R - 1| < 0.02$. Like in the previous subsection, we oversample the 21 nuisance parameters only in the traditional run and in the initial training phase of the OLÉ run. As explained in appendix B.2, we perform the OLÉ run using default precision, except for two parameters. First, for a faster run, we loosen our requirement on the log-likelihood error $\Delta \ln \mathcal{L}$ below which the emulator is considered valid. However, to get an accurate emulation of the physical effect of the neutrino mass, we enhance the accuracy setting controlling the number of PCA components.

**Accuracy.** Figure 5 shows an extremely good level of agreement between the posterior distributions obtained in the two cases. In particular, the non-Gaussian nature of the posteriors for additional parameters is well captured by both MCMCs. The normalized shifts of the means $D_x$ between the posteriors computed with or without OLÉ never exceeds 0.048, which is comparable to the results of the previous subsection. Note that the physical observables of the extended model feature four additional physical effects compared to the ΛCDM case. We find that the OLÉ emulator is able to capture these additional effects and to maintain the same final level of accuracy on parameter inference. This illustrates the robustness and the stability of the automatic error estimation algorithm implemented in OLÉ.

**Efficiency.** In this case, the runtime of CAMB on 4 cores is ~15 seconds,<fnref n="18" /> while the evaluation of the likelihoods takes about ~0.05 seconds. Thus, equation (2.8) suggests a maximum achievable speed-up of about ~300, again, neglecting the cost of building up the cache and

<footnote n="18">¹⁸This runtime is increased because we ask for high precision in the CMB spectra at high $\ell$, as required by the SPT likelihood. Note that using lower accuracy settings could result in at most 4 points shift in $\chi^2$ (see appendix A.3 of [79]).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Posteriors of the extended cosmology model on *Planck* and 2024 DESI BAO data as outlined in example 3.2. Posterior means obtained using **OLÉ** and **CAMB** differ by $D_x < 0.048$, indicating good emulator accuracy.</figcaption></figure>

testing the emulator accuracy in the **OLÉ** run. In practice, we observe that **OLÉ** increases the overall sampling rate from 0.5 ESS/CPUh to 173.1 ESS/CPUh, boosting the entire MCMC run by a factor of 354. This is particularly impressive considering the fact that the **OLÉ** emulator does not require any pre-training.

**Adding nuisance parameters.** To test **OLÉ**'s performance in even more extreme conditions, we stick to the same extended model, but we add CMB temperature and polarization data from SPT-3G 2018 [5] to that from *Planck* and 2024 DESI BAO. The SPT-3G 2018 likelihood includes several nuisance parameters, resulting in a final 64-dimensional parameter space (10 cosmological parameters, 21 *Planck* nuisance parameters, and 33 SPT-3G nuisance parameters). With such high dimensionality, MCMC runs become prohibitively expensive.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Posterior distributions of the extended cosmology model on Planck, SPT-3G, and 2024 DESI BAO data as outlined in example 3.2. Posterior means obtained using OLÉ and CAMB differ by $D_x < 0.02$, indicating excellent emulator accuracy.</figcaption></figure>

Using CAMB, the MH sampler of Cobaya and its default oversampling option, we find that by running 8 chains (on 4 cores each) we reach $|R-1| = 0.02$ after 13 days and 11 hours. Instead, an equivalent run performed using OLÉ (with only 2 cores per chain and no oversampling) converges at the level of $|R-1| = 0.02$ in just 39 hours.<fnref n="19" /> The accuracy of the OLÉ run remains excellent, with $D_x < 0.02$.

The examples of this section suggest that the accuracy and efficiency of OLÉ scale extremely well with the complexity of the physical model and with the number of nuisance parameters. Some runs that would have a considerable computational cost with a traditional MCMC approach become easily feasible using the on-the-fly OLÉ emulator, without compromising accuracy.

<footnote n="19">¹⁹We don’t report the effective sampling rate in this case, because the CPU efficiency was not reported by the HPC cluster.</footnote>
## 3.3 Stage-IV LSS forecast for an extended cosmology using MontePython and CLASS

**Data.** As a further test of the emulator’s applicability to different types of likelihoods and data, we run MCMC chains on synthetic data accounting for a Stage-IV galaxy survey, using a version of the mock Euclid likelihoods for MontePython described in references [80–83]. More precisely, like in reference [83], we use a set of likelihoods describing the power spectrum of weak lensing maps, photometric galaxy clustering maps, the cross-correlation between the two previous observables, and spectroscopic galaxy clustering maps (abbreviated as WL + GC<sub>ph</sub> + XC<sub>ph</sub> + GC<sub>sp</sub>).<sup>20</sup> Both of the runs of this section are based on the MH sampler implemented in MontePython, with 12 parallel chains running for 65 hours. We run the standard CLASS chains with 8 cores per chain oversampling the nuisance parameters by a factor of five and the OLÉ chains on a single core per chain. We perform the OLÉ run using precision parameters very similar to those of the previous subsection, with a looser-than-default requirement on the log-likelihood error $\Delta \ln L$ below which the emulator is considered valid, and an enhanced accuracy setting for the number of PCA components (see appendix B.3 for details).

**Model.** Like in the analysis of reference [83], we use an extended cosmology with four free parameters on top of the five standard $\Lambda$CDM parameters.<sup>21</sup> Our model includes the CPL parametrization of DE, see equation (3.2), massive neutrinos parameterized with the summed mass $\sum m_\nu$, as well as additional ultra-relativistic free-streaming relics parameterized through $\Delta N_{\text{eff}}$. We adopt fiducial parameter values close to the Planck $\Lambda$CDM best-fit, with a pure cosmological constant, $\sum m_\nu = 0.06$ eV, and no relativistic relics.

**Accuracy.** Visually, we find excellent agreement between the CLASS-only posteriors and those computed with the help of OLÉ, with normalized shifts in the means, $D_x \lesssim 0.1$ for all parameters, with $D_x \sim 0.01$ for the majority. A close inspection of figure 7 shows that in all our runs, there is a small shift between the mean posterior value and the fiducial value for $h$ and $n_s$. This is caused by the non-Gaussian nature of the posterior, as already found in reference [83] (see figure 13 therein).

**Efficiency.** In this case, we measure $T_\tau \sim 19$ s and $T_\ell \sim 0.32$ s on 4 cores. The maximum achievable speed-up according to equation (2.8) is around $\sim 59$, slightly smaller than in previous cases because the mock Euclid likelihoods are significantly slower than the Planck, supernovae and BAO likelihoods.<sup>22</sup> The ESS and the effective sampling rate ESS/CPUh

<footnote n="20"><sup>20</sup>We made some very minor changes to the existing public Euclid likelihoods, since unlike the default CMB and BAO likelihoods, these are not well-optimized for OLÉ’s automatic emulation. These changes are explained in appendix A.2, and the modified likelihoods will be included in a public MontePython release. Note that this does not change the fact that OLÉ works out of the box with both Cobaya and MontePython, requiring no modifications to existing likelihoods in nearly all cases.</footnote>

<footnote n="21"><sup>21</sup>In absence of CMB data, we do not include the optical depth to reionization in the list of free $\Lambda$CDM parameters.</footnote>

<footnote n="22"><sup>22</sup>CLASS is also significantly slower in this case because the Euclid likelihoods require a calculation of the matter power spectrum up to a large wavenumber and over a fine grid. However, in the ratio, the effect of likelihoods being slower wins.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Expected posteriors of the $w_0w_a$ CDM + $\sum m_\nu + N_{\text{eff}}$ model according to the mock Euclid WL + GC$_{ph}$ + XC$_{ph}$ + GC$_{sp}$ likelihoods described in section 3.3, with a fiducial model close to the Planck ΛCDM best-fit. Dashed gray lines indicate the fiducial parameter values. We recover the same posterior means with or without the OLÉ emulator to within $D_x \lesssim 0.1$. Details on the precision parameters can be found in the appendix B.3.</figcaption></figure>

found in runs with and without OLÉ are presented in table 1. They show that OLÉ speeds up the MCMC by a factor of 49, close to the theoretical speed-up factor. Once more, we find that OLÉ has a decisive impact on the time and energy consumption of the runs.

## 3.4 Early dark energy with current data using MontePython and TriggerCLASS

**Model.** As an example of a cosmology with strongly non-Gaussian posteriors, we use a model of New Early Dark Energy (NEDE), which has been studied recently in the context
<table><caption>**Table 1.** Mean effective sample size (ESS) across MCMC chains and effective sampling rate (ESS/CPUh) for our mock Euclid forecasts, in a traditional run and in an OLÉ run. For details on the OLÉ precision parameters see appendix B.3.</caption><thead><tr><td>Run</td><td>ESS (mean)</td><td>ESS/CPUh</td></tr></thead><tbody><tr><td>CLASS+MontePython</td><td>2204</td><td>0.35</td></tr><tr><td>same+OLÉ</td><td>13541</td><td>17</td></tr></tbody></table>

of cosmological tensions [84, 85]. This analysis relies on the theory code TriggerCLASS.<fnref n="23" /> In addition to the standard ΛCDM parameters, the NEDE model introduces four extra free parameters, {$f_{\text{NEDE}}, \Omega_{\phi}, \log_{10} z_{*}, w_{*}$}, accounting respectively for the fraction of NEDE at decay time, the density of the trigger field, the redshift at decay time, and the equation of state after the decay. We fix the summed neutrino mass to $\sum m_{\nu} = 0.06$ eV.

**Data.** Our analysis includes Planck 2018 data for temperature, polarization [4] and lensing [72] anisotropies, as well as Supernova luminosity data from the Pantheon+ sample [86] and BAO data from the SDSS DR7 [74] and DR12 [75] samples. We use the MH sampler of MontePython. We run 16 chains of 400,000 points each, on 4 cores each in the traditional run and 1 core each in the OLÉ run.

**Accuracy.** On top of experimenting the NEDE model, we use this section to test the impact of a few precision settings in OLÉ. In appendix B.4, we define two sets of OLÉ accuracy parameters that we call high precision (HP) and low precision (LP). These parameters control the criterion for considering the emulator as accurate enough in a given point and the number of PCA components. The LP settings degrade the default requirements on the emulator accuracy and on the number of PCA components. In the HP settings, the requirements on the emulator accuracy are stronger than in the default settings, but the number of PCA components is in between the default and LP settings. With both settings, we find excellent agreement between the posteriors obtained with and without OLÉ, as illustrated in figure 8. We do not compute the shifts in the means $D_x$ for this model, since normalizing to the standard deviation is less meaningful when the posterior shapes are so far from Gaussian. Instead we wish to stress the very high level of agreement of the posteriors, which is evident on visual inspection of the contours. This shows that the OLÉ emulator is able to train over a non-trivially-shaped region of the parameter space and to accurately capture all the subtle physical effects of the NEDE model. Additionally, the stability of the results against different accuracy settings (HP and LP) shows that our error auto-evaluation scheme is robust and that OLÉ does not require a fine-tuning of its precision parameters in order to remain accurate.

**Performance.** In this example, we find $T_T \sim 4$ s<fnref n="24" /> and $T_L \sim 0.08$ s on 4 cores. The running time of the likelihood is increased compared to section 3.1 because the Pantheon+ likelihood is slower than the older Pantheon one. The maximum achievable speed-up according to equation (2.8) is 50, smaller than in section 3.1 due to the slower likelihoods. The OLÉ

<footnote n="23">²³https://github.com/NEDE-Cosmo/TriggerCLASS.</footnote>

<footnote n="24">²⁴In section 3.1 we quoted a runtime of 5 s for CLASS. The small difference just comes from using slightly faster processors in the runs of this section.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Posteriors of NEDE parameters along with $\Omega_{cdm}$ and $h$. Despite the strongly non-Gaussian shapes of the posteriors, OLÉ does an excellent job of recovering the results of the standard MCMC analysis.</figcaption></figure>

<table><caption>**Table 2.** Mean effective sample size (ESS) across MCMC chains and effective sampling rate for our NEDE analysis, in a traditional run and in an OLÉ run using LP. For details on the OLÉ precision parameters see appendix B.4.</caption><thead><tr><td>Run</td><td>ESS (mean)</td><td>ESS/CPUh (norm.)</td></tr></thead><tbody><tr><td>CLASS+MontePython</td><td>24949.04</td><td>1.2</td></tr><tr><td>same+OLÉ (LP)</td><td>66670.37</td><td>31.2</td></tr></tbody></table>
run with LP is accelerated by a factor ~31. Even in this extreme case, OLÉ is able to save between one and two orders of magnitude in runtime and energy consumption while providing accurate results.

# 4 Staged NUTS sampling with differentiable likelihoods

In this section, we interface OLÉ with the CMB likelihood library candl.<fnref n="25" /> While OLÉ’s emulators calculate CMB spectra from a set of cosmological parameters, candl modifies the CMB spectra to account for foregrounds, calibration, and other systematic effects, and calculates the likelihood value for a chosen data set. Crucially, both codes are written with JAX support, which features an automatic differentiation algorithm; OLÉ and candl together form a fully differentiable pipeline with quick and easy access to derivatives of the log-likelihood value with respect to cosmological and nuisance parameters. This allows for the use of more efficient minimization and MCMC sampling algorithms [see e.g. 87–93].

We develop the following inference strategy to exploit the differentiability of OLÉ models and candl likelihoods:

1. Run MH sampling using a full Boltzmann solver with multiple chains in parallel, saving the predictions to the cache.

2. Once sufficient points have been collected, train and switch over to the emulator.

3. Minimize the posterior to find the best-fit point using the gradient-based truncated Newton algorithm implemented in Scipy [87, 89, 94].

4. Perform gradient-based No U-Turns (NUTS) [90] sampling of the posterior.

We refer to this prescription as *staged NUTS sampling* and visualize the procedure in figure 9. Note that during NUTS sampling we can still fall back to occasional evaluations of the Boltzmann solver and re-training of the emulator if its accuracy is insufficient. In fact, the testing of the emulator accuracy is improved when using a differentiable likelihood as the emulator uncertainty can simply be propagated through to the posterior value, without having to draw multiple realizations of the observables and evaluating the corresponding posterior values individually.

To assess the performance of staged NUTS sampling compared to regular MH sampling with OLÉ we compare the following two cases: (1) MH sampling with eight parallel chains with eight CPUs each using Cobaya, no oversampling of nuisance parameters and (2) staged NUTS sampling as described above across eight CPUs. In both cases we explore ΛCDM, using CLASS as the Boltzmann solver and the differentiable SPT-3G 2018 TT/TE/EE lite likelihood [5, 95] available in candl [65]. Both chains start with the same sub-optimal proposal matrix.

We find that NUTS sampling explores the posterior distribution more efficiently: the ESS/CPUh after the burn-in phase is a factor of 4.0 higher for our staged NUTS sampler than for the traditional MH chains. This indicates more informative samples are generated at a faster rate and intuitively means that staged NUTS chains with the same information

<footnote n="25">²⁵https://github.com/Lbalkenhol/candl.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Example chain using the staged NUTS sampler. We explore ΛCDM using the SPT-3G 2018 *TT/TE/EE* lite likelihood and show *h* samples. We indicate the phases of the staged sampler along the horizontal axis; from left to right: MH sampling (gray), minimization (orange), NUTS warmup sampling (light blue), NUTS sampling (blue). Points retained for emulator training have been highlighted with open black circles. Note that the x-axis scale is arbitrary and has been chosen for visualization purposes; we indicate the number of sample points in each stage along the top. This inference strategy exploits the differentiability of *OLÉ* and *cand1*; burn-in is effectively complete after the minimization step and the subsequent NUTS sampling explores the posterior more efficiently than traditional MH sampling.</figcaption></figure>

content as MH chains can be run in one quarter of the time, provided the same hardware and energy allocation. This is similar to what was found when comparing NUTS and MH sampling by [65] for CMB data from the Atacama Cosmology Telescope [2, 96] and by [97] for tomographic 3 × 2 point analyses. However, we stress that in these cases the emulators used were trained ahead of time using significant computational resources, whereas here we generate them on the fly with comparatively little cost.

Furthermore, we note that the time it takes to arrive at posterior-representative sampling, i.e. the burn-in, is efficiently used and elegantly defined for our staged NUTS sampler. Burn-in is completed with step (3) in the description above and since minimization is fast when gradient information is available, the burn-in time is effectively confined to how long it takes to gather training points and deploy the emulator. When this can be done depends on the complexity of the model being explored, though we recommend 80 points as a default option. This is lower than what is typical for MH-based samplers, though we note that in this case the burn-in definition varies depending on the specific implementation. For *Cobaya*, it is usually recommended to throw away the initial 20–30% of points, which scales unfavorably with the runtime due to the proportionality with the chain length.<fnref n="26" /> For *MontePython*, burn-in is regarded as complete once chains have reached the vicinity of the best-fit point, which is sensitive to the quality of the starting point and the initial proposal matrix as well as the shape of the posterior distribution.

<footnote n="26">²⁶See https://cobaya.readthedocs.io/en/stable/sampler_mcmc.html.</footnote>
# 5 Conclusion

In this work, we have presented the emulator code OLÉ, including its design goals and architecture. We have demonstrated its performance and accuracy, achieving comparable results to traditional Einstein-Boltzmann solvers CLASS and CAMB at a fraction of the computational cost across a range of cosmological models and datasets. OLÉ relies on active sampling and online learning to efficiently acquire training points with high information content on the fly, requiring no pre-training and outperforming similar emulators in data efficiency. It is released with its own built-in sampler as well as user-friendly interfaces for both the MontePython and Cobaya samplers, allowing for easy integration in existing inference pipelines. We showcased the performance on a variety of cosmological data sets including measurements from Planck, SPT-3G, Euclid mock data, BAO and supernovae data. The tested cosmological models range from the 6 parameter ΛCDM model over different extensions to a New Early Dark Energy model with 10 parameters with highly non-Gaussian posterior distributions. We find a computational speed-up in the parameter inference process ranging from one to two orders of magnitude, depending on the considered model and data set. It is close to the maximum achievable speed-up expected when considering the evaluation of CAMB or CLASS as the bottleneck of the inference process. Lastly, we use OLÉ together with the CMB likelihood library candl. We exploit the differentiability of both codes to construct a staged NUTS sampler that uses the gradient of the likelihood to explore the posterior distribution approximately four times more efficiently than traditional MH sampling. OLÉ is publicly available and we release a series of examples that show how to use the code.<fnref n="27" />

The application of machine-learning techniques is becoming wide-spread in cosmology and allows us, as demonstrated here, to accelerate common computational tasks. This facilitates the broadening and deepening of cosmological analyses by exploring constraints placed by many variations of the data across a broad range of models. OLÉ thus represents an advancement in numerical methods demanded by the precision of contemporary cosmological data. As the evaluation time of the theory code becomes subdominant, it is pertinent to improve the efficiency of our samplers and likelihood codes in the future. This development is particularly timely given the environmental impact of high-performance computing.

# Acknowledgments

This work uses JAX [63] and the scientific python stack [94, 101–103]. Triangle plots have been produced using GetDist [70]. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 101001897). This work has received funding from the Centre National d’Etudes Spatiales and has made use of the Infinity Cluster hosted by the Institut d’Astrophysique de Paris. RKS thanks the Alexander von Humboldt Foundation for their support. SG and JL acknowledge support from DFG grand LE 3742/6-1. JL and MM acknowledge support from the DFG grant LE 3742/8-1. Calculations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt. Computations were performed with computing resources granted by RWTH Aachen University

<footnote n="27">²⁷https://github.com/svenguenther/OLE.</footnote>
under project 'rwth1661'. We would like to thank Lasse Ausborm for testing OLÉ and Santiago Casas, Jesús Torrado, Andreas Nygaard, Emil Brinch Holm, Steen Hannestad, Will Handley, Thomas Tram, Marco Bonici, Nils Schöneberg and Silvia Galli for useful discussions and comments. SG is first author as the main author of the OLÉ code, remaining authors are in alphabetical order. The authors gratefully acknowledge the computing time provided to them at the NHR Center NHR4CES at RWTH Aachen University (project number p0021792). This is funded by the Federal Ministry of Education and Research, and the state governments participating on the basis of the resolutions of the GWK for national high performance computing at universities (https://www.nhr-verein.de/unsere-partner).

# A Using OLÉ

The OLÉ code has been designed to be easy to use, either as a full standalone pipeline, or as a supplement to accelerate the theory code in existing pipelines. For this purpose, OLÉ includes its own samplers, as well as interfaces for the Cobaya and MontePython samplers.

As OLÉ automatically copies the functionality of the theory code, it should work out of the box with any modified version of CLASS or CAMB that do not change the structure of calls to the wrapper required to compute the observables. For example, we make use of TriggerCLASS in section 3.4, which is one such modified version of CLASS.

## A.1 Implemented sampler

In addition to the emulator, OLÉ includes a set of integrated sampling methods. Alongside an ensemble sampler based on the Python package `emcee` [98], it features a minimizer utilizing Scipy [99] and a gradient-based NUTS sampler. An example demonstrating the combination of these methods is provided in section 4.

The implementation of theory codes is inspired by Cobaya. A theory code is implemented as a child class of the OLÉ-internal Theory class. It has an attribute requirements, which represents the set of quantities required by the likelihood. This attribute is dynamically determined by all likelihoods within the sampler.

The theory class also provides the function compute, which takes a dictionary of parameters and populates it with the required quantities. The likelihoods are implemented as child classes of the Likelihood class. Each likelihood child class must define the loglike function, which takes the dictionary of computed quantities and evaluates the log-likelihood. Additionally, it includes the function update_theorySettings, which allows enforcing specific settings onto the theory code. If the loglike function is implemented using JAX and is differentiable, gradient-based methods such as NUTS or a gradient-based minimizer can be employed.

Currently, we provide interfaces for CLASS and CAMB, as well as likelihood implementations for candl and BAOs.<fnref n="28" />

## A.2 MontePython interface

The MontePython interface allows OLÉ to (partially) replace CLASS when using MontePython, relying on the OLÉ emulator while letting MontePython handle the sampling and likelihood

<footnote n="28">²⁸Examples can be found at https://github.com/svenguenther/OLE/tree/main/OLE/examples/candl.</footnote>
calls. The interface allows OLÉ to work out of the box with existing likelihoods. It only requires setting the path to your existing MontePython installation in the MP_PATH configuration file inside your OLÉ installation, and executing the interface as you would normally execute MontePython,

python /path/to/OLE/interfaces/montepython_interface.py -p your_MCMC.param

including any additional MontePython arguments at the end. No changes to existing input parameter files are needed, but the emulator settings to pass to OLÉ can be set in the input parameter file via the new dictionary data.emulatorSettings.

The interface operates by inserting OLÉ as an intermediate layer between MontePython and the cosmology code, in this case usually CLASS. It creates copies of all functions in the Python wrapper classy and constructs emulators for those functions utilized in the likelihoods specified in the input parameter file. These emulators then approximate the output of these functions and provide the results to each likelihood code. To ensure accurate emulation, the interface classifies the classy wrapper functions based on their input and output types, allowing the emulator to correctly reproduce their behavior. The most commonly used functions suitable for OLÉ are listed in montepython_interface.py. If any new functions suitable for emulation are added to the classy wrapper, they should be added to this list, however, most CLASS modifications do not require the addition of new functions to the wrapper, in which case this step is not necessary.

Even if the interface allows OLÉ to work out of the box with existing likelihoods, it is worth considering whether a given likelihood is written in the optimal way in regards to emulation. Many existing MontePython likelihoods utilize functions that return both the cosmological results on the grid used internally by CLASS and the grid itself. For example, there are function returning both P<sub>m</sub>(k, z) and a (k, z)-grid, or the CMB C<sub>ls</sub> along with a one-dimensional l-grid. In such cases, OLÉ would, by default, emulate both the desired quantity and the coordinate values of the grid points. In addition to the computational overhead of emulating these outputs, this introduces extra uncertainty from the emulator, and a priori OLÉ does not have any protection against non-nonsensical coordinate values, e.g. multiple identical k-values or non-integer l-s.

In the case of C_l, this can be solved trivially by passing the names of outputs we are not interested in via the 'skip_emulation quantities' parameter.

In the case of large scale structure functions, the easiest solution is to replace in the likelihood codes any call to a function returning *k* and *z* grids in addition to power spectra by equivalent functions returning *P(k, z)* on an input grid, which can be defined in the likelihood's **.data** file. This requires only minor modifications to existing likelihoods, and the equivalent functions are already defined in *classy*. We provide examples of this in the form of the *euclid_photometric_alm_OLE* and *euclid_spectroscopic_OLE* likelihoods used in section 3.3, which are equivalent to the existing public *MontePython* likelihoods without the _OLE suffix as used in reference [83]. Using these likelihoods as an example, the concrete changes that were made to the likelihoods and the *classy* wrapper are as follows:
* We replaced calls to *classy* functions of the type `get_pk_and_k_and_z()` with calls of the type `get_pk()`. In practice, this means the interpolation in the $P(k, z)$-grid computed by *CLASS* happens in the wrapper, rather than in the likelihood.

* The existing public likelihoods get the reduced Hubble parameter from a call to the *classy* function `h()`. Thus, `OLÉ` may automatically try to emulate this function. Instead, we implement a check to see if this parameter is in the list of MCMC parameters, and can thus be obtained directly from the sampler. The function `h()` is used and emulated by `OLÉ` only when this is not the case. Such checks can be implemented for other background cosmological parameters needed in a given likelihood.

* The existing likelihoods call two *classy* functions to get the total matter spectrum $P_m(k, z)$ and the baryon plus CDM power spectrum $P_{cb}(k, z)$, before taking their ratio. It is more efficient for `OLÉ` to emulate directly the ratio instead of the numerator and the denominator. For this purpose, we added to *classy* a function `get_Pk_cb_m_ratio()` returning directly the ratio. For this specific case, an additional function, `get_tk()`, was also implemented, which returns a given transfer function on a provided $(k, z)$-grid, following the conventions of the `get_pk()`-type functions.

Finally, as explained in section 2.4, after the initial emulator training phase, `OLÉ` automatically disables oversampling and the fast/slow parameter decomposition by setting the `MontePython` argument `-jumping` to `global` instead of the default `fast`.

## A.3 Cobaya interface

The Cobaya interface of `OLÉ` allows the use of the `OLÉ` algorithm inside the Cobaya framework with all included likelihoods, theory codes and samplers. The implementation is such that the Cobaya code files and the overall usage of Cobaya remain unchanged. `OLÉ` only modifies parts of the code once they are loaded in memory at the time of execution.

The usage of the emulator can be turned on by giving the `theory-element` in the initial `yaml/dict` the keyword `emulate` with the value `true`. The keyword `emulatorSettings` takes a dictionary of settings to change the default `OLÉ` settings.<fnref n="29" /> When implementing an inference pipeline with the Cobaya internal `Theory` classes, it is sufficient to import them from `OLE.interfaces.cobaya_interface` while all other components can remain as plain Cobaya.<fnref n="30" /> However, there remain a few points to consider when using the Cobaya-OLÉ interface. These are:

* The interface works when Cobaya calls one theory code only (e.g., either `CAMB` or `CLASS`). The user can involve multiple theory codes in their pipeline (e.g. a Boltzmann code and a Nucleosynthesis code) only if they are nested into each other (e.g., `Primat` [100] could still be called from inside `CLASS`).

* As outlined in section 2.4 oversampling is useful in the early *burn-in* stage to converge to the relevant samples more quickly. However, once the emulator is trained, the

<footnote n="29">²⁹Examples are available at https://github.com/svenguenther/OLE/tree/main/OLE/examples/Cobaya.</footnote>

<footnote n="30">³⁰For detailed documentation on building `Theory` and `Likelihood` classes, see https://cobaya.readthedocs.io/en/latest/.</footnote>
evaluation time of the emulator becomes typically shorter than that of the likelihood. Then, it is useful to switch oversampling off. This can be done by providing the MCMC sampler with the options `'callback_function': OLE_callback_function` and `'callback_every': 1.

* Due to the particular implementation of CAMB,<fnref n="31" /> it is required to manually set a blocking of the cosmological and nuisance parameters, as in the example script shown below. Further examples are available on the OLÉ GitHub page.

* The interface has been tested for a variety of likelihoods, including CMB, BAO, lensing and SNIa likelihoods. When using likelihoods with an unusual input/output format, we strongly recommend to validate the results first before using OLÉ. Note that emulating observables with varying size is not implemented in OLÉ.

Eventually, the script to start an inference run barely changes:

```python
# import OLE-cobaya interface for your Cobaya version
from OLEYG interfaces.cobaya_interface import *

info = {
    'theory': {
        'classy': {
            'emulate': true # necessary line to request emulation
            'emulator_settings': {...}, # optional change of default settings
            ...},
            'likelihood': {
                'planck_2018_highl_plik.TTTEE': {},
                ... },
            'params': {
                ... },
            'sampler': {
                'mcmc': {
                    'callback_function': OLE_callback_function, # see 'callback_every'
                    'callback_every': 1, # to disable oversampling once trained
                    'blocking': [1, ['HO', ..., 'ns']], [5, ['A_planck', ...]], # CAMB only
                }, }, }

# run MCMC
updated_info, sampler = cobaya.run(info)
```

## B Precision settings

The accuracy of the OLÉ emulator is determined by a set of precision parameters that can either be specified by the user or take default values. For most of the examples presented

<footnote n="31">³¹The theory code of CAMB is split between `camb` and `camb_transfers`. This allows oversampling for some of the cosmological parameters. However, this behavior is not supported by OLÉ and can lead to errors.</footnote>
in section 3, we have conducted internal tests with multiple sets of OLÉ precision settings to examine the relationship between these settings, accuracy, and computational cost. For completeness, this appendix lists the precision parameters that were set to non-default values in the presented examples.

## B.1 Default precision parameters

We present an overview of the most important precision parameters, followed by their default value in bold. For a complete list, please refer to the online documentation.

* `cache_size`: 500

Maximum number of stored training data points. If more data points are to be added, the one with the smallest log-likelihood is removed.

* `min_data_points`: 80

Minimum number of points in the cache before the emulator can be trained. If this number is too small, the emulator will require many re-trainings. If it is too large, the initial data-gathering phase of OLÉ is unnecessarily long.

* `delta_loglike`: 50

This parameter sets the threshold between relevant data points (to be cached) and outliers. All points in the cache whose log-likelihood differs from the maximum log-likelihood by more than `delta_loglike` are classified as outliers and removed. If `Nsigma` and `dimensionality` (see below) are set, this parameter is ignored.

* `dimensionality`: None

As an alternative to `delta_loglike`, we can infer the threshold between relevant points and outliers by computing the $\Delta \log \mathcal{L}$ of a Gaussian distribution (with a dimensionality set by the parameter `dimensionality`) from its best fit point to `Nsigma` standard deviations. Thus, if the posterior is approximated by a Gaussian, all points in the cache lay inside the `Nsigma` contour, while points outside are classified as outliers. This parameter should be set to the total number of MCMC parameters (cosmology and nuisance parameters). If no dimensionality is specified, `delta_loglike` is used, but in general it is recommended to pass `Nsigma` and `dimensionality`.

* `Nsigma`: 3.0

This parameter sets the range around the best-fit point within which points should be cached rather than classified as outliers. To be used in conjunction with `dimensionality`.

* `min-variance-per-bin`: $5 \times 10^{-6}$

The level of compression of each observable depends on the number of PCA components. To determine this number, OLÉ increases the number of PCA components until the explained variance per bin times the bin size exceeds this parameter. Setting this variance to $\sigma^2 = 10^{-4}$ would mean that for each observable, the truncation of the expansion into PCA components introduces a typical relative error on the normalized observable of $\sigma \sim 10^{-2}$. Thus, this parameter accounts for the maximal achievable
precision of the emulator. If it is too large, results might be biases. In presence of highly correlated parameter effects, it is advisable to reduce this number by one or two orders of magnitude.

* `quality_threshold constant : 0.1`

In addition to the compression of the PCA, the second major source of error in OLÉ is the precision of the GP emulation. This is controlled by the following three precision parameters. First, `quality_threshold constant` is the constant term for determining the maximum allowed error on the log-likelihood, as in equation (2.7).

* `quality_threshold linear : 0.05`

This is the linear term for determining the maximum allowed error on the log-likelihood in equation (2.7).

* `quality_threshold quadratic : 10⁻⁴`

This is the quadratic term for determining the maximum allowed error on the log-likelihood in equation (2.7). When using `dimensionality` and `Nsigma`, this parameter is ignored and the quadratic term is computed automatically in such way that points beyond `Nsigma` are considered as outliers.

* `N质量Samples : 5`

The number of samples drawn from the emulator to estimate its own accuracy. Large values result in repeated likelihood evaluations and slow the code down. Small values increase the degree of randomness in the evaluation of the emulator accuracy, which may lead to adding unnecessary points to the cache.

## B.2 Extended cosmology example

In section 3.2, in order to get faster convergence and a smaller data set, we relax the precision parameter `quality_threshold linear` to 0.1. For parameters such as $\Omega_\nu$, whose posterior intersects the prior boundary $\Omega_\nu > 0$ and whose effect on CMB observables is far from linear near $\Omega_\nu = 0$, we want to have a very precise reconstruction of the observable. Thus, we set the sensitivity for the PCA compression to `min-variance-per-bin` to $10^{-6}$.

## B.3 Stage-IV LSS forecast example

For the example in section 3.3 we use once more `min-variance-per-bin` : $10^{-6}$ and relax the precision parameter `quality_threshold constant` to 1.0 for faster convergence.

## B.4 Early dark energy example

## B.5 Staged NUTS sampling example

For the example in section 4 we set `min-variance-per-bin` : $10^{-6}$ and `N-sigma` : 4.0.

# C Cosmological parameters

Table 4 summarizes all cosmological parameters referenced throughout this work.
<table><caption>**Table 3.** OLÉ precision settings different from default values for the NEDE example of section 3.4.</caption><thead><tr><th></th><th>High Precision</th><th>Low Precision</th></tr></thead><tbody><tr><td>quality_threshold constant</td><td>0.1</td><td>0.4</td></tr><tr><td>quality_threshold linear</td><td>0.01</td><td>0.1</td></tr><tr><td>min variance per bin</td><td>10<sup>-5</sup></td><td>10<sup>-4</sup></td></tr></tbody></table>

<table><caption>**Table 4.** Brief definitions of all cosmological parameters that appear in the manuscript.</caption><thead><tr><th>Parameter</th><th>Definition</th></tr></thead><tbody><tr><td><i>h</i></td><td>Expansion rate today <i>H</i><sub>0</sub> in units of 100 km/s/Mpc</td></tr><tr><td><i>Ω</i><sub>b</sub></td><td>Baryon fractional density</td></tr><tr><td><i>Ω</i><sub>cdm</sub></td><td>Cold dark matter fractional density</td></tr><tr><td>ln(10<sup>10</sup><i>A</i><sub>s</sub>)</td><td>Amplitude of the power spectrum of initial scalar fluctuations</td></tr><tr><td><i>n</i><sub>s</sub></td><td>Tilt of the power spectrum of initial scalar fluctuations</td></tr><tr><td><i>τ</i><sub>reio</sub></td><td>Optical depth to reionization</td></tr><tr><td><i>w</i><sub>0</sub>, <i>w</i><sub>a</sub></td><td>Dark Energy equation of state parameters in equation (3.2)</td></tr><tr><td><i>Ω</i><sub>k</sub></td><td>Mean spatial curvature parameter</td></tr><tr><td><i>σ</i><sub>8</sub></td><td>R.m.s. of matter fluctuations in a sphere of comoving radius of 8 <i>h</i><sup>-1</sup> Mpc</td></tr><tr><td><i>Ω</i><sub>ν</sub></td><td>Neutrino fractional density</td></tr><tr><td><i>Ω</i><sub>m</sub></td><td>Total matter fractional density</td></tr><tr><td>Δ<i>N</i><sub>eff</sub></td><td>Deviation of the number of effective neutrino species from the ΛCDM prediction</td></tr><tr><td>Σ<i>m</i><sub>ν</sub></td><td>Sum of neutrino masses</td></tr><tr><td><i>f</i><sub>NEDE</sub></td><td>Fraction of NEDE at the time of decay</td></tr><tr><td>log<sub>10</sub> <i>z</i><sub>*</sub></td><td>Logarithm of the redshift at NEDE decay time</td></tr><tr><td><i>w</i><sub>*</sub></td><td>Equation of state of the NEDE fluid after decay</td></tr><tr><td><i>Ω</i><sub>ϕ</sub></td><td>Density of the NEDE trigger field</td></tr></tbody></table>

# D Introduction to Gaussian processes

Gaussian Processes (GPs) are a non-parametric framework for modeling distributions over functions. They are especially well-suited for regression, classification, and optimization tasks where quantifying uncertainty is important.

## D.1 Overview

A Gaussian Process is a collection of random variables, any finite number of which have a joint Gaussian distribution. Formally, a GP is defined as:

$$f(x) \sim \mathcal{GP}(m(x), k(x, x'))$$

where:

* $m(x) = \mathbb{E}[f(x)]$ is the mean function

* $k(x, x') = \mathbb{E}[(f(x) - m(x))(f(x') - m(x'))]$ is the covariance (kernel) function

In practice, the mean function is often assumed to be zero, $m(x) = 0$, without loss of generality. The choice of kernel $k(x, x')$ determines the smoothness and generalization properties of the model.
## D.2 Gaussian process regression

Given data $\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n$ with observations $y_i = f(x_i) + \epsilon_i$, where $\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$, Gaussian Process regression involves placing a GP prior on $f(x)$ and computing the posterior distribution over functions. We write the set of arguments $x_i$ of the data set $\mathcal{D}$ as $X$. The set of points at which we want to predict the function $f$ is denoted as $X_*$.

The joint prior over the training outputs $\mathbf{f}$ and test outputs $\mathbf{f}_*$ is:

$$ \begin{bmatrix} \mathbf{f} \\ \mathbf{f}_* \end{bmatrix} \sim \mathcal{N} \left( \mathbf{0}, \begin{bmatrix} K(X, X) + \sigma_n^2 I & K(X, X_*) \\ K(X_*, X) & K(X_*, X_*) \end{bmatrix} \right) $$

The posterior predictive distribution for test points $X_*$ is:

$$ \mathbf{f}_* | X, y, X_* \sim \mathcal{N}(\mu_*, \Sigma_*) $$

with:

$$ \mu_* = K(X_*, X)[K(X, X) + \sigma_n^2 I]^{-1}y $$

$$ \Sigma_* = K(X_*, X_*) - K(X_*, X)[K(X, X) + \sigma_n^2 I]^{-1}K(X, X_*) $$

## D.3 Kernel functions

The kernel function $k(x, x')$ encodes prior assumptions about the emulated function. One common choice is the RBF kernel:

$$ k(x, x') = \sigma_f^2 \exp \left( - \frac{(x - x')^2}{2\ell^2} \right) $$

which describes a function that is smooth with a characteristic length scale $\ell$. Kernels often include hyperparameters such as $\ell$, which are typically optimized by maximizing the marginal likelihood.

## D.4 Advantages and limitations

Gaussian Processes provide uncertainty estimates for predictions while being highly flexible. However, their computational cost scales as $\mathcal{O}(n^3)$ for $n$ training points, making them inefficient for large datasets. Their performance depends heavily on kernel choice and hyperparameter tuning, so having some prior knowledge of the behavior of the modeled data is advantageous.

Various sparse GP methods, such as inducing point approximations, have been proposed to reduce computational complexity for large datasets.

## References

[1] LSST SCIENCE and LSST PROJECT collaborations, LSST science book, version 2.0, arXiv:0912.0201 [INSPIRE].

[2] ACT collaboration, The Atacama Cosmology Telescope: a measurement of the cosmic microwave background power spectra at 98 and 150 GHz, JCAP 12 (2020) 045 [arXiv:2007.07289] [INSPIRE].
[3] DESI collaboration, DESI 2024 VI: cosmological constraints from the measurements of baryon acoustic oscillations, JCAP 02 (2025) 021 [arXiv:2404.03002] [INSPIRE].

[4] PLANCK collaboration, Planck 2018 results. V. CMB power spectra and likelihoods, Astron. Astrophys. 641 (2020) A5 [arXiv:1907.12875] [INSPIRE].

[5] SPT-3G collaboration, Measurement of the CMB temperature power spectrum and constraints on cosmology from the SPT-3G 2018 TT, TE, and EE dataset, Phys. Rev. D 108 (2023) 023510 [arXiv:2212.05642] [INSPIRE].

[6] EUCLID THEORY WORKING GROUP collaboration, Cosmology and fundamental physics with the Euclid satellite, Living Rev. Rel. 16 (2013) 6 [arXiv:1206.1225] [INSPIRE].

[7] SKA collaboration, Cosmology with phase 1 of the Square Kilometre Array: red book 2018. Technical specifications and performance forecasts, Publ. Astron. Soc. Austral. 37 (2020) e007 [arXiv:1811.02743] [INSPIRE].

[8] A. Spurio Mancini et al., CosmoPower: emulating cosmological power spectra for accelerated Bayesian inference from next-generation surveys, Mon. Not. Roy. Astron. Soc. 511 (2022) 1771 [arXiv:2106.03846] [INSPIRE].

[9] EUCLID collaboration, Euclid preparation: II. The EuclidEmulator — a tool to compute the cosmology dependence of the nonlinear matter power spectrum, Mon. Not. Roy. Astron. Soc. 484 (2019) 5509 [arXiv:1809.04695] [INSPIRE].

[10] G. Aricò et al., The BACCO simulation project: a baryonification emulator with neural networks, Mon. Not. Roy. Astron. Soc. 506 (2021) 4070 [arXiv:2011.15018] [INSPIRE].

[11] G. Aricò, R.E. Angulo and M. Zennaro, Accelerating large-scale-structure data analyses by emulating Boltzmann solvers and Lagrangian perturbation theory, arXiv:2104.14568 [DOI:10.12688/openreseurope.14310.2] [INSPIRE].

[12] K.R. Moran et al., The Mira–Titan universe — IV. High-precision power spectrum emulation, Mon. Not. Roy. Astron. Soc. 520 (2023) 3443 [arXiv:2207.12345] [INSPIRE].

[13] M. Mohri, A. Rostamizadeh and A. Talwalkar, Foundations of machine learning, second edition, The MIT Press, U.S.A. (2018).

[14] S. Günther, Uncertainty-aware and data-efficient cosmological emulation using Gaussian processes and PCA, arXiv:2307.01138 [INSPIRE].

[15] A.R. Khalife et al., Review of Hubble tension solutions with new SH0ES and SPT-3G data, JCAP 04 (2024) 059 [arXiv:2312.09814] [INSPIRE].

[16] D. Blas, J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: approximation schemes, JCAP 07 (2011) 034 [arXiv:1104.2933] [INSPIRE].

[17] J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) IV: efficient implementation of non-cold relics, JCAP 09 (2011) 032 [arXiv:1104.2935] [INSPIRE].

[18] A. Lewis, A. Challinor and A. Lasenby, Efficient computation of CMB anisotropies in closed FRW models, Astrophys. J. 538 (2000) 473 [astro-ph/9911177] [INSPIRE].

[19] C. Howlett, A. Lewis, A. Hall and A. Challinor, CMB power spectrum parameter degeneracies in the era of precision cosmology, JCAP 04 (2012) 027 [arXiv:1201.3654] [INSPIRE].

[20] K. Heitmann et al., The Coyote universe I: precision determination of the nonlinear matter power spectrum, Astrophys. J. 715 (2010) 104 [arXiv:0812.1052] [INSPIRE].

[21] K. Heitmann et al., The Coyote universe II: cosmological models and precision emulation of the nonlinear matter power spectrum, Astrophys. J. 705 (2009) 156 [arXiv:0902.0429] [INSPIRE].
[22] E. Lawrence et al., *The Coyote universe III: simulation suite and precision emulator for the nonlinear matter power spectrum*, Astrophys. J. 713 (2010) 1322 [arXiv:0912.4490] [INSPIRE].

[23] S. Bhattacharya, S. Habib, K. Heitmann and A. Vikhlinin, *Dark matter halo profiles of massive clusters: theory vs. observations*, Astrophys. J. 766 (2013) 32 [arXiv:1112.5479] [INSPIRE].

[24] S. Agarwal et al., *PkANN — I. Non-linear matter power spectrum interpolation through artificial neural networks*, Mon. Not. Roy. Astron. Soc. 424 (2012) 1409 [arXiv:1203.1695] [INSPIRE].

[25] J. Kwan, S. Bhattacharya, K. Heitmann and S. Habib, *Cosmic emulation: the concentration-mass relation for wCDM universes*, Astrophys. J. 768 (2013) 123 [arXiv:1210.1576] [INSPIRE].

[26] J. Kwan et al., *Cosmic emulation: fast predictions for the galaxy power spectrum*, Astrophys. J. 810 (2015) 35 [arXiv:1311.6444] [INSPIRE].

[27] S. Agarwal et al., *pkann — II. A non-linear matter power spectrum interpolator developed using artificial neural networks*, Mon. Not. Roy. Astron. Soc. 439 (2014) 2102 [arXiv:1312.2101] [INSPIRE].

[28] K. Heitmann et al., *The Coyote universe extended: precision emulation of the matter power spectrum*, Astrophys. J. 780 (2014) 111 [arXiv:1304.7849] [INSPIRE].

[29] K. Heitmann et al., *The Mira–Titan universe: precision predictions for dark energy surveys*, Astrophys. J. 820 (2016) 108 [arXiv:1508.02654] [INSPIRE].

[30] E. Lawrence et al., *The Mira–Titan universe II: matter power spectrum emulation*, Astrophys. J. 847 (2017) 50 [arXiv:1705.03388] [INSPIRE].

[31] J. DeRose et al., *The Aemulus project I: numerical simulations for precision cosmology*, Astrophys. J. 875 (2019) 69 [arXiv:1804.05865] [INSPIRE].

[32] T. McClintock et al., *The Aemulus project II: emulating the halo mass function*, Astrophys. J. 872 (2019) 53 [arXiv:1804.05866] [INSPIRE].

[33] Z. Zhai et al., *The Aemulus project III: emulation of the galaxy correlation function*, Astrophys. J. 874 (2019) 95 [arXiv:1804.05867] [INSPIRE].

[34] T. McClintock et al., *The Aemulus project IV: emulating halo bias*, arXiv:1907.13167 [INSPIRE].

[35] EUCLID collaboration, *Euclid preparation: IX. EuclidEmulator2 — power spectrum emulation with massive neutrinos and self-consistent dark energy perturbations*, Mon. Not. Roy. Astron. Soc. 505 (2021) 2840 [arXiv:2010.11288] [INSPIRE].

[36] S. Bocquet et al., *The Mira–Titan universe. III. Emulation of the halo mass function*, Astrophys. J. 901 (2020) 5 [arXiv:2003.12116] [INSPIRE].

[37] M.-F. Ho, S. Bird and C.R. Shelton, *Multifidelity emulation for the matter power spectrum using Gaussian processes*, Mon. Not. Roy. Astron. Soc. 509 (2021) 2551 [arXiv:2105.01081] [INSPIRE].

[38] D. Jamieson et al., *Field-level neural network emulator for cosmological N-body simulations*, Astrophys. J. 952 (2023) 145 [arXiv:2206.04594] [INSPIRE].

[39] J. Kwan et al., *Galaxy clustering in the Mira–Titan universe. I. Emulators for the redshift space galaxy correlation function and galaxy–galaxy lensing*, Astrophys. J. 952 (2023) 80 [arXiv:2302.12379] [INSPIRE].
[40] A. Upadhye et al., Cosmic-Ev: an-emulator for the non-linear neutrino power spectrum, Mon. Not. Roy. Astron. Soc. 530 (2024) 743 [arXiv:2311.11240] [INSPIRE].

[41] S.T. Brown et al., ARTEMIS emulator: exploring the effect of cosmology and galaxy formation physics on Milky Way-mass haloes and their satellites, Mon. Not. Roy. Astron. Soc. 532 (2024) 1223 [arXiv:2403.11692] [INSPIRE].

[42] Y. Yang, S. Bird and M.-F. Ho, Ten-parameter simulation suite for cosmological emulation beyond ΛCDM, Phys. Rev. D 111 (2025) 083529 [arXiv:2501.06296] [INSPIRE].

[43] J.R. Gair et al., The hitchhiker's guide to the galaxy catalog approach for dark siren gravitational-wave cosmology, Astron. J. 166 (2023) 22 [arXiv:2212.08694] [INSPIRE].

[44] S. Bird et al., An emulator for the Lyman-α forest, JCAP 02 (2019) 050 [arXiv:1812.04654] [INSPIRE].

[45] C. Pedersen et al., An emulator for the Lyman-α forest in beyond-ΛCDM cosmologies, JCAP 05 (2021) 033 [arXiv:2011.15127] [INSPIRE].

[46] R. Kugel et al., FLAMINGO: calibrating large cosmological hydrodynamical simulations with machine learning, Mon. Not. Roy. Astron. Soc. 526 (2023) 6103 [arXiv:2306.05492] [INSPIRE].

[47] M. Schaller et al., The FLAMINGO project: baryon effects on the matter power spectrum, Mon. Not. Roy. Astron. Soc. 539 (2025) 1337 [arXiv:2410.17109] [INSPIRE].

[48] M. Walther et al., Emulating the Lyman-α forest 1D power spectrum from cosmological simulations: new models and constraints from the eBOSS measurement, JCAP 05 (2025) 099 [arXiv:2412.05372] [INSPIRE].

[49] M. Kaplinghat, L. Knox and C. Skordis, Rapid calculation of theoretical cmb angular power spectra, Astrophys. J. 578 (2002) 665 [astro-ph/0203413] [INSPIRE].

[50] R. Jimenez, L. Verde, H. Peiris and A. Kosowsky, Fast cosmological parameter estimation from microwave background temperature and polarization power spectra, Phys. Rev. D 70 (2004) 023005 [astro-ph/0404237] [INSPIRE].

[51] W.A. Fendt and B.D. Wandelt, Pico: parameters for the impatient cosmologist, Astrophys. J. 654 (2006) 2 [astro-ph/0606709] [INSPIRE].

[52] J. Albers et al., CosmicNet. Part I. Physics-driven implementation of neural networks within Einstein-Boltzmann solvers, JCAP 09 (2019) 028 [arXiv:1907.05764] [INSPIRE].

[53] S. Günther et al., CosmicNet II: emulating extended cosmologies with efficient and accurate neural networks, JCAP 11 (2022) 035 [arXiv:2207.05707] [INSPIRE].

[54] A. Nygaard, E.B. Holm, S. Hannestad and T. Tram, CONNECT: a neural network based framework for emulating cosmological observables and cosmological parameter inference, JCAP 05 (2023) 025 [arXiv:2205.15726] [INSPIRE].

[55] M. Bonici, F. Bianchini and J. Ruiz-Zapatero, Capse.jl: efficient and auto-differentiable CMB power spectra emulation, Open J. Astrophys. 7 (2024) [arXiv:2307.14339] [INSPIRE].

[56] D. Piras and A. Spurio Mancini, CosmoPower-JAX: high-dimensional Bayesian inference with differentiable cosmological emulators, Open J. Astrophys. 6 (2023) [arXiv:2305.06347] [INSPIRE].

[57] H.T. Jense et al., A complete framework for cosmological emulation and inference with CosmoPower, RAS Tech. Instrum. 4 (2025) rzaf002 [arXiv:2405.07903] [INSPIRE].

[58] M. Bonici, G. D'Amico, J. Bel and C. Carbone, Effort: a fast and differentiable emulator for the effective field theory of the large scale structure of the universe, arXiv:2501.04639 [INSPIRE].
[59] M.D. Schneider, O. Holm and L. Knox, Intelligent design: on the emulation of cosmological simulations, Astrophys. J. 728 (2011) 137 [arXiv:1002.1752] [INSPIRE].

[60] C.M. Bishop, Pattern recognition and machine learning, first edition, Springer (2007).

[61] T. Brinckmann and J. Lesgourgues, MontePython 3: boosted MCMC sampler and other features, Phys. Dark Univ. 24 (2019) 100260 [arXiv:1804.07261] [INSPIRE].

[62] J. Torrado and A. Lewis, Cobaya: code for Bayesian analysis of hierarchical physical models, JCAP 05 (2021) 057 [arXiv:2005.05290] [INSPIRE].

[63] J. Bradbury et al., JAX: composable transformations of Python+NumPy programs, http://github.com/google/jax (2018).

[64] J.-E. Campagne et al., JAX-COSMO: an end-to-end differentiable and GPU accelerated cosmology library, Open J. Astrophys. 6 (2023) 1 [arXiv:2302.05163] [INSPIRE].

[65] L. Balkenhol, C. Trendafilova, K. Benabed and S. Galli, candl: cosmic microwave background analysis with a differentiable likelihood, Astron. Astrophys. 686 (2024) A10 [arXiv:2401.13433] [INSPIRE].

[66] D. Piras et al., The future of cosmological likelihood-based inference: accelerated high-dimensional parameter estimation and model comparison, Open J. Astrophys. 7 (2024) [arXiv:2405.12965] [INSPIRE].

[67] T. Pinder and D. Dodd, GPJax: a Gaussian process framework in JAX, J. Open Source Softw. 7 (2022) 4455.

[68] A. Gelman and D.B. Rubin, Inference from iterative simulation using multiple sequences, Statist. Sci. 7 (1992) 457 [INSPIRE].

[69] S. Brooks, A. Gelman, G. Jones and X.L. Meng, Handbook of Markov chain Monte Carlo, CRC Press, U.S.A. (2011).

[70] A. Lewis, GetDist: a python package for analysing Monte Carlo samples, JCAP 08 (2025) 025 [arXiv:1910.13970] [INSPIRE].

[71] R. Kumar, C. Carroll, A. Hartikainen and O. Martin, ArviZ a unified library for exploratory analysis of Bayesian models in Python, J. Open Source Softw. 4 (2019) 1143.

[72] PLANCK collaboration, Planck 2018 results. VIII. Gravitational lensing, Astron. Astrophys. 641 (2020) A8 [arXiv:1807.06210] [INSPIRE].

[73] PAN-STARRS1 collaboration, The complete light-curve sample of spectroscopically confirmed SNe Ia from Pan-STARRS1 and cosmological constraints from the combined Pantheon sample, Astrophys. J. 859 (2018) 101 [arXiv:1710.00845] [INSPIRE].

[74] A.J. Ross et al., The clustering of the SDSS DR7 main galaxy sample — I. A 4 per cent distance measure at z = 0.15, Mon. Not. Roy. Astron. Soc. 449 (2015) 835 [arXiv:1409.3242] [INSPIRE].

[75] BOSS collaboration, The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, Mon. Not. Roy. Astron. Soc. 470 (2017) 2617 [arXiv:1607.03155] [INSPIRE].

[76] E.V. Linder, Exploring the expansion history of the universe, Phys. Rev. Lett. 90 (2003) 091301 [astro-ph/0208512] [INSPIRE].

[77] M. Chevallier and D. Polarski, Accelerating universes with scaling dark matter, Int. J. Mod. Phys. D 10 (2001) 213 [gr-qc/0009008] [INSPIRE].
[78] W. Fang, W. Hu and A. Lewis, *Crossing the phantom divide with parameterized post-Friedmann dark energy*, Phys. Rev. D **78** (2008) 087303 [arXiv:0808.3125] [INSPIRE].

[79] J.C. Hill et al., *Atacama Cosmology Telescope: constraints on prerecombination early dark energy*, Phys. Rev. D **105** (2022) 123536 [arXiv:2109.04451] [INSPIRE].

[80] B. Audren, J. Lesgourgues, K. Benabed and S. Prunet, *Conservative constraints on early cosmology: an illustration of the Monte Python cosmological parameter inference code*, JCAP **02** (2013) 001 [arXiv:1210.7183] [INSPIRE].

[81] T. Sprenger et al., *Cosmology in the era of Euclid and the Square Kilometre Array*, JCAP **02** (2019) 047 [arXiv:1801.08331] [INSPIRE].

[82] EUCLID collaboration, *Euclid: validation of the MontePython forecasting tools*, Astron. Astrophys. **682** (2024) A90 [arXiv:2303.09451] [INSPIRE].

[83] EUCLID collaboration, *Euclid preparation — LIV. Sensitivity to neutrino parameters*, Astron. Astrophys. **693** (2025) A58 [arXiv:2405.06047] [INSPIRE].

[84] A. Chatrchyan, F. Niedermann, V. Poulin and M.S. Sloth, *Confronting cold new early dark energy and its equation of state with updated CMB, supernovae, and BAO data*, Phys. Rev. D **111** (2025) 043536 [arXiv:2408.14537] [INSPIRE].

[85] M. Garny, F. Niedermann, H. Rubira and M.S. Sloth, *Hot new early dark energy bridging cosmic gaps: supercooled phase transition reconciles stepped dark radiation solutions to the Hubble tension with BBN*, Phys. Rev. D **110** (2024) 023531 [arXiv:2404.07256] [INSPIRE].

[86] D. Brout et al., *The Pantheon+ analysis: cosmological constraints*, Astrophys. J. **938** (2022) **110** [arXiv:2202.04077] [INSPIRE].

[87] S.G. Nash, *Newton-type minimization via the Lanczos method*, SIAM J. Numer. Anal. **21** (1984) 770.

[88] S. Duane, A.D. Kennedy, B.J. Pendleton and D. Roweth, *Hybrid Monte Carlo*, Phys. Lett. B **195** (1987) 216 [INSPIRE].

[89] J. Nocedal and S.J. Wright, *Numerical optimization*, Springer, New York, NY, U.S.A. (2006).

[90] M.D. Hoffman and A. Gelman, *The no-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo*, arXiv:1111.4246 [INSPIRE].

[91] D.P. Kingma and J. Ba, *Adam: a method for stochastic optimization*, arXiv:1412.6980 [INSPIRE].

[92] DeepMind, I. Babuschkin et al., *The DeepMind JAX ecosystem*, http://github.com/deepmind (2020).

[93] J. Robnik, G.B. De Luca, E. Silverstein and U. Seljak, *Microcanonical Hamiltonian Monte Carlo*, arXiv:2212.08549 [INSPIRE].

[94] P. Virtanen et al., *SciPy 1.0 — fundamental algorithms for scientific computing in python*, Nature Meth. **17** (2020) 261 [arXiv:1907.10121] [INSPIRE].

[95] L. Balkenhol, *Compressed 'CMB-lite' likelihoods using automatic differentiation*, Open J. Astrophys. **8** (2025) [arXiv:2412.00826] [INSPIRE].

[96] ACT collaboration, *The Atacama Cosmology Telescope: DR4 maps and cosmological parameters*, JCAP **12** (2020) 047 [arXiv:2007.07288] [INSPIRE].
[97] A. Mootoovaloo, J. Ruiz-Zapatero, C. García-García and D. Alonso, *Assessment of gradient-based samplers in standard cosmological likelihoods*, Mon. Not. Roy. Astron. Soc. 534 (2024) 1668 [arXiv:2406.04725] [INSPIRE].

[98] D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, *emcee: the MCMC hammer*, Publ. Astron. Soc. Pac. 125 (2013) 306 [arXiv:1202.3665] [INSPIRE].

[99] P. Virtanen et al., *SciPy 1.0 — fundamental algorithms for scientific computing in python*, Nature Meth. 17 (2020) 261 [arXiv:1907.10121] [INSPIRE].

[100] C. Pitrou, A. Coc, J.-P. Uzan and E. Vangioni, *Precision Big Bang nucleosynthesis with the new code PRIMAT*, JPS Conf. Proc. 31 (2020) 011034 [arXiv:1909.12046] [INSPIRE].

[101] E. Jones et al., *SciPy: open source scientific tools for python*, https://scipy.org/ (2001), accessed 22 October 2014.

[102] J.D. Hunter, *Matplotlib: a 2D graphics environment*, Comput. Sci. Eng. 9 (2007) 90 [INSPIRE].

[103] S. van der Walt, S.C. Colbert and G. Varoquaux, *The NumPy array: a structure for efficient numerical computation*, Comput. Sci. Eng. 13 (2011) 22 [arXiv:1102.1523] [INSPIRE].