Cosmo-Learn: code for learning cosmology using different methods and mock data

TL;DR

Cosmo-Learn is an open-source Python package that simulates cosmological data and applies various statistical and machine learning methods for inference, aiding research and education in cosmology.

Contribution

It introduces a flexible, modular framework for benchmarking and comparing different cosmological inference techniques using realistic simulated data.

Findings

Demonstrates internal consistency of simulated data with input cosmology

Supports multiple inference methods including MCMC, neural networks, and Gaussian processes

Provides a standardized platform for method comparison and educational purposes

Abstract

We present cosmo_learn, an open-source python-based software package designed to simulate cosmological data and perform data-driven inference using a range of modern statistical and machine learning techniques. Motivated by the growing complexity of cosmological models and the emergence of observational tensions, cosmo_learn provides a standardized and flexible framework for benchmarking cosmological inference methods. The package supports realistic noise modeling for key observables in the late Universe, including cosmic chronometers, supernovae Ia, baryon acoustic oscillations, redshift space distortions, and gravitational wave bright sirens. We demonstrate the internal consistency of the simulated data with the input cosmology via residuals and parameter recovery using a fiducial $w$CDM model. Built-in learning and inference modules include traditional Markov Chain Monte Carlo, as well as more recent approaches such as genetic algorithms, Gaussian processes, Bayesian ridge regression, and artificial neural networks. These methods are implemented in a modular and extensible architecture designed to facilitate comparisons across inference strategies in a common pipeline. By providing a flexible and transparent simulation and learning environment, cosmo_learn supports both educational and research efforts at the intersection of cosmology, statistics, and machine learning.

Tables

Table 1. Table 1: Key components of a GA.

Fitness Function	Quantifies how well an individual performs relative to the target objective.
Selection	Determines which individuals are chosen to reproduce. Common techniques include roulette-wheel or tournament selection, which favor individuals with higher fitness.
Crossover	Combines parent chromosomes to produce offspring, enabling exploration of new regions in the parameter space by mixing genetic material.
Mutation	Introduces random genetic changes, to avoid local optima and maintain genetic diversity.
Elitism	Preserves top-performing individuals across generations.

Table 2. Table 2: Comparison on SNe Ia mock data (reference universe 14000605; processor: Apple M4); t t is training time in seconds, Mean, RMS are the residuals mean and root-mean-square values over the training and the test sets, χ red 2 \chi^{2}_{\mathrm{red}} is the reduced chi-squared, χ red 2 = ∑ ( ( y rec − y data ) / σ res ) 2 / N \chi^{2}_{\rm red}=\sum\left((y_{\rm rec}-y_{\rm data})/\sigma_{\rm res}\right)^{2}/N , where y rec y_{\rm rec} , y data y_{\rm data} are the means of the reconstruction and the data at each redshift, and σ res 2 \sigma_{\rm res}^{2} is the combined variance of the reconstruction and the data.

		Training set			Test set
Method	$t$ (s)	Mean	RMS	${\chi }_{\mathrm{red}}^2$	Mean	RMS	${\chi }_{\mathrm{red}}^2$
GP	30.593	$-$0.0016	0.2387	0.9614	$-$0.0185	0.2223	1.0030
BRR	0.001	$-$0.0017	0.2382	0.9608	$-$0.0171	0.2212	0.9978
ANN	67.362	$-$0.0018	0.2384	0.4939	$-$0.0179	0.2216	0.4841
MCMC	30.605	$-$0.0046	0.2381	0.9600	$-$0.0189	0.2200	0.9876
GA-Fisher	9.836	$-$0.0076	0.2383	0.9619	$-$0.0226	0.2210	1.0037

Equations

$$\mu (z_i,\Theta )=m-M=5{\log }_{10}[d_L(z_i,\Theta )]+25,$$

$$d_L(z_i,\Theta )=c(1+z_i){\int }_0^{z_i}\frac{d{z}^{\prime }}{H(z^{\prime },\Theta )},$$

$${\delta }^{\prime \prime }(a)+\left(\frac{3}{a}+\frac{H^{\prime }(a)}{H(a)}\right){\delta }^{\prime }(a)-\frac{3}{2}\frac{{\Omega }_{m0}{G}_{\mathrm{eff}}(a,k)/G_{\mathrm{N}}}{a^{5}H(a{)}^2/H_0^2}\delta (a)=0,$$

$$\delta (a)=a{}_2{F}_1\left(\frac{1}{3},1,\frac{11}{6};a^3\left(1-{\Omega }_{m0}^{-1}\right)\right),$$

$$f{\sigma }_8(z)=-{\sigma }_8(1+z)\frac{{\delta }^{\prime }(z)}{\delta (z=0)}.$$

$$D_M(z)=c{\int }_0^z\frac{d{z}^{\prime }}{H(z^{\prime })},$$

$$D_V(z)=(1+z{)}^{-1}{\left[D_M(z{)}^2\frac{cz}{H(z)}\right]}^{\frac{1}{3}}.$$

$$r_{\mathrm{D}}={\int }_{z_{\mathrm{D}}}^{\infty }\frac{c_s(z)}{H(z)}dz,$$

$$h(\tau )\propto \frac{{\mathcal{M}}^{5/3}}{D_L(z,\Theta )}f(\tau {)}^{2/3},$$

$$\frac{{\sigma }_{\text{lens}}(z,\Theta )}{d_L(z,\Theta )}\sim 0.066{\left(\frac{1-(1+z{)}^{-0.25}}{0.25}\right)}^{1.8}\left[1-\frac{0.3}{\pi /2}\arctan \left(\frac{z}{z^{\ast }}\right)\right],$$

$$\frac{{\sigma }_v(z,\Theta )}{d_L(z,\Theta )}\sim \left[1+\frac{c(1+z{)}^2}{H(z,\Theta )d_L(z,\Theta )}\right]\frac{\sqrt{⟨v^2⟩}}{c}.$$

$$\frac{{\sigma }_{\text{LISA}}(z,\Theta )}{d_L(z,\Theta )}\sim 0.05\left(\frac{d_L(z,\Theta )}{36.6\text{Gpc}}\right).$$

$${\sigma }_{\text{photo}}(z,\Theta )\sim 0.03(1+z).$$

$$\sigma (z,\Theta )=\sqrt{{\sigma }_{\text{lens}}^2(z,\Theta )+{\sigma }_v^2(z,\Theta )+{\sigma }_{\text{LISA}}^2(z,\Theta )+{\sigma }_{\text{photo}}^2(z)}.$$

$$p_{\text{DE}}=w{\rho }_{\text{DE}},$$

$$H^2(z)=H_0^2\left[{\Omega }_m(1+z{)}^3+{\Omega }_r(1+z{)}^4+{\Omega }_{\text{DE}}(1+z{)}^{3(1+w)}\right],$$

$$p(\Theta |X_{\mathrm{Obs}})=\frac{L(\Theta |X_{\mathrm{Obs}})p(\Theta )}{\int L({\Theta }^{\prime }|X_{\mathrm{Obs}})p({\Theta }^{\prime })d{\Theta }^{\prime }},$$

$$H(z)\sim \mathcal{G}\mathcal{P}\left(\mu (z),k(z,\tilde{z})\right),$$

$$k(z,\tilde{z})={\sigma }_H^2\left[\frac{l_H}{(z-\tilde{z}{)}^2+l_H^2}\right],$$

$$H^{\ast }\sim \mathcal{N}\left({\mu }^{\ast },K(Z^{\ast },Z^{\ast })\right),$$

$$\overline{H}(z_i)\sim \mathcal{N}\left(\mu (z_i),k(Z,Z)+C\right),$$

$$\left(\begin{array}{c}\overline{H}\\ H^{\ast }\end{array}\right)\sim \mathcal{N}\left(\left(\begin{array}{c}\mu \\ {\mu }^{\ast }\end{array}\right),\left(\begin{array}{cc}k(Z,Z)+C& k(Z,Z^{\ast })\\ k(Z^{\ast },Z)& k(Z^{\ast },Z^{\ast })\end{array}\right)\right).$$

$$p\left(\overline{H}|Z,{\sigma }_H,l_H\right)=\int p\left(\overline{H}|H,Z\right)p\left(H|Z,{\sigma }_H,l_H\right)dH,$$

$$\ln \mathcal{L}$$

$$=-\frac{1}{2}{\left(\overline{H}-\mu \right)}^T{\left[k(Z,Z)+C\right]}^{-1}\left(\overline{H}-\mu \right)-\frac{1}{2}\ln |k(Z,Z)+C|-\frac{n}{2}\ln 2\pi ,$$

$$R=\alpha ||w|{|}^2$$

$$p(w|\lambda )=\mathcal{N}(0,{\lambda }^{-1}I).$$

$$p(y|X,w,\alpha )=\mathcal{N}(Xw,{\alpha }^{-1}I).$$

$$f(x)=\{\begin{array}{ll}x& \text{if }x>0\\ \alpha (e^x-1)& \text{if }x\le 0\end{array},$$

$$L1=\sum _i|{{\rm Y}}_{\mathrm{train}}(z_i)-{{\rm Y}}_{\mathrm{pred}}(z_i)|,$$

Full text

1 Introduction

The search for new physics has taken on central importance in meeting the growing challenges of modern cosmology. Before the precision cosmological era, supernovae were found fainter than expected in the cosmological model without dark energy, and the fit is improved by introducing $\Lambda$, which supported the observed accelerated expansion. Now, precision measurements of the expansion of the Universe [1, 2] show a faster rate than predicted by the $\Lambda$CDM concordance model [3], when early Universe data sets are considered. As a consequence, there has been an increased effort to produce realistic potential resolutions to this open question. Moreover, the standard concordance model is undergoing a reassessment of the constituent components in its makeup, including the direct observation of cold dark matter (CDM) [4, 5, 6], the theoretically problematic cosmological constant [7, 8], and observational measurements of cosmic inflation [9, 10]. The latest addition to the family of open questions is that of the Hubble tension [11, 12], which expresses the discordance between estimates of the Hubble constant in the local Universe and concordance model based measurements using early-time observations. This effect is compounded by the possibility of tensions in the rate of growth of large scale structures [11, 13, 14, 15, 16, 17]. There are several competing interpretations of the question of cosmic tensions [11, 18], but it is progressively becoming a characteristic between direct measurements in the late Universe [19, 20, 21, 22], and inferred cosmological parameter constraints using measurements from the early Universe [23, 24, 25, 26, 27].

As the wealth of data releases continues to grow, the likelihood of the Hubble tension question being the result of a single systematic that permeates through the data pipelines in question becomes increasingly less likely. In this context, there has been a spectrum of potential directions beyond $\Lambda$CDM in which to build potential cosmological solutions. These have taken the form of a number of promising expressions including modifications in the physics beyond recombination with early dark energy models [28] being the most promising in this regime, additional extra relativistic degrees of freedom [29], as well as the plethora of suggestions for modified gravitational physics [30, 31, 32, 33, 34, 35] in different sectors of phenomenology. This spectrum of different physics at all scales of phenomenology has produced a deluge of potential directions in which to build potentially viable physical models, which has only increased the need for a more standardized benchmark on which to compare cosmologies.

The vast array of cosmological frameworks on which models are built are in direct competition with each other, and in many cases produce models that are degenerate in terms of their constraints through observational measurements. This leads to an environment where having precisely comparable tests has become an issue of central importance.

An additional direction is cosmological reconstruction, where the focus shifts from constraining cosmological models to directly inferring the underlying expansion history and growth rate of structure from data. Such approaches aim to remain agnostic to the detailed form of the cosmological model, allowing for a more data-driven view of the Universe. Reconstructions of the Hubble parameter, the dark energy equation of state, or the growth rate of structure provide a flexible way of assessing whether observational data favor the $\Lambda$CDM concordance model [36, 37, 38, 39]. Cosmological reconstructions serve as a bridge between theoretical model-building and empirical investigations, and they highlight where tensions may arise without relying on model priors. In this sense, reconstructions complement both traditional parameter estimation and the newer statistical and machine learning tools [40].

It is in this context that we aim to homogenize some aspects of potential benchmark testing for model comparisons. We do this through a new open and transparent pipeline for mock data generation, where mock data is generated as a representation of different real observational measurements for different classes of cosmological models. In this way, mock data around an underlying fiducial cosmology can be built by mimicking the real data noise distribution and assuming a mean from a noisy realization of the underlying cosmology. This will be important for producing standardized processes for various cosmological scenarios. As a preliminary approach, we include data sets in the local Universe in this first iteration of release. This includes a cosmic chronometer data set release [41, 42, 43, 44, 45, 46], the Pantheon+ compilation [19, 47, 48], redshift space distortions measurements [49] related to baryonic acoustic oscillations [50], and gravitational wave bright sirens. This gives a combination of relatively standard and new data sets under which to consider the resulting mock data generation pipeline. Moreover, this gives a healthy product that probes deep into relatively high redshift space. Through this implementation, we hope to provide a community toolkit that can more easily tie into processing pipelines requiring a mock data generation implementation.

The establishment of an accessible and transparent mock data generator is important to build toolkits for the next generation of cosmological probes given their nuances and individual complexities. As these surveys build momentum toward their eventual data releases, collaborations such as Euclid [51] and the Large Synoptic Survey Telescope [52], among others, stand to benefit from standardized mock data pipelines that allow methods to be developed, tested, and compared in a controlled setting. Moreover, having an easy-to-use mock data generator specific to cosmology will also accelerate efforts to develop toolkits related to other cosmological probes which may involve more nuanced information such as in the cases of weak lensing, galaxy cluster counts or halo occupation distributions, among others.

Another aspect of the issue of constraining cosmological models in a statistically robust manner is related to the method by which the constraint analysis takes place. Traditionally, Bayesian inference is performed using Markov chain Monte Carlo (MCMC) methods [53, 54, 55], in which a Markov chain is constructed to approximate the posterior distribution of the model parameters. By generating a large number of correlated samples from this chain, one can approximate credible intervals, and other summary statistics of the posterior. While effective, this method is computationally expensive both when large data sets are considered, and similarly when complex physical models with multiple nuisance parameters and numerically solved parameters are incorporated. In recent years, there has been a concerted effort to meet these challenges through the development of novel statistical tools based on machine learning in order to understand cosmological models, as well as to constrain them. One popular approach has been to use simulation-based inference (SBI) techniques where the outputs of a full numerical simulator are emulated using different methods. These techniques offer an important way to accelerate the running of complex cosmological simulations while retaining numerically precise model predictions. SBIs primarily depend on different neural network configurations which aim to mimic the (1) Bayes posterior direction [56, 57], (2) the individual likelihood output for each iteration [58, 59], or (3) the marginalized outputs for model parameters [60, 61]. Other approaches exist to accelerate the constraint of complex cosmological models. However, SBIs form a key part of this spectrum of new methods to assess cosmological models. Another important technique is that of genetic algorithms (GAs) [62, 63, 64, 17] where a potentially unbounded space of model realizations can be searched using techniques based on natural selection-inspired approaches. Here, the parameter values define the infinite space of potential realizations, and a potentially faster approach to the best-fit constraint is approximated using the GA method. This may be more effective for the increasingly large space of highly complex cosmological models.

Another key technique in the toolkit of novel methods is that of Gaussian processes (GPs) [65, 66] where an underlying data set is used to train a kernel function, which is a representation of the relationships between the constituents of the data set. The goal of this approach is to produce a kernel that best represents the data set combination. On the other hand, this produces a powerful way to generate data sets that change the data set distribution while also being genuine representations of the measurements. In the same spirit, Bayesian ridge regression (BRR) [67, 68] offers a parametric approach to estimating data points and their associated uncertainties for a redshift range. Here, the parametric approach is dynamic in that it can be self-modulated in relation to the complexity of the underlying data sets under consideration. These two methods are closely related to one another.

These methods are powerful in their predictive power but they take on a number of statistical assumptions in their estimation method. In order to counter this potential issue, neural network architectures [69, 70, 71, 72, 73] offer a stronger statistically robust base on which to produce representative data sets that mimic the underlying real measurements. Neural network methods are based on a biologically-inspired approach in which neurons are organized into layers with connections between neurons across different layers. These connections and the neuron responses are controlled by non-physical hyperparameters. Through this approach, a realistic estimation of the data point values and associated uncertainties can be directly obtained in a natural way.

In this work, we explore the question of how these methods compare against each other. We do this using different combinations of mock generated data in order to give a clear benchmark on which to probe these techniques. In this way, we aim to extend the discussion in the literature on possible novel methods on which to rely beyond traditional approaches. In Section 2, the mock data generation approach is outlined, as are the observed data sets and the cosmological model on which this is based. The validation approach is then outlined in Section 3 where $w$CDM is used to estimate cosmological parameters, which acts as the underlying fiducial model. Then, in Section 4 the methods outlined here are applied to the mock data sets in order to reconstruct the cosmic evolution deep into intermediate redshift space. This reconstruction of the expansion and large scale structure formation may provide key insights in disentangling the degeneracies between cosmological models and understanding nuances in the comparison of their underlying physics.

The cosmo_learn package can be found at https://github.com/reggiebernardo/cosmo_learn; which come with tutorial notebooks for an interactive and easier familiarization of the modules. This uses the modules astropy [74, 75, 76], scikit-learn [77], reFANN [78] (ANN based on pytorch [79]), gwcatalog [80, 81, 82], emcee [83], geneticalgorithm, as well as numpy [84], scipy [85], matplotlib [86], and corner [87].

2 Simulating cosmological data

In this section, we present the default setup used by cosmo_learn to generate a mock data realization based on a provided set of authentic samples and a ‘true’ cosmology. Then, we describe in brief each of the mock data sets that can be simulated with cosmo_learn: SNe Type Ia distance moduli, CC, BAO, RSD growth rate measurements, and GW standard sirens’ luminosity distances.

2.1 Mock data generation

In cosmo_learn, mock data $\left(z_{\rm mock},y_{\rm mock}\pm\Delta y_{\rm mock}\right)$ is by default generated based on actual observed samples $\left(z_{\rm obs},y_{\rm obs}\pm\Delta y_{\rm obs}\right)$ and a ‘true’ cosmology $y_{\rm true}[H_{0},{\mathbf{\Omega}}]$ in the following steps:

Fix the redshifts to match actual ones, $z_{\rm mock}=z_{\rm obs}$; 2. 2.

Draw the mean from normal distribution around the true cosmology with a variance fixed by observation, $y_{1}={\cal N}\left(y_{\rm true},\left(\Delta y_{\rm obs}\right)^{2}\right)$; 3. 3.

Final mock data $\left(z_{\rm mock},y_{\rm mock}\pm\Delta y_{\rm mock}\right)$ is $\left(z_{\rm obs},y_{1}\pm\Delta y_{\rm{obs}}\right)$.

This default setting leads to a conservative mock data realization to avoid possible spurious forecasts. For the purposes of demonstration and entry for students, we found that this is sufficient. Experts on the other hand can modify the function how_to_mock in cosmo_learn.py to play with the mock data generation, such as consider a full covariance $\Sigma$ and draw the redshift samples $z_{\rm mock}$ using an empirical Gamma or Beta distributions [70, 80].

The present version ‘true’ cosmology is based on $w$CDM (Section 2.3). Thus, an instance of cosmo_learn can be initiated by providing a value of the Hubble constant $H_{0}$, nonrelativistic matter density parameter $\Omega_{m0}$, the smoothed amplitude of matter fluctuations $\sigma_{8}$, and the dark energy equation of state $w$. It is useful to provide a seed specific to an instance or realization of mock data (in a simulated universe) for reproducibility. A minimal example to start an instance of cosmo_learn with mock data generated by the default setting is

my_cl=CosmoLearn([H0, Om0, w0, s8], seed=14000605) my_cl.make_mock(mock_keys=[‘CosmicChronometers’, ‘SuperNovae’, ‘BrightSirens’])

where the $H_{0}$ is in units of km s Mpc*-1* and $\Omega_{m0},w_{0},S_{8}$ are dimensionless; e.g., [50, 23]

H0, Om0, w0, s8 = 67.74, 0.3095, -0.997, 0.834

The above two lines generate an instance of cosmo_learn in universe 14000605 with a true cosmology given by the DESI Year 1 Early Data Release results [50], with cosmic chronometers, supernovae, and GW bright sirens mock data. Obviously, the choice of true cosmology based on a set of cosmological parameters $(H_{0},\Omega_{m0},w,S_{8})$ is arbitrary for learning purposes and completely up to the user; where $S_{8}=\sigma_{8}\sqrt{\Omega_{m0}/0.3}$. Keys

‘BaryonAcousticOscillations’, ‘RedshiftSpaceDistorsions’

for the method make_mock are similarly available to generate baryon acoustic observation samples and redshift space distortion growth rate measurements.

Mock data and their corresponding residuals across 100 realizations are shown in Figure 1. This shows that the generated data across multiple realizations are centered on the true cosmology, i.e., residuals are white, as expected since they were drawn from a Gaussian distribution. Mock data on a single realization, e.g., reference universe 14000605, can manifest $\sim$1-2$\sigma$ deviations due to the noise in the actual data. Later tests (Section 3) will show that cosmological parameters can be inferred with $\sim$1-2$\sigma$ deviations relative to the true cosmology, consistent with expectations.

Mock data generated in cosmo_learn are divided into a training set (90% of the full data) and a test set (10%) that are specific to the random state assigned per realization. This is used for training the methods integrated into the package (Section 4).

2.2 Standard cosmological data

We describe the data sets considered in cosmo_learn for generating mock data. The package is built to be versatile and allow additional and other data sets to be added as needed. Below, we discuss the built-in data sets that the package uses, which are centered in the local Universe, but which can be extended.

Cosmic Chronometers

Cosmic Chronometers (CC) provide a practical way to directly constrain the Hubble expansion rate $H(z,\Theta)$ at different redshifts due to the differential aging method on which these data sets are built [88]. In our implementation, we adopt 31 data points as compiled from Ref. [43, 89, 45, 90, 42, 41, 44]. The spectroscopic dating technique used in this method expresses the differential aging of cosmological local galaxies with their redshifts, $\Delta z/\Delta t$, so that the Hubble expansion at those redshifts can be determined $H_{\mathrm{obs}}(z)=-(1+z)^{-1}\Delta z/\Delta t$. In this way, CC data is a reliable source of data points that are free of cosmological models, and moreover their direct nature renders them independent of any complex calculations. Their only dependence is on the absolute aging of galaxies which is a fairly well understood feature in astrophysics [91].

Supernovae Type Ia

The Pantheon+ ($\mathrm{PN}^{+}$ & SH0ES) sample [19, 47, 48] contains one of the largest samples of type Ia supernovae (SNe) with 1701 measurements spanning a redshift range of $0.0$ to $2.26$. The Pantheon+ analysis uses Cepheid hosts as distance anchors together with geometric calibrators. In this approach, the degeneracy between the absolute magnitude $M$ and $H_{0}$ is effectively addressed through the use of SH0ES prior information. The expansion rate at different redshifts, $H(z,\Theta)$, is measured through the observed apparent magnitude, $m$, through the relation

[TABLE]

where $d_{L}(z_{i},\Theta)$ represents the luminosity distance to each redshift expressed as

[TABLE]

where $c$ is the speed of light. We consider the subset $z>0.01$ to be conservative of potential issues with very low redshift samples [92, 93].

RSD Growth Data

Redshift-space distortions (RSD) offer a powerful observational handle on the large-scale structure of the Universe by measuring anisotropies in the galaxy distribution; caused by peculiar velocities superimposed on the Hubble flow. These anisotropies carry crucial information about the growth of cosmic structures and allow us to track the evolution of the linear matter density contrast, $\delta(z)$, across redshift. Importantly, $\delta(z)$ is not only a tracer of structure formation but also a sensitive probe of gravity on cosmological scales. At the linear perturbation level, the evolution of the density contrast is governed by the differential equation

[TABLE]

where $G_{\rm eff}(a,k)$ characterizes potential modifications to the gravitational strength on large scales [94, 95, 96, 97, 98, 49]. In GR, $G_{\rm eff}(a,k)/G_{\rm N}=1$, while departures from unity can signal deviations from Einsteinian gravity. Hence, by constraining the evolution of $\delta(z)$ through RSD data, one effectively places bounds on gravity itself at linear cosmological scales.

Observationally, this information is captured through measurements of $f\sigma_{8}(z)$, the product of the linear growth rate $f(z)=d\ln\delta(z)/d\ln a(z)$ and the root-mean-square matter fluctuation amplitude $\sigma_{8}$ on scales of $8\,h^{-1}\,\text{Mpc}$. This quantity effectively encapsulates both the rate of structure growth and the normalization of the matter power spectrum.

In the standard cosmological model, the linear matter overdensity evolves as [98]:

[TABLE]

where $\,{}_{2}F_{1}(a,b,c;x)$ is the hypergeometric function. This expression admits generalization to models with a constant dark energy equation of state, a.k.a. $w$CDM [98].

The growth rate observable is then given by [16]:

[TABLE]

Here, $\sigma_{8}(z)=\sigma_{8}\delta(z)/\delta(z=0)$ reflects the redshift evolution of the fluctuation amplitude.

For our analysis, we employ the compilation of RSD measurements from [49].

Baryon Acoustic Oscillations

Baryonic acoustic oscillations (BAO) map out the fluctuations of the density profile of visible matter. These repeated oscillations can be used as a standard ruler to precisely measure cosmological distances. In this work, the Dark Energy Spectroscopic Instrument (DESI) data release 1 (DR1) [50] was used, which consists of BAO distance scales binned in several redshifts using different tracers, namely the bright galaxy sample (BGS) at $0.1<z<0.4$, the luminous red galaxy (LRG) sample at $0.4<z<0.6$ and $0.6<z<0.8$, a combination of LRG and emission line galaxies (ELG) at $0.8<z<1.1$, the ELG sample at $1.1<z<1.6$, the quasar sample at $0.8<z<2.1$, and the Lyman-alpha forest sample at $1.77<z<4.16$. These seven redshift bins make up the DESI DR1 sample.

For the BAO DESI DR1 sample, we compute the Hubble distance $D_{H}(z)$ through the expression $D_{H}(z)=c/H(z)$, which gives the effective distance along the line-of-sight. The effective distance transverse to the line-of-sight is given by

[TABLE]

which is also related to the angular diameter distance by $D_{A}=(1+z)^{-1}D_{M}(z)$. The volume average distance is then given by

[TABLE]

BAO measurements are inherently degenerate with the comoving sound horizon at the drag epoch, $r_{\rm D}$, which acts as the standard ruler. This scale corresponds to the distance a sound wave can travel in the primordial photon-baryon fluid before recombination, and is given by

[TABLE]

where $z_{\rm D}\simeq 1059.94$ is the redshift of the drag epoch [99], and $c_{s}(z)$ is the sound speed of the baryon-photon fluid. BAO observables are typically reported as ratios such as $D_{H}(z)/r_{\rm D}$, $D_{M}(z)/r_{\rm D}$, or $D_{V}(z)/r_{\rm D}$, ensuring that model predictions can be directly compared with observational data.

GW Bright Sirens

GWs are tidal ripples in spacetime caused by the dynamics of massive astrophysical objects and can manifest as strains that are captured by large-scale interferometers such as the ground-based LIGO-Virgo-Kagra (LVK) detectors and the anticipated space-borne Laser Interferometer Space Antenna (LISA). GWs may originate from the coalescence of compact binaries, rapid rotation of astrophysical bodies, and even phase transitions from the early universe [100, 101, 102, 103, 104, 105, 106, 107, 108]. The signals from these events carry direct information on the luminosity distance $D_{L}(z,\Theta)$ to the source, making them powerful standard sirens [109]. For coalescing compact binaries, the luminosity distance can be inferred from the strain measurements captured by interferometric GW detectors. If the merger event is accompanied by an electromagnetic (EM) counterpart (i.e. quasar-like luminous flares, radio jets, flares), the redshift $z$ of the source can be independently determined, either spectroscopically or photometrically [105]. These joint detections, where both GW and EM signals are observed, are referred to as bright standard sirens, and they offer a direct probe of the cosmological expansion through the distance-redshift relation (2).

At leading Newtonian order, the GW strain amplitude $h(\tau)$ is related to the luminosity distance by

[TABLE]

where $\mathcal{M}$ is the chirp mass of the binary system, $f(\tau)$ is the instantaneous GW frequency, $h(\tau)$ is the strain, and $\tau$ is the merger time of the binary [102, 109]. Note that merger time scales are typically much shorter than cosmological time scales, even for massive black hole binaries. Since all these quantities are measurable, the luminosity distance can be inferred directly. When the redshift is also available through EM observations, one can construct the Hubble diagram using bright standard sirens in an analogous manner to SNe Ia, following (2).

To construct a mock catalog of standard siren events, we adopt the redshift distribution shown in Figure 2, corresponding to one of the most promising eLISA mission configurations studied in [110]. The figure reflects the expected number of detectable standard sirens with EM counterparts, based on simulations incorporating the eLISA sensitivity curve [110, 106]. The redshift distribution varies according to the underlying binary population model.

To simulate the events, we first interpret the redshift histogram as a probability distribution function that is nonzero only in the interval $z\in(0.1,9)$. For each event, we draw a trial redshift from a uniform distribution, and then draw a trial event count from another uniform distribution between the minimum and maximum values of a predetermined expected number of events or mission lifetime; as a default, we consider a 5-year lifetime. The redshift is accepted if the sampled count is less than the maximum expectation; otherwise, the process is repeated. This rejection sampling is continued until the desired number of events for a given binary model is reached. Each accepted redshift is then associated with a source luminosity distance via (2).

To make the mock catalog more realistic, we account for anticipated observational uncertainties to each event [80, 81, 82]. We incorporate weak lensing uncertainties using the fit from [110, 111],

[TABLE]

with $z^{*}=0.073$ as used in [112] to control the growth of lensing-induced scatter. We also include uncertainties due to peculiar velocities $v$ following [110, 113]:

[TABLE]

Then, the instrumental noise from LISA is approximated to scale linearly with distance as described in [114]:

[TABLE]

Finally, we account for redshift measurement uncertainties from the EM counterpart. For redshifts below $z\sim 2$, which are typically measured spectroscopically, the uncertainties are negligible. For $z\gtrsim 2$, where photometric methods dominate, we use the estimate from [115]:

[TABLE]

These sources of error are added in quadrature to obtain the total uncertainty:

[TABLE]

We sample the luminosity distance from a Gaussian distribution $\mathcal{N}(d_{L}(z,\Theta),\sigma(z,\Theta))$ centered at the theoretical $d_{L}(z,\Theta)$ with standard deviation $\sigma(z,\Theta)$. This procedure yields a full mock standard siren catalog, which can be combined with simulated data from other cosmological probes. Note that the errors explicitly depend on the cosmological parameters. A sample realization is shown in Figure 1.

This module in cosmo_learn is based on our modified version of the gwcatalog codebase [80, 81, 82]. We consider a default Pop III binary population model; users are free to change this via the attribute pop=‘Pop III’ to the cosmo_learn instance; e.g., key in my_cl.pop_model=‘Delay’ or my_cl.pop_model=‘No Delay’ for the the delay and no delay models, respectively. The mission lifetime can be changed using option years=5 when calling make_mock.

2.3 $w$CDM cosmology

The $w$CDM model is the simplest phenomenological extension of the standard $\Lambda$CDM model, obtained by generalizing the equation of state of dark energy from a fixed value ($w=-1$) to a free parameter $w$,

[TABLE]

where $w\in\mathbb{R}$ is constant in time. The value $w=-1$ recovers the $\Lambda$CDM scenario, while $w<-1$ and $w>-1$ correspond to phantom and quintessence-like behaviors [116], respectively. The Hubble expansion rate in this model is given by

[TABLE]

where $\Omega_{\text{DE}}=1-\Omega_{m}-\Omega_{r}$ under the assumption of spatial flatness, and the usual scaling laws are assumed for matter ($\propto(1+z)^{3}$) and radiation ($\propto(1+z)^{4}$).

The $w$CDM model is frequently used in cosmological analyses as a minimal benchmark to assess deviations from $\Lambda$CDM, particularly in testing the robustness of parameter estimation methods and exploring the sensitivity of current and future surveys to the properties of dark energy. Because it involves only one additional parameter at the background, it provides a controlled framework for quantifying the impact of relaxing the cosmological constant assumption without invoking a full scalar field model or modified gravity. Furthermore, for linearized perturbations, $w$CDM admits an analytical solution for the density contrast of nonrelativistic matter [96, 98], allowing for numerically efficient computation of observables depending on linear perturbations such as the growth rate $f\sigma_{8}(z)$ [49]. Few other cosmological models have similar analytical expressions available, making the $w$CDM model a suitable basis for studying cosmology through numerical analysis.

The $w$CDM model has been extensively employed in data analyses from major cosmological surveys, such as Planck [99], the Dark Energy Survey (DES) [24], and the Pantheon compilation of type Ia supernovae [47]. These studies have placed tight constraints on $w$, with current data being consistent with $w=-1$ at the few-percent level, though future surveys are expected to improve these bounds significantly.

cosmo_learn uses the analytical expressions for the background and linear perturbation in $w$CDM in [98]. By default, DE perturbations are turned off through the attribute de_model=‘no pert’. This can be changed to one with DE perturbations but with the sound speed, $c_{s}$, fixed to zero as de_model=‘static’. Alternatively, DE perturbations can be turned on with a nonzero sound speed by fixing de_model=‘dynamic’ and configuring the attribute k2cs2=1e-10, corresponding to the product $k^{2}c_{s}^{2}$ for modes with a wavenumber $k$. For further details, see [98].

3 Testing mock data

To validate the mock data generation beyond the appearance of white residuals, we demonstrate that the input cosmology can be statistically recovered through parameter inference. Specifically, we use the built-in MCMC engine in cosmo_learn which leverages emcee to estimate the posterior distributions of cosmological parameters under a $w$CDM model. The resulting constraints, aggregated over 100 mock realizations, are shown in Figure 3, where all available observables (Section 2.2) are utilized.

The constraints from individual realizations such as the reference universe demonstrate that the input cosmological parameters are recoverable within $\sim$1–2$\sigma$, depending on the noise level. Furthermore, when posteriors from all realizations are stacked, the resulting ensemble distribution closely approximates a Gaussian centered on the true cosmology. This supports the conclusion that the residuals are effectively white when analyzed over many mock data sets (see Figure 1).

A particularly noteworthy result is the precise constraint on the sound horizon at the drag epoch, $r_{\rm D}$, which shows sub-$1\sigma$ deviations on average. This parameter remains fixed for a given instance of cosmo_learn and is unaffected by noise, even though it must still be sampled during BAO-based inference.

These MCMC-based results confirm that the mock data generated by cosmo_learn are statistically consistent with the input cosmology. In the next section, we describe additional inference tools beyond MCMC that are integrated into the package for flexibility and benchmarking.

4 Learning Methods

We describe learning methods in cosmo_learn (MCMC, GA, GP, BRR, ANN) that users can readily play with to seamlessly perform parameter estimation and cosmological reconstruction.

4.1 MCMC

Markov chain Monte Carlo methods (MCMC) [53, 54, 55] has become a core tool in estimating posterior parameter distributions for cosmological models in Bayesian inference studies. In this framework, a Monte Carlo approach is taken wherein a large number of random samples are drawn from a normal distribution. The parameters that make up each draw are connected to physical predictions through a system of representative equations such as the Friedmann or Boltzmann equations. The MCMC pipeline also necessitates a set of observational data to be used in order to compare the system predictions in the Monte Carlo segment. As the number of random draws for a sample increases, the parameter means contained in the sample gradually asymptote to the underlying means. This is the benefit of the approach since calculating the mean directly from a system of equations can be cumbersome for most scenarios, while the Monte Carlo method can be used more directly in a vast array of scenarios. The addition of the Markov chain property to the analysis pipeline embeds the idea that the random draws for each sample iteration are connected together in a sequence such that successive draws will either be improved samples in the sequence or discarded in the chain. Notice the important property that this process only requires information on the previous sample iteration and not any of the elements before that.

This approach is particularly useful for Bayesian inference [117, 118, 119] analyses where posterior distributions are often times difficult to determine. Using Bayes’ theorem, for a data set $X_{\rm Obs}$ the inferred posterior distribution will be defined by

[TABLE]

where $L(\Theta|X_{\rm Obs})$ represents the likelihood of predicting the data given a physical model with parameter $\Theta$, called the likelihood function, $p(\Theta)$ is the prior distribution which normally represents a Crenel function on the parameters that make up the underlying physical model. The denominator is called the marginal likelihood and denotes the likelihood produced by the data itself. This probability of evidence likelihood integral is at the core of the MCMC estimator approach. Naturally, one would like to find the posterior distribution of the parameters given the data. To achieve this, the MCMC pipeline takes successive iterations of $\Theta^{\prime}$ to estimate the averaged sum approximation of this integral. The end result is that the marginal likelihood is approximated by the mean for $n$ samples drawn giving $(1/n)\Sigma_{i}L(\Theta^{\prime}|X_{\rm Obs})$ for each drawn $\Theta^{\prime}$. Through these successive iterations, the final goal is to determine the posterior distributions of the parameters in $\Theta$ given the observational data $X_{\rm Obs}$.

The results of MCMC over a 100 realizations are shown in Figure 3. In cosmo_learn, MCMC in a realization can be performed as follows (after initializing a cosmo_learn instance my_cl and mock data):

initialize wcdm likelihood

llprob=lambda x: my_cl.llprob_wcdm(x, prior_dict=prior_dict, rd_fid_prior=rd_fid_prior)

setup and perform mcmc

p0=[70, 0.3, -1, 0.8, 147] # initial point dres=[0.05, 0.005, 0.01, 0.01, 0.005] # resolution nwalkers=15; nburn=100; nmcmc=2000 # number of walkers, burn in, steps per walker my_cl.get_mcmc_samples(nwalkers, dres, llprob, p0=p0, nburn, nmcmc)

This results in MCMC samples stored in my_cl.mcmc_samples that can be plotted using standard packages such as corner or getdist. cosmo_learn has a built-in method show_param_posterior that uses the former. In the reference realization (seed 14000605) the MCMC samples are shown in Figure 4 together with the corresponding results of GA-Fisher (Section 4.2).

In addition to corner plots, best fits with error bands can be easily visualized using a method show_bestfit_curve, as shown in Figure 5, for MCMC and GA-Fisher inferred parameters. These are constructed by resampling the posterior as a Gaussian distribution with a mean and covariance provided by the sampled parameters, and for each sampled point drawing the corresponding observable, making a distribution of functions. Note that this method and the default corner plotter can be used for both MCMC and GA-Fisher, which directly uses the samples stored in the cosmo_learn instance.

It is also worth noting that the likelihood initialization comes with a default set of priors for the parameters, that users can freely play with by modifying the corresponding python dictionaries, prior_dict and rd_fid, as follows:

priors (flat for H0, Om0, w0, s8, Gaussian for rd_fid)

prior_dict={‘H0_min’: 0, ‘H0_max’: 100, ‘Om0_min’: 0, ‘Om0_max’: 1,
‘w0_min’: -10, ‘w0_max’: 10, ‘s8_min’: 0.2, ‘s8_max’: 1.5} rd_fid_prior={‘mu’: 147.46, ‘sigma’: 0.28} llprob=lambda x: my_cl.llprob_wcdm(x, prior_dict=prior_dict, rd_fid_prior=rd_fid_prior)

The likelihood can be extended beyond $w$CDM, should users wish to. A tailored likelihood, e.g., for CPL or axion-like early dark energy, can be setup and given as input llprob to the method get_mcmc_samples.

4.2 GA-Fisher: Genetic Algorithm + Fisher Forecast

Genetic algorithm (GA) is a nature-inspired, population-based optimization strategy that emulates the biological principles of natural evolution [62, 63, 64, 17]. It is classified as a metaheuristic method, meaning it does not rely on gradients or derivatives to search for optima. Instead, it heuristically explores the solution space, with the goal of identifying the fittest solution within a population of candidate solutions. GAs are especially powerful in high-dimensional or multimodal problems, where local optima may obscure the global minimum. Under suitable conditions, GAs are known to converge to the global optimum [120].

GA has been applied in various areas of science, including high-energy physics [121], gravitational wave astronomy [122], and increasingly in cosmology. They were introduced to cosmological analysis in [123] as a model-independent reconstruction tool to alleviate the arbitrariness of dark energy parametrizations. Follow-up studies such as [124, 125] established GA as a flexible tool for cosmological parameter estimation and uncertainty quantification. A comprehensive and modern overview of GA for cosmological inference is provided in [62]. We refer interested readers to [64, 63] for a pedagogical introduction to GA and nature-inspired methods for cosmological parameter estimation.

In the context of cosmology, an individual in the GA population represents a realization of cosmological parameters, e.g., $\left(H_{0},\Omega_{m0},w_{0},w_{a}\right)$. Each individual is encoded as a chromosome, consisting of a set of genes, each corresponding to a parameter (see Figure 2 of [64]). The GA evolves this population over generations via selection, crossover, and mutation, with the aim of improving the fitness; typically defined via a cost or likelihood function related to observational data.

The key ingredients and steps of GA are shown in Table 1 and Figure 6.

At its core, GA ranks individuals based on fitness, promotes the fittest to reproduce through crossover and mutation, and ensures survival of the best candidates via elitism. The iterative cycle continues until convergence is achieved or a maximum number of generations is reached.

In cosmological applications, each gene typically represents a parameter such as $H_{0}$, $\Omega_{m0}$, $w_{0}$, or $w_{a}$, while the chromosome is a vector of these parameters. For example, in Chevallier–Polarski–Linder (CPL) dark energy parametrization, the chromosome could be $(H_{0},\Omega_{m0},w_{0},w_{a})$ [62].

Beyond traditional parameter estimation, GA has been deployed in more sophisticated frameworks. In the pioneering works [123, 124, 125], grammatical evolution (GE)—an extension of GA—was used to reconstruct cosmological functions without a fixed parametric form. GA has also been used to eliminate the arbitrariness in kernel choice for Gaussian Process regression [126], optimize neural network architectures [127], and improve spectroscopic modeling and model selection [128, 129, 130, 131].

In cosmo_learn, GA can be performed using a method get_gaFisher_samples, e.g., in $w$CDM:

llprob=lambda x: my_cl.llprob_wcdm(x) # log likelihood function fitness_func=lambda x: -2*llprob(x) # fitness function prior_ga=[[0, 100], [0, 1], [-10, 10], [0, 1.5], [145, 149]] my_cl.get_gaFisher_samples(fitness_func, prior_ga, llprob)

This scans the same parameter space as the built-in MCMC, $(H_{0},\Omega_{m0},w,S_{8},r_{\rm D})$, and so the priors [min, max] for each parameter must be given in the same order. Note that our hybrid method applies a Fisher matrix uncertainty estimate on top of optimization, thus requiring an input likelihood probability. Mainly, GA is used to search the global optimum of the likelihood in parameter space, then the likelihood is used to forecast the uncertainty given the Fisher matrix approach. The resulting method is referred to as GA-Fisher.

Figures 4-5 illustrate a representative result of cosmological parameter estimation using the GA-Fisher method in the reference realization 14000605. In this case, GA correctly identifies the best fit parameters corresponding to the maximum likelihood point in the parameter space which agree with the MCMC-derived best fit within $2\sigma$. This agreement extends to the confidence regions: the contours produced by GA-Fisher closely match those from MCMC, both in shape and size. This similarity arises because GA-Fisher uses the likelihood function to estimate uncertainties via the Fisher matrix. Such consistency between the two approaches is expected primarily when the data set is sufficiently constraining, i.e., when multiple cosmological probes such as CC, SNe, RSD, BAO, and GWs are combined. However, it is important to emphasize that the level of agreement between GA-Fisher and MCMC is realization-dependent, driven by the underlying noise in the mock data. While realizations such as the reference case yield good agreement, others may exhibit more significant discrepancies.

The package comes with the following default settings for GA: population size 100, number of generations 1000, parent selection rate 0.3, mutation probability 0.3, uniform crossover rate 0.5. These can be changed to user-preferred values as follows,

my_cl.ga_params[‘max_num_iteration’]=100 my_cl.ga_params[‘population_size’]=50 my_cl.ga_params[‘mutation_probability’]=0.3 my_cl.ga_params[‘elit_ratio’]=0.01 my_cl.ga_params[‘crossover_probability’]=0.5 my_cl.ga_params[‘parents_portion’=0.3 my_cl.ga_params[‘crossover_type’]=’uniform’

Lastly the convergence curve (fitness function as a function of generations) can be shown by adding the option convergence_curve=False when calling the method get_gaFisher_samples.

cosmo_learn uses the geneticalgorithm python library to perform GA.

4.3 Gaussian Processes

A Gaussian process (GP) is a generalization of a single Gaussian distribution to a series of Gaussian distributions [65, 66]. This is normally based on a set of Gaussian points in a training data set. Here, some reconstructed function $H\equiv H(z)$ is defined as a Gaussian process with mean $\mu(z)$ and variance $\sigma^{2}(z)$ for any point $z$ within the function domain. The reconstructed profile is formed by considering the relationship that neighbouring points $\tilde{z}$ have on the functional value at $z$. This is quantified through the covariance function ${\rm cov}\left(H(z),H(\tilde{z})\right)=k\left(z,\tilde{z}\right)$, which expresses the strength and separation magnitude that neighbouring points will have on each other. Then, the function can be written as a Gaussian process

[TABLE]

formed with only the relationship between points in the data set and does not include any physical parameters of some underlying model.

The covariance function can be defined through a kernel which depends on the relationship between neighbouring points of the functional $H$. As an example, we consider the Cauchy kernel; but there are many other possible ways to relate neighbouring points along the function. This is defined through

[TABLE]

where $\sigma_{H}$ and $l_{H}$ are constants that quantify the amplitude and distance at which pairs of points are correlated with each other. For the Cauchy kernel, the point square distance in the denominator means that the further points are from each other the weaker their covariance will be. This decay is polynomial which can be more realistic than other options such as exponential scenarios. However, in all cases, it is the non-parametric constants that characterize the resulting functional behavior for $H(z)$. It is important to emphasize that while the $\sigma_{H}$ and $l_{H}$ parameters do relate an expected statistical behavior of the points along the $H(z)$, there is no information embodied about any underlying physical model.

Taking some set of points $Z^{*}=\{z^{*}\}$, one can define a mean $\mu^{*}$ and kernel $\left[K(Z,Z)\right]_{ij}=k(z_{i}^{*},z_{j}^{*})$ from which the function values $H^{*}=H(z_{i}^{*})$ are described by

[TABLE]

where $\mathcal{N}$ refers to the GP as evaluated at the specific points in $Z^{*}$. While the values of $H^{*}$ form a normal distribution, its values at specific points $z_{i}^{*}$ are determined by the correlations embodied in the kernel function and respective mean $\mu^{*}$. The function $H^{*}$ is very arbitrary and is only set by setting values for the mean and the kernel constants, or statistical hyperparameters. Another way to make this determination is to use observational data to set these parameter values and then reconstruct the function at other points where there is no data.

One can consider observations information $Z=\{z_{i}\}$ which have known functional values $(z_{i},\overline{H}(z_{i}))$, where $\overline{H}(z_{i})$ must have Gaussian uncertainties. These points form part of a normal distribution or GP evaluated at specific points at which data is known

[TABLE]

where $C$ is the covariance matrix which accounts for correlations between the data. For scenarios where there are no correlations, the covariance matrix takes on the form $C={\rm Diag}\left(\sigma_{i}^{2}\right)$. Together with GP, this data can be used to reconstruct information at the unknown points $Z^{*}$ using observations at $Z$ points. This is achieved by first combining both distributions

[TABLE]

It is through this combination of observational data and target reconstruction points, that the reconstructed function’s profile can be obtained. This is achieved through a supervised training process in which the marginal likelihood

[TABLE]

is maximized, where the $H$ represents GP function values, and $Z$ are points where observational points exist. Notice that $Z^{*}$ does not feature in this integral, meaning that the points at which the reconstruction takes place are not part of the marginalization process. Once the GP hyperparameters are constrained, the locations at which reconstruction is to take place can be performed.

Given a Gaussian prior $H|Z,\sigma_{H},l_{H}\sim\mathcal{N}\left(\mu,k(Z,Z)\right)$ and taking the observational data to be a Gaussian distribution at every point $\overline{H}|H\sim\mathcal{N}\left(H,C\right)$, results in the log marginal likelihood

[TABLE]

which accounts for covariances between the observational data set, and provides a likelihood through which to optimize the nonphysical hyperparameters. For most scenarios, the mean can be set to zero without loss of generality provided this is done before the hyperparameters are fit. The GP hyperparameters will be valued accordingly to account for all information obtained from the training data [37, 132, 133]. In this way, cosmological data can be used to realize reconstructed diagrams for various cosmological parameters, including the Hubble diagram [132, 134, 135, 37, 133, 14, 136, 137, 138, 139, 140, 16, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 40].

In cosmo_learn, a GP can be trained with a single command following the initialization of a mock data set:

my_cl.train_gp()

This function fits a separate GP model for each cosmological probe present in mock_data. The training uses a default of n_restarts_optimizer = 10 to help avoid local minima in the hyperparameter space.

By default, the package employs a radial basis function (RBF or Gaussian) kernel, but alternative kernels can be specified via the kernel_key option in train_gp. Available choices include Matern, RationalQuadratic (which generalizes the Cauchy kernel), ExpSineSquared, and DotProduct. In all cases, the hyperparameter initial values and bounds follow scikit-learn’s defaults.

GP regression is non-parametric and highly flexible, allowing for smooth interpolation of the data with uncertainty estimates. For CC, BAO, RSD, and GW data sets, the model is trained in linear redshift space, while for SNe data, it operates in $\log_{10}(z)$ to better accommodate the steep slope at low redshift. A comparison of default GP and BRR reconstructions using the mock reference cosmology is shown in Figure 7.

cosmo_learn performs GP regression using scikit-learn’s GaussianProcessRegressor module [77].

4.4 Bayesian Ridge Regression

Bayesian regression techniques offer a probabilistic approach to linear modeling, where regularization parameters are inferred from the data rather than fixed a priori. This allows the model to adapt its complexity to the data at hand. Unlike classical ridge regression, where the regularization strength $\alpha$ is manually set, Bayesian regression treats it as a hyperparameter drawn from a prior distribution [67, 68].

Ridge regression, also known as L2 regularization, is a linear regression technique used to mitigate multicollinearity in multiple regression data. Multicollinearity occurs when two or more independent variables are highly correlated, leading to unstable coefficient estimates. In ridge regression, a regularization term $R$ is added to the ordinary least squares (OLS) cost function:

[TABLE]

where $\alpha$ is the regularization hyperparameter and $w$ is the coefficient vector. The regularization term penalizes large coefficients, helping to reduce model variance and prevent overfitting.

Bayesian Ridge Regression builds on this by placing a Gaussian prior over the coefficients $w$, with zero mean and precision $\lambda$:

[TABLE]

The likelihood is modeled as a Gaussian distribution of the targets $y$ around the linear model $Xw$ with noise precision $\alpha$:

[TABLE]

To make the model fully Bayesian, Gamma priors are placed over both $\alpha$ and $\lambda$, the conjugate priors for Gaussian precision. The hyperparameters are estimated by maximizing the log marginal likelihood (also known as evidence approximation). This leads to an automatic trade-off between model complexity and data fidelity.

As with GP regression, Bayesian Ridge Regression (BRR) training in cosmo_learn can be performed with a single command:

my_cl.train_brr()

This routine trains an independent BRR model for each cosmological probe available in mock_data, using polynomial basis functions of degree n_order, which defaults to 3. The polynomial order determines the model’s capacity to capture non-linear trends in the data: higher orders yield more flexible fits but may risk overfitting.

BRR extends linear regression by introducing priors over the model coefficients, allowing for automatic regularization. The model internally applies inverse-variance weighting based on the provided uncertainties and initializes the hyperparameters for the coefficient and noise precisions with alpha_init = 1.0 and lambda_init = $10^{-3}$, respectively. These priors are iteratively updated during training to maximize the marginal likelihood.

For the SNe data set, training is carried out in $\log_{10}(z)$ instead of $z$ to better resolve the structure of the Hubble diagram at low redshift. Reconstructions using the default BRR settings are shown alongside Gaussian Process regression in Figure 7 for direct comparison.

cosmo_learn uses scikit-learn’s BayesianRidge estimator [77] to perform BRR, following the formulation of Refs. [153, 154].

4.5 Neural Nets

Artificial neural networks (ANN) can be used as another method by which to reconstruct cosmological evolution parameters in an accurate and fast manner [155]. ANN structures are composed of input and output layers through which interface with the input data and output predictions can be respectively made. A series of parameters connecting hidden layers can then be optimized to mimic real training data relationships. Each layer is composed of a number of neurons that connect to neurons in the layers preceding or proceeding that layer.

As an illustrative example, we show Figure 8 where a series of three hidden layers structure the NN neurons and which has three inputs together with two outputs. In this way, an input signal would traverse the entirety of the network architecture in order to produce an output. The NN is constructed such that a linear transformation (composed of linear weights and biases) is employed for each of the different layers.

The neurons in the NN are constructed out of an activation function which can be used to model complex relationships in the data when there are many neurons. In this work, and in many works, the Exponential Linear Unit (ELU) [155] is used which is specified by

[TABLE]

where $\alpha$ is a strictly positive hyperparameter that moderates the value at which the ELUs saturate for negative inputs, which are then set to unity. The action of this function is to shift negative inputs closer to unity while leaving positive ones unchanged. Through the entirety of the NN, the combinations of the linear transformations and activation functions produce a vast array of hyperparameters that specify the NN responses. These can be optimized through training when their values are optimized for a specific set of training data. For each iteration during the process of training, the difference between the NN predictions and the ground truth in the training data is characterized by a so-called loss function which should be minimized. In the NN implementations that follow, we adopt Adam’s algorithm [156] which employs an improvement of the traditional gradient descent method which has been modified to accelerate convergence.

The simplest of these loss functions is the absolute difference between the training and predicted outputs which is then summed for each redshift. This is called the L1 loss function, and is defined as

[TABLE]

where $\Upsilon_{\rm train}(z)$ and $\Upsilon_{\rm pred}(z)$ are observed and NN predicted values of cosmological parameter $\Upsilon$ at observation redshifts $z$. This is the most similar to the log-likelihood function for uncorrelated data for MCMC methods. Other popular choices of loss functions include the mean squared error (MSE) loss function, which minimizes the squared differences between the training and predicted outputs, as well as the smoothed L1 (SL1) loss function, which softens the L1 absolute difference for similar values of outputs. It is at the level of the loss function that the complexities inherent in the training data are adopted by the NN. In our work, we adopt an L1 loss function since it is the most versatile and has shown itself to be the most efficient in training NN architectures in the cosmological regime.

Trained NN can be useful in reconstructing the evolution of various cosmological parameters such as reconstructing the expansion velocity and growth of large scale structures as in Refs. [157, 158, 70]. This removes the reliance on a fiducial model to extract points that do not coincide with observational redshifts. In Ref. [71], NN methods were used to test whether the absolute magnitude of type 1a supernovae had a significant aspect of evolution. NN implementations were also used to reconstruct specific classes of cosmological models such as in Refs. [159, 69] in which a Monte Carlo approach was taken for the derivatives of the reconstructed parameters. Recently, there has also been work on addressing the known problem of incorporating covariance matrices into the training process of NN implementations as in Ref. [160].

ANN training can be performed with cosmo_learn in two lines, initializing the architecture for each data set and training using a default architecture:

my_cl.init_ann()
my_cl.train_ann()

The results are shown in Figures 9 and 10. The convergence of the Loss function does appear to be slow for the individual data sets, which may be caused by several factors including the nuances of the underlying data, the complexity of the calculated cosmological parameters, as well as the configuration of the Universe in question. On the other hand, the convergence rate is not far beyond other comparable universe settings.

In our implementation, we employ the neural network package ReFANN [78] for function reconstruction. As with GP and BRR approaches, each data set is treated independently. For all reconstructions, we adopt the default ReFANN architecture, which consists of a single hidden layer with 4096 neurons and uses the hyperparameter configuration ‘rec_1’, characterized by the absence of batch normalization and dropout; use the option hp_model=‘rec_2’ for batch normalization. The number of epochs (default is 30000) used for training can be changed using the iteration option.

All models are compiled using the Adam optimizer [156] and trained with the L1 loss function. The regularization terms of ReFANN were specifically optimized for cosmological function reconstruction [78], and we found that no further tuning was required to reach our illustrative purposes. However, the effect of regularization terms on reconstruction is non-trivial and may be data- and probe-dependent. This is an avenue for data-driven exploration through cosmo_learn, with results potentially deserving of a dedicated, systematic treatment, which we reserve for a separate work.

For the SNe data set, as in the GP and BRR methods, redshifts are logarithmically rescaled before training the ANN to better capture features across a broad dynamic range. The full architecture used for each reconstruction can be printed upon initialization by setting the option show_summary=‘True’. To display training outputs, the option print_info=‘True’ can be used during initialization.

After training, the final architectures along with the evolution of training losses are stored in a dictionary accessible via my_cl.ANN_dict. For instance, the method show_ann_loss (which produced Figure 9) retrieves the training history from this dictionary for each data set. This modular structure enables flexible integration of heterogeneous cosmological data sets within a unified ANN-based reconstruction framework.

The results of ANN reconstruction on the reference mock data is shown in Figures 9-10. The ANN architecture can be customized for specific purposes by modifying the dictionary ANN_dict after initialization via init_ANN and before training via train_ANN.

The package cosmo_learn internally leverages reFANN [78], which is built on pytorch [79], to carry out ANN training and cosmological reconstructions.

4.6 Benchmarks

A timely effort afforded by a singular unifying framework for cosmological simulation, inference and reconstruction, such as cosmo-learn, is comparing the performance of old and each new statistical framework that is introduced into cosmological analysis. This is work in progress. Nonetheless, we can assess the performance of the different methods in cosmo-learn with simplified metrics, based on the computation time and the residuals obtained after subtracting the best fit or mean reconstruction to the data points. For this we rely to the SNe Ia mock data set as the densest of our data sets and the reference universe 14000605; this avoids having to complicate this initial benchmark with dealing with the residuals carrying different units for different observables. The result is presented in Table 2.

This shows that the effectiveness of each method in minimizing the residuals during training, but also when predicting cosmological realizations in the test set. Firstly, the training times of the different methods vary significantly; in order of increasing computation speed, ANN takes the most time, GP and MCMC take up roughly the same time, GA-Fisher, and BRR. We emphasize that these were obtained with the default configurations of each method in cosmo-learn and that the performance of each can be improved with further tuning hyperparameters. In MCMC this depends mostly on the number of dimensions of the parameter space, in GA-Fisher the mutation rate and number of function evaluations, in GP mostly on the number of data points since matrix inversion is an ${\cal O}(N^{3})$ procedure. On the other hand, ANN’s training efficiency depends on its architecture (number of layers, neurons, loss function) and BRR on the degree and type of basis functions used to construct the design matrix. BRRs speed advantage is derived through analytic tractability: hyperparameters are iteratively estimated with expectation-maximization and the computational cost is ${\cal O}(Np^{2}+p^{3})$ for $p$ basis functions (with $p=4$ utilized as the default in cosmo-learn’s BRR module). Clearly, scalability is an advantage for analyzing the rapidly increasing space of cosmological data.

However, scalability is only next in priority to the actual performance in terms of reconstruction of a statistical method. The mean, RMS and reduced chi-squared values in Table 2 show us a glimpse of the reconstruction performances reflected in residuals with respect to the SNe Ia mock data in the reference universe. Overall, all methods achieve comparatively small residuals mean and RMS, and $\chi^{2}_{\rm red}$ of order unity on both the training and test sets, indicating statistically consistent reconstructions. The performance of ANN in terms of the reduced chi-squared metric is noteworthy. It is worth noticing that the generally negative mean obtained across methods is a feature of the underlying realization (in the reference seed 14000605); in any realization, it is expected that the data may sit a hair above or below the fiducial curves. A rigorous comparison of these methods should be carried out across multiple mock realizations and observational probes, using a set of unbiased metrics; this is left for future work.

5 Conclusions

In this paper, we have introduced the python-based package cosmo_learn, a flexible and user-friendly toolkit for simulating cosmological data and performing data-driven inference. The software supports realistic noise modeling for various cosmological probes, including cosmic chronometers, supernovae, baryon acoustic oscillations, redshift space distortions, and gravitational waves. We have verified the consistency of the simulated data with the input cosmological model through residuals and parameter inferences. Additionally, we have demonstrated the core functionalities of several built-in learning and inference tools, including MCMC, genetic algorithm, Gaussian processes, Bayesian ridge regression, and artificial neural networks. Emphasis was placed on maintaining simplicity and adaptability throughout the development process to offer users an intuitive and extensible platform for exploring cosmology through simulation and machine learning.

Future developments of the package aim to incorporate additional cosmological data sets, such as cosmic microwave background measurements and big bang nucleosynthesis priors, to allow for more consistent and comprehensive modeling of initial conditions. Existing modules will be updated as improved data or models become available. While the current implementation uses a constant dark energy equation of state for which both background evolution and linear perturbations have analytic solutions, the framework can be readily extended to support dynamical dark energy models [161, 162, 163, 164]. The current learning framework is also designed to be extensible, making it straightforward to implement methods such as Approximate Bayesian Computation [126, 165, 166] and other emerging techniques that offer complementary statistical perspectives on cosmological tensions. We encourage contributions from the community and welcome modifications, extensions, and forks of our open-source codebase, which will remain publicly available and free to use.

With cosmo_learn, we hope to support both teaching and research efforts at the interface of cosmology, simulation, and machine learning.

\bmhead

Acknowledgements

RCB was supported by an appointment to the JRG Program at the APCTP through the Science and Technology Promotion Fund and Lottery Fund of the Korean Government, and was also supported by the Korean Local Governments in Gyeongsangbuk-do Province and Pohang City. RCB thanks the Institute of Physics, Academia Sinica, Taiwan for the hospitality that enabled the completion of this work. This article is also based upon work from COST Action CA21136 Addressing observational tensions in cosmology with systematics and fundamental physics (CosmoVerse) supported by COST (European Cooperation in Science and Technology). DG and VHC acknowledge partial support from Centro de Física Teórica de Valparaíso CEFITeV. JLS would also like to acknowledge funding from “Xjenza Malta” as part of the “Technology Development Programme” DTP-2024-014 (Cosmic Learning) Project. RCB and GCB thank participants of the external program APCTP-GW2025 [or APCTP2025-E05] held at Academia Sinica, Taipei, Taiwan for fruitful discussions.

\bmhead

Data availability No new data were generated in this work.

\bmhead

Code availability The python code cosmo-learn is publicly available in GitHub: https://github.com/reggiebernardo/cosmo_learn.