Cosmological Parameter Estimation with Genetic Algorithms

Abstract

Genetic algorithms are a powerful tool in optimization for single and multimodal functions. This paper provides an overview of their fundamentals with some analytical examples. In addition, we explore how they can be used as a parameter estimation tool in cosmological models to maximize the likelihood function, complementing the analysis with the traditional Markov chain Monte Carlo methods. We analyze that genetic algorithms provide fast estimates by focusing on maximizing the likelihood function, although they cannot provide confidence regions with the same statistical meaning as Bayesian approaches. Moreover, we show that implementing sharing and niching techniques ensures an effective exploration of the parameter space, even in the presence of local optima, always helping to find the global optima. This approach is invaluable in the cosmological context, where an exhaustive space exploration of parameters is essential. We use dark energy models to exemplify the use of genetic algorithms in cosmological parameter estimation, including a multimodal problem, and we also show how to use the output of a genetic algorithm to obtain derived cosmological functions. This paper concludes that genetic algorithms are a handy tool within cosmological data analysis, without replacing the traditional Bayesian methods but providing different advantages.

1. Introduction

Genetic algorithms (GAs), established for decades, are tools from evolutionary computation [1,2,3,4,5] that solve many function optimization problems. Evolutionary computation is focused on algorithms exploiting randomness to solve search and optimization problems using operations inspired by natural evolution [6]. It includes several methods for stochastic or metaheuristic optimization [7,8]; notable examples are Particle Swarm Optimization (PSO) [9] based on the social behavior of organisms of the same species such as birds, the Giant Trevally Optimizer (GTO) [10,11,12] inspired by the hunting behavior of predatory fish, and Artificial Rabbits Optimization (ARO), drawing inspiration from social interactions among rabbits [13,14]. Within evolutionary computation, the most relevant methods are genetic algorithms [15,16], genetic programming [17], and evolutionary strategies [18]; their success is due to their ability to navigate intricate, non-linear, and high-dimensional search spaces.

In particular, genetic algorithms stand out as powerful tools for optimization problems because they mathematically always guarantee, under certain conditions, to find the best solution. Despite challenges posed by local optimum values [19], this property puts them at an advantage over other techniques. Rooted in the emulation of natural selection and evolution, the iterative process of GAs involves generating a population, subjecting it to fitness-based selection, and applying genetic operators, such as crossover and mutation. This iterative approach drives the evolution of increasingly optimal solutions over generations. GAs thrive in situations with multiple optima, irregular landscapes, or where an analytical solution is difficult to achieve. Its adaptability allows for the simultaneous exploration of numerous candidate solutions, making them effective in various optimization challenges. Unlike traditional optimization methods, GAs have the advantage of not relying on derivatives, providing excellent robustness in high-dimensional or more complex problems. Inspired by natural evolution, these algorithms efficiently explore vast and unknown search spaces [20]. Their ability to solve complex and dynamic projects makes them valuable in diverse fields, including medicine [21,22,23], epidemic dynamical systems [24,25], geotechnics [26], market forecasts [27], and industry [28], among others. A particularly successful application in the Deep Learning era is the optimization of neural networks, which are huge computational models in which genetic algorithms help to find optimal combinations of hyperparameters [29,30,31].

With the accelerated development of computational resources, genetic algorithms and other machine learning algorithms have been exploited in several scientific fields in recent years. Remarkably, they have resulted in significant advances in understanding particle physics [32,33,34], astronomical information [35,36,37,38], and cosmological phenomena [39,40,41,42,43,44].

Genetic programming, another method from evolutionary computation, has been widely used in astrophysics and cosmology [45,46,47,48,49,50], which allows for symbolic regression for a given dataset, treating regression as a search problem to find the best combination of mathematical operators generating an expression fitting the data. Although genetic programming and genetic algorithms solve different tasks, they use similar operators to find solutions. In this work, we focus on genetic algorithms, mentioning genetic programming for reference, assuming the astrophysical community may be more familiar with it. Moreover, genetic algorithms are the most fundamental and successful evolutionary computation technique, and understanding them is useful for studying other evolutionary computation methods, including genetic programming.

On the other hand, parameter estimation in cosmology is a very relevant task that finds a combination of values for parameters describing a cosmological model based on observational data. The goal is to refine theoretical models to align with observations for a more precise understanding of the universe. In cosmological parameter estimation, the most robust and successful algorithms are Markov chain Monte Carlo; however, these methods sometimes are computationally expensive, and recent advancements try to attack this issue with new statistical or machine learning techniques, including the iterative Gaussian emulation method [51], adaptive importance sampling, parallelizable Bayesian algorithms [52], Bayesian inference accelerated with machine learning [53,54,55], or likelihood-free methods [56,57].

This paper aims to achieve two primary objectives: firstly, to provide a comprehensive introduction to genetic algorithms and elucidate their application in cosmological parameter estimation, and secondly, to demonstrate the complementarity of GAs with traditional Bayesian inference methods. We include illustrative examples of optimization problems and their applications in cosmology. Particularly, we delve into using genetic algorithms to constrain the parameter space of dark energy models based on observational data. It is pertinent to mention that GAs cannot perform the same tasks as MCMC methods, and we do not try to replace them; we only perform parameter estimation with GAs by optimizing the likelihood function, whereas MCMC methods sample the posterior probability function. However, we analyze their relevance as an alternative and complementary method, as discussed in Section 4.1.

The structure of this paper is as follows: In Section 2, we present the basics of genetic algorithms and an insight into their functionality. In Section 3, we provide some examples of the optimization of analytical functions by applying genetic algorithms. Section 4.1 describes the path to perform cosmological parameter estimation using these algorithms. Section 4.2 contains examples of multimodal problems in cosmology, and in Section 4.3, we justify how to obtain cosmological-derived parameters from a likelihood optimization. Finally, Section 5 summarizes our final remarks.

2. Fundamentals of Genetic Algorithms

2.1. Biological Fundamentals

Bioinspired computing is a field of computer science based on observing and imitating natural processes and phenomena to develop algorithms and computational systems [58]. These algorithms seek to solve complex problems. The bioinspired computation is classified into three main categories [58]: evolutionary algorithms (such as genetic algorithms), particle swarm intelligence (imitating collective behaviors) [7,59,60,61], and computational ecology (inspired by ecological phenomena) [8,62].

Genetic algorithms solve optimization [1,2,3,4,5] and search problems inspired by fundamental concepts of genetics and evolution [8,63,64]; some of its key points are as follows:

• Natural selection—This is the central principle in the theory of evolution. Just as better-adapted organisms are more likely to survive and reproduce in nature, GAs favor the fittest or most promising solutions from a population of candidate solutions. In nature, over several generations, the most promising characteristics of individuals survive to be inherited by the new generations. This is what genetic algorithms seek to do to have better solutions as more generations pass by.
• Crossing—Also called recombination, it is a process in which genes from two parents are combined to create offspring with characteristics inherited from both parents. GAs apply the idea of crossover by combining partial solutions from two individuals in the population to generate new solutions that can inherit desirable characteristics from both parents.
• Mutation—A mutation is recognized as the stochastic alterations in an organism’s genetic material. In the GAs, a mutation introduces random changes in a small part of the candidate solutions, e.g., it may change the value of a bit, which increases the diversity of possible solutions and improves the exploration of the search space.
• Reproduction and inheritance—In the same sense as in nature, in genetic algorithms, these operations allow for the transmission of some characteristics of the parent solutions to the solutions of the next generation (offspring).

2.2. Genetic Algorithm Operations

John Holland was the first to introduce the genetic algorithm in 1975 in his book Adaptation in Natural and Artificial Systems [3,15]. According to the GA context, a population is a set of possible solutions to a given problem. Each individual has a genotype encoded in bits, which is then expressed as a phenotype in the problem context. The way to encode the possible solutions is fundamental to attacking a problem with GAs, and there are several options to perform it, for example, with binary, integer, or real encoding, among others [65].

Alternatively, assessing an individual’s quality or a potential solution involves employing a metric or target function, which is ideally expected to approach its optimal value in the final generations. For the analogy of natural selection, this target function, or objective function, is called the fitness function. In practice, in GAs, the fitness function is directly the function to be optimized. This is unlike genetic programming, where the fitness function is a measure of the error between the algebraic expression found and the dataset used due to the regression task that genetic programming addresses.

The continuous evaluation of all the individuals (possible solutions) of a population with this fitness function and the applications of genetic operations to produce new generations allow for GAs to find the optimal value of this function. In the following list, we describe the fundamental procedures of genetic algorithms [66]:

• Selection—It is the method of choosing the best solutions to play the role of parents and improve the quality of the offspring. Several selection methods include the roulette [67], random [68], ranking [69], tournament [70], and Boltzmann entropy selections [71].
• Crossover—It is also called recombination, which generates a new possible solution given two previously selected parents. There are several crossover methods, such as one point, two points, N points, uniform, three parents, random, and order. The crossover operation has an associated probability ($P_c$) that determines how many individuals recombine given the population, with $P_c=1$ indicating that all the products come from the recombination and $P_c=0$, meaning they are exact copies of the parents.
• Mutation—After crossover, mutations make it possible to maintain diversity in the population and prevent it from stagnating at the local optima [72]. There are several types of mutation operators, such as flipping a gene if it is in the same position as in the parent; swapping values at random positions; flipping values from left to right, or in a random sequence; and shuffling random positions. Mutation also has a probability associated with it that indicates how likely it is to randomly change a gene (bit) of a possible solution. The mutation value must be low for an efficient search within the genetic algorithm1.
• Replacement—The last step is the replacement, which keeps the population size constant by eliminating individuals after recombination. There are three methods: strong replacement (random), weak replacement (the two fittest), and replacing both parents (the children replace both parents).
• Elitism and Hall-of-Fame—The elitism method ensures that the best individuals are not discarded but transferred directly to the next generation. Hall-of-Fame is an integer that indicates how many individuals are considered under elitism to be retained in the next generation. Elitism is necessary to ensure that genetic algorithms always find the best solution [19]. Elitism and Hall-of-Fame are often considered distinct from the general replacement strategy. While the replacement strategy primarily focuses on selecting individuals for reproduction and forming the next generation, the elitism and Hall-of-Fame mechanisms specifically address preserving the best-performing individuals.
• Stopping criteria—A mechanism is needed to finalize the execution of the genetic algorithm. Some ways to perform it are to stop after a fixed number of generations, after a specific time-lapse, to finish the process if the best fitness does not change for several generations (steady fitness), or to stop it if there are no improvements in the objective function for several consecutive generations (generation stagnation).

In this way, we can summarize that genetic algorithms are a process that involves some crucial steps: the initialization of a population form of solutions, selection of parents according to their fitness, recombination of genes by crossing, introduction of variability by mutation, substitution of individuals, and running the algorithm until the stopping criterion is satisfied. The operations described above are repeated within a loop, generation after generation, until a satisfactory solution or convergence criterion is reached.

2.3. Schema Theorem

The heuristic search of genetic algorithms is based on Holland’s schema theorem, which states that the chromosomes have patterns called schemas. This schema theorem deals with the decomposition of chromosomes into schemas and their influence on the evolutionary dynamics of the population.

A schema is a binary string of fixed length representing a chromosome pattern. For example, in a chromosome of length 6, the schema 001X00 defines a string that starts with 001, has an unknown bit X, and ends with 00.

The fitness of a schema refers to how many individuals in the population contain that specific schema. It can be represented as a fitness function $F\left(S\right)$ that denotes the fitness of the schema S.

The schema theorem states that high-fitness schemas are more prevalent in future generations. This is because schemas with high fitness are more likely to be selected and recombined, leading to population improvement in terms of fitness. Mathematically, we can express this as:

$$F\left(S_{t+1}\right)\ge (1-p_m)·F\left(S_t\right),$$

(1)

where $F\left(S_{t+1}\right)$ is the fitness of the schema S at the next generation $(t+1)$, $F\left(S_t\right)$ is the fitness of the schema S in the current generation $\left(t\right)$, and finally, $p_m$ is the mutation probability.

This equation indicates that the fitness of the schema in the next generation is at least equal to the current fitness, modulated by the mutation probability. If $p_m$ is low, schemas with high fitness will likely survive and propagate in future generations, contributing to population improvement.

3. Genetic Algorithm Application

In this section, we implement a genetic algorithm to optimize univariate functions and extend its application to higher-dimensional problems. The general structure of a genetic algorithm is provided in the Algorithm 1.

Algorithm 1 Simple Genetic Algorithm

Parents ← {randomly generated population} While not (termination) Calculate the fitness of each parent in the population Children $\leftarrow \varnothing$ while |Children| < |Parents| Use fitness to probabilistically select a pair of parents for mating Mate the parents to create children $c_1$ and $c_2$ Children ← Children $\cup \{c_1,c_2\}$ Loop Randomly mutate some of the children Parents ← Children Next generation

Several libraries incorporate genetic algorithms, such as Distributed Evolutionary Algorithms (DEAP) [73], Karoo GP [74], Tiny Genetic Programming [75], and Symbiotic Bid-Based GP [76]. These libraries simplify the implementation of genetic algorithms. In this paper, we have utilized the DEAP library, which boasts comprehensive documentation.

3.1. Single Variable Functions

Considering the following three functions, we aim to use a custom genetic algorithm to find their global maxima:

• $f_1\left(x\right)=(x^2+x)cos\left(2x\right)+x^2$;
• $f_2\left(x\right)={sin}^2(3x+45)+0.9{sin}^3\left(9x\right)-sin(15x+50)cos(2x-30)$;
• $f_3\left(x\right)=-x^6/60-x^5/50+x^4/2+2x^3/3-3.2x^2-6.4x$.

In Figure 1, it can be seen how the above functions are optimized by a genetic algorithm, using a population size of 100 individuals, with a Hall-of-Fame size equal to 1, a mutation probability of 0.2, and a crossover probability of 0.5, over 50 generations. Note that as the generations progress, the individuals are closer to the global maxima. Another interesting feature is that, despite the local optima, the genetic algorithm in all the functions can find the global optima, as it is mentioned in the Introduction and Ref. [19].

3.2. Multimodal Functions

Genetic algorithms can also address problems with multiple dimensions and maxima by modifying the representation of candidate solutions and the operators used to generate new solutions. They can explore complex search spaces efficiently and identify global or local optima by appropriately designing crossover and mutation operators and analyzing different encoding techniques.

We use the Himmelblau function to demonstrate how genetic algorithms can be used to optimize these types of multimodal functions. We use the DEAP library, a robust Python framework for evolutionary computation, to achieve our goal. The following equation defines Himmelblau’s function:

$$f(x,y)={(x^2+y-11)}^2+{(x+y^2-7)}^2.$$

(2)

The niching and sharing technique is employed to identify all the global optima within a single genetic algorithm run. This concept draws inspiration from nature, where regions are divided into sub-environments or niches, enhancing population efficiency and survival. Individuals compete for resources in these niches independently of those in other niches. By integrating a sharing mechanism into the genetic algorithm, individuals are incentivized to explore new niches, discovering multiple optimal solutions, each considered as a niche. Typically, this is achieved by dividing an individual’s fitness value by the sum of distances from all other individuals. This approach penalizes overpopulated niches by distributing the local rewards among their individuals [77].

Niching involves dividing the population into subpopulations, each assigned to explore a specific region in the solution space. This encourages diversity by allowing genetically engineered individuals to compete for fitness locally. Conversely, sharing ensures a fair distribution of the fitness resources among individuals within the same niche. An individual’s fitness is influenced not only by its performance but also by the performance of its neighbors, preventing overemphasis on a specific region and promoting a balanced exploration. This approach prevents premature convergence to a local maximum, allowing for the simultaneous exploration of different regions and ultimately facilitating the identification of the global maximum. Applying this technique effectively requires a larger population size and more generations than a simple genetic algorithm. This is essential to spread the population across the sample space, targeting different niches and, consequently, identifying multiple optimal maxima. In our experiment, we executed the algorithm with 200 individuals and 200 generations, and the outcomes are summarised in Figure 2 and

Table 1. A comparison is made among the four real global optima of Himmelblau’s function [77] and those found by the genetic algorithm using niching and sharing.

Real Optimum	Optimum Found by GA
$(3.000,2.000)$	$(3.010,1.998)$
$(-2.805,3.131)$	$(-2.802,3.133)$
$(-3.779,-3.283)$	$(-3.774,-3.292)$
$(3.584,-1.848)$	$(3.585,-1.847)$

.

As can be seen in Table 1, these results are remarkably similar to the real values. Improving these results is possible by increasing the number of individuals and generations. It should also be noted that this technique is not limited to three dimensions but can be generalized to N dimensions and can support the search for global M optima. However, it is important to remember that as the number of dimensions increases, more computational resources are required to search effectively.

3.3. Statistical Analysis

Genetic algorithms are handy tools in statistical applications for optimizing likelihood functions, thereby determining the parameters of a scientific model (which is precisely what this article aims to demonstrate). However, reporting a confidence interval for the output of a genetic algorithm can be more complex than in classical statistical methods. The most rigorous technique relies on having a mathematical model of the genetic algorithm’s convergence that extends beyond Holland’s schema theory for the simple genetic algorithm published in 1975.

Because the state of the population in a genetic algorithm depends solely on the previous state in a probabilistic manner, Markov chains have been studied as suitable models for specific applications, and more recently, others have been modeled as martingales [78,79].

However, it is possible to resort to less rigorous techniques. One approach is to assume a distribution for the optimized parameters. For instance, assuming the parameters follow a normal distribution, the confidence interval can be calculated based on standard deviations, and the confidence ellipses can be computed using Fisher matrices. This is the procedure employed in this article. Another procedure involves using the Bootstrap method or other re-sampling techniques [80].

4. Application in Cosmology

In observational cosmology, one of the fundamental tasks is to determine the values of the free parameters for a given theoretical model based on observational measurements. This involves creating a function that captures the discrepancies between observed data and theoretical predictions and using it to obtain a parameter estimate that fits the data well. The likelihood function is typically used to represent the data’s conditional probability given the theory and its parameters. Although Bayesian inference is the most robust method for parameter estimation in cosmology, as it allows for sampling the posterior probability of parameters given the data, it can be computationally intensive (see the nomenclature under the Bayesian formalism of Bayes’ theorem [81,82]); instead of sampling the posterior probability function to estimate parameter values efficiently, optimization algorithms can be used to find the maximum likelihood function. In Reference [82], there is an exciting overview of the difference between sampling and optimization, and it can be seen that they are two different tasks that can be complementary. This section presents three applications that show how genetic algorithms can be applied to analyze cosmological data. First, we offer parameter estimation in three cosmological models: $\mathsf{\Lambda }$CDM, CPL, and PolyCDM. We then discuss how genetic algorithms can be used in a cosmological model with multiple maximum values, such as the graduated dark energy model presented in Ref. [83].

The datasets utilized in this section comprise 31 cosmic chronometers [84,85,86,87,88,89,90,91], Baryon Acoustic Oscillation measurements (BAO) [92,93,94,95,96,97], 1048 Type Ia supernovae (SNeIa) sourced from the Pantheon compilation [98], and binned data from the Joint Light Analysis SNeIa compilation [99].

Considering the datasets mentioned above, we employ the following log-likelihood functions for the Bayesian inference and optimization methods:

$$log{\mathcal{L}}_i=-\frac{1}{2}{(D_{\mathrm{th}}^i-D_{\mathrm{obs}}^i)}^T·C_i^{-1}·(D_{\mathrm{th}}^i-D_{\mathrm{obs}}^i),$$

(3)

where the index i ranges from 1 to 3, corresponding to the three datasets: cosmic chronometers [$D^{i=1}=H\left(z\right)$] and BAO [$D^{i=2}=D_A\left(z\right)$], where $D_A\left(z\right)$ represents the Hubble, volume averaged, and angular distance; and SNeIa [$D^{i=3}=\mu \left(z\right)$], where $\mu \left(z\right)$ denotes the distance modulus. In this context, $D_{\mathrm{obs}}$ represents the observed measurements, while $D_{\mathrm{th}}$ represents the theoretical values for the cosmological models. The matrices $C_i$ encompass the covariance information, accounting for systematic and statistical errors.

We implemented a module to work with the DEAP genetic algorithms within the SimpleMC2 code for our cosmological parameter estimation [100]. In some of the subsequent results, we compare the genetic algorithm’s outcomes with those of Bayesian inference obtained using the nested sampling algorithms, a specialized type of Markov chain Monte Carlo (MCMC) technique [81,101]. Additionally, we utilize the Fisher matrix formalism described in Refs. [102,103] to calculate the confidence intervals and generate error plots for the genetic algorithm-based parameter estimation. It is important to emphasize that genetic algorithms are not employed to generate posterior samples; instead, they are used to explore maximum likelihood estimation, which can yield similar and quicker results than parameter estimation. However, they cannot replace the robustness of MCMC methods. Furthermore, we conducted maximum likelihood estimation using a classical optimization method, specifically the L-BFGS algorithm [104], for comparison purposes and to assess the advantages of genetic algorithms.

4.1. Cosmological Parameter Estimation

As previously mentioned, we employ genetic algorithms to evaluate their effectiveness in parameter estimation. As a proof of the concept, and for simplicity, we consider three cosmological models, $\mathsf{\Lambda }$CDM, CPL, and PolyCDM, which are described below:

• $\mathsf{\Lambda }$CDM. The $\mathsf{\Lambda }$CDM model serves as the standard cosmological model and comprises two primary components: cold dark matter (CDM), which plays a pivotal role in the universe’s structure formation, and dark energy, which exhibits a counter-gravitational behavior, leading to the universe’s accelerated expansion. The cosmological constant, denoted by $\mathsf{\Lambda }$, is the simplest and most straightforward representation of dark energy, which exerts a pressure equal in magnitude but opposite in sign to the universe’s energy density ($p=-\rho$). For a flat universe in the late stages of its evolution, the equation governing its expansion is given by $H^2\equiv {\left(\frac{\dot{a}}{a}\right)}^2={\rho }_m\left(t\right)+{\rho }_{\mathsf{\Lambda }}\left(t\right)$, where a represents the scale factor, the dot denotes the derivative with respect to time, ${\rho }_m$ signifies the density of dark matter and baryons, and ${\rho }_{\mathsf{\Lambda }}$ accounts for the dark energy content in the form of a cosmological constant. These two parameters describe the evolution of the universe’s content. Incorporating their initial conditions denoted with a subscript 0, this equation can be re-expressed in terms of the redshift $1+z=1/a$ as follows:

  $$H^2=H_0^2[{\mathsf{\Omega }}_{\mathrm{CDM},0}{(1+z)}^3+{\mathsf{\Omega }}_{\mathsf{\Lambda },0}],$$

  (4)

  where $H_0$ denotes the Hubble constant, providing the present rate of expansion of the universe. The parameters ${\mathsf{\Omega }}_{\mathrm{CDM},0}$ and ${\mathsf{\Omega }}_{\mathsf{\Lambda },0}$ are specific to the $\mathsf{\Lambda }$CDM model. The former represents the current dimensionless density of dark matter (plus baryons), while the latter signifies the dimensionless density of dark energy. These parameters are subject to the constraint ${\mathsf{\Omega }}_{\mathrm{CDM},0}+{\mathsf{\Omega }}_{\mathsf{\Lambda },0}=1$; when this equality holds, we have a flat universe [105]. Consequently, for this model, we effectively have two free parameters, namely, h and ${\mathsf{\Omega }}_{\mathrm{CDM},0}$, which we simplify by denoting ${\mathsf{\Omega }}_{\mathrm{CDM}}$ as ${\mathsf{\Omega }}_m$ for brevity.
• CPL model. One can discern dark energy’s characteristics by investigating its state equation, denoted as $w\left(z\right)$, where p and $\rho$ represent the pressure and dark energy density, respectively [106]. Chevallier, Polarski, and Linder introduced the following parametrization for the equation of state, $w\left(z\right)=w_0+w_a\frac{z}{1+z}$, where $w_0$ signifies the current value of the equation of state. In contrast, $w_a$ represents its rate of change over time [106]. This equation of state leads to the following derivation:

  $$\begin{array}{cc}\hfill H{\left(z\right)}^2& =H_0^2[{\mathsf{\Omega }}_{m,0}{(1+z)}^3+\hfill \\ \hfill & (1-{\mathsf{\Omega }}_{m,0}){(1+z)}^{3(1+w_0+w_a)}{e}^{-\frac{3w_{a}z}{1+z}}].\hfill \end{array}$$

  (5)
  Now, the parameter estimation consists of finding the free parameters $H_0$, ${\mathsf{\Omega }}_{m,0}$, and $w_0$ and $w_a$.
• PolyCDM. We can consider an extension of dynamical dark energy by introducing spatial curvature, ${\mathsf{\Omega }}_1$, which adapts to the evolution of dark energy at low redshifts [41]. By performing a Taylor series expansion of the Equation (4) [107], we arrive at the PolyCDM model:

  $$\begin{array}{cc}\hfill H^2=& H_0^2({\mathsf{\Omega }}_{m,0}{(1+z)}^3+\hfill \\ & {\mathsf{\Omega }}_{1,0}{(1+z)}^2+{\mathsf{\Omega }}_{2,0}(1+z)\hfill \\ \hfill & +(1-{\mathsf{\Omega }}_{m,0}-{\mathsf{\Omega }}_{1,0}-{\mathsf{\Omega }}_{2,0})),\hfill \end{array}$$

  (6)

  where ${\mathsf{\Omega }}_{m,0}$ represents the dark matter; and baryon, contribution, and ${\mathsf{\Omega }}_{2,0}$ can be interpreted as the “lost matter” [107]. PolyCDM can be considered a parametrization of the Hubble parameter  [108].

For all the models mentioned above, we use a genetic algorithm with elitism, using 50 generations, a mutation probability of 0.2, a crossover probability of 0.7, a population comprising 100 individuals, and a Hall-of-Fame size of 2 to maximize the likelihood probability function.

Table 2. Parameter estimation via genetic algorithms for the $\mathsf{\Lambda }$CDM, CPL, and PolyCDM models utilizing cosmic chronometers, BAO, and SNeIa datasets. The $-2log\mathcal{L}$ value represents the optimal fitness value.

		Data: CC + BAO + SNeIa
Model	Parameters	L-BFGS Optimizer	Genetic	Nested
$\mathsf{\Lambda }$CDM	$h_0$	$0.6972\pm 0.0170$	$0.6964\pm 0.0170$	$0.6963\pm 0.0160$
	${\mathsf{\Omega }}_m$	$0.2950\pm 0.0133$	$0.2958\pm 0.0133$	$0.2960\pm 0.0134$
	$-2log\mathcal{L}$	$1049.2424$	$1049.2476$	$1049.2445$
CPL	$h_0$	$0.6864\pm 0.0259$	$0.6916\pm 0.0258$	$0.6901\pm 0.0240$
	${\mathsf{\Omega }}_m$	$0.2853\pm 0.0221$	$0.2919\pm 0.0218$	$0.2892\pm 0.0211$
	$w_0$	$-1.0082\pm 0.0840$	$-0.9803\pm 0.0912$	$-0.9909\pm 0.0861$
	$w_a$	$0.2556\pm 0.5188$	$0.0330\pm 0.6035$	$0.0679\pm 0.5296$
	$-2log\mathcal{L}$	$10483.9018$	$1049.0778$	$1048.9415$
PolyCDM	$h_0$	$0.6913\pm 0.0283$	$0.6916\pm 0.0283$	$0.6916\pm 0.0250$
	${\mathsf{\Omega }}_m$	$0.2899\pm 0.0290$	$0.2931\pm 0.0294$	$0.2945\pm 0.0198$
	${\mathsf{\Omega }}_{1,0}$	$0.0150\pm 0.4254$	$0.0947\pm 0.4271$	$0.1232\pm 0.1795$
	${\mathsf{\Omega }}_{2,0}$	$0.0136\pm 0.1995$	$-0.0147\pm 0.2007$	$-0.0298\pm 0.0903$
	${\mathsf{\Omega }}_k$	$-0.0013\pm 0.0703$	$-0.0076\pm 0.0702$	$-0.0004\pm 0.0117$
	$-2log\mathcal{L}$	$1049.0688$	$1049.0660$	$1049.1286$

and Figure 3 present the parameter estimation results obtained throughout the three methods outlined earlier. It is noticeable that, in most cases, the genetic algorithm results closely align with the parameter estimations derived from the nested sampling. Consequently, although they are slower than optimization methods like the L-BFGS method, genetic algorithms offer greater precision while remaining faster than MCMC algorithms. It is important to note that genetic algorithms maximize the likelihood function rather than sampling the posterior distribution. This distinction can be computationally advantageous compared to Bayesian inference procedures in specific scenarios. However, GAs lack the assignment of weights to individuals, as found in Bayesian inference samples, and their exploration of parameter space differs from MCMC methods, which rely on Markov chains and probabilistic conditions. Genetic algorithms, instead, focus on achieving improved solutions in each generation.

4.2. Multimodal Models

Parameter inference in some models can lead to the identification of multiple optima, meaning that posterior probability functions can have multimodal distributions. To address this complexity, Bayesian nested inference algorithms, such as multinest [109], are a sampling method designed to deal with multimodal distributions, allowing for the effective sampling of the parameter space. In contrast, classical optimization algorithms are limited to finding a single maximum. Genetic algorithms, thanks to niche and sharing techniques (see Section 3.2), have the ability to exhaustively explore the parameter space, even in the presence of local maxima. An example of a model with multiple maxima in its posterior distribution is the case of graduated dark energy [83], which is governed by the following Friedmann equation:

$$\begin{array}{cc}\hfill H^2=& H_0^2[{\mathsf{\Omega }}_{r,0}{(1+z)}^{-4}+{\mathsf{\Omega }}_{m,0}{(1+z)}^{-3}+\hfill \\ \hfill & {\mathsf{\Omega }}_{\mathrm{DE},0}{sgn[1-\psi lna]|1-\psi lna|}^{\frac{1}{1-\lambda }}],\hfill \end{array}$$

(7)

where ${\mathsf{\Omega }}_{\mathrm{DE},0}$ is the dimensionless density parameter of the dark energy with $\psi <0$ and $\lambda =0,-2,-4,\dots$. Also, $\psi$ is defined in terms of $\lambda$ and another parameter $\gamma$ in the following way: $\psi \equiv -3\gamma (\lambda -1)$. One maximum value corresponds to the $\mathsf{\Lambda }$CDM model, whereas the other is present to alleviate the Hubble tension. This model resembles a rapid transition of the universe from anti-de Sitter vacua to de Sitter vacua; see the details of the model in the references [83,110,111,112,113].

For the genetic algorithm with elitism used in this case, we set 20 generations, 200 individuals for the population, crossover and mutation probabilities of 0.5 and 0.2, respectively, and a Hall-of-Fame of size 2. Therefore, the free parameters for the graduated dark energy model are ${\mathsf{\Omega }}_{m,0}$, $h_0$, $\lambda$, and $\gamma$. For this example, to appreciate the multimodality in the graduated DE model, we use the same data as that in the original work (Ref. [83]), i.e., cosmic chronometers, BAO, and SNeIa (binned data from the Joint Light Analysis compilation [99]), but for simplicity, we do not use the Planck information. We also fix $\lambda =-20$. Performing Bayesian inference on this model, the posterior distribution for the $\gamma$ parameter is shown in Figure 4, in which two modes exist. In

Table 3. Parameter estimation with nested sampling (sampling the posterior probability distribution function), L-BFGS, and genetic algorithm. In these cases, we only consider the maximum likelihood found in the three methods and their corresponding parameter values.

	Nested Sampling	L-BFGS	Genetic
${\mathsf{\Omega }}_m$	0.3264	0.2991	0.2959
h	0.6947	0.6760	0.6765
$\gamma$	−0.0129	0.0000	−0.0127
$-2log\mathcal{L}$	55.8700	60.5781	61.6997

, we can analyze the outputs of the parameter estimation using nested sampling through posterior distribution sampling, the L-BFGS optimization method, and a genetic algorithm maximizing the likelihood distribution function; we can notice that the results maximizing the likelihoods are roughly consistent with the parameter estimation with Bayesian inference; however, for the $\gamma$ value, the L-BFGS method is unable to find a value different to zero, and it is far from the estimation of this parameter using the same data.

As mentioned above, some algorithms for Bayesian inference, such as multinest nested sampling, could explore the regions with these two maxima; however, most MCMC methods cannot achieve this task. Using genetic algorithms with the niching and sharing techniques, we can quickly find and explore the parameter space with these two optima without performing a Bayesian inference process; we can notice them in the histograms of Figure 5, in which the GAs explore the regions of both modes of the $\gamma$ parameter. Therefore, we can have more confidence in the results of a genetic algorithm than a classical optimization method.

To conclude this section, it is worth noting that there are other multimodal cosmological models, mainly involving neutrinos and spatial curvature, documented in the literature [114,115,116,117,118], and it is worth exploring in future works where these techniques could prove valuable for conducting efficient and rapid assessments.

4.3. Derived Functions

As an additional application, taking advantage of the genetic algorithms’ nature, we can use the saved individuals along generations to maximize the likelihood function and calculate the derived functions to analyze their phenomenological behavior. This technique is usually used with the samples of the posterior probability with Bayesian inference algorithms, mapping the sampling of an estimated parameter to another derived one. For example, the library fgivenx [119] allows for this mapping. In the case of the individuals of likelihood optimization using genetic algorithms, the statistical meaning of the plots is not directly related to the posterior probability function; however, it can provide an idea of the behavior of the derived functions given the estimated parameters.

In Figure 6, we compare the equation of state reconstructed from the outputs of Section 4.1 for the CPL model, and we use the samples for the $w_0$ and $w_a$ from the nested sampling and the values of these same parameters from the history of the individuals of the genetic algorithm population. We can notice that the behavior of the equation of state, analyzing the darkest regions, is similar in both cases, and it suggests that for a quick test, we can use this technique with genetic algorithms. Regarding the confidence regions, because we are only optimizing the likelihood function with the genetic algorithms, we cannot have a formal way to estimate them correctly.

5. Conclusions

In this study, we have leveraged genetic algorithms as an effective tool to estimate the free parameters of four cosmological models. Individuals generated in each genetic algorithm population have demonstrated the ability to achieve faster parameter estimates than those obtained using MCMC methods, thus reducing the number of likelihood function evaluations required. In addition, these genetic algorithms allow for a rapid computation of derived parameters, which adds flexibility and efficiency to the estimation process.

However, it is important to note that genetic algorithms differ from Bayesian approaches in their sampling process. While MCMC methods fully sample the posterior probability function, genetic algorithms focus on maximizing the likelihood function. This distinction implies that genetic algorithms cannot directly provide confidence regions with the same statistical significance as Bayesian inference procedures. However, they offer significant advantages, such as a faster speed and better results than other optimization methods, such as the L-BFGS algorithm.

Additionally, we have explored the usefulness of sharing and niche techniques in genetic algorithms, ensuring practical parameter space exploration, even in local or global optima. These features may be especially valuable in cosmology as a prior analysis to maximize the likelihood function before undertaking more computationally expensive Bayesian parameter estimation.

Throughout this paper, we can understand why genetic algorithms have been a very promising field of research over the last decades. Their flexibility allows for their application in diverse tasks, such as optimization, combinatorics, statistics, and even to speed up computational algorithms. The potential future applications of genetic algorithms in cosmological research are vast. With the presented study, we show the prospect of using them as a complement within cosmological data analysis. This is in agreement and complementary with the existing research that also focuses on the statistical applications of evolutionary computation [39,120]. In our case, we have not proposed a novel method or algorithm; however, we have analyzed how to use GAs so that they can complement a traditional analysis of cosmological data and be an alternative to optimize the likelihood function. We are convinced that genetic algorithms are a great technique with diverse cosmological and statistical applications. For example, in a parallel work, we have explored their usefulness to improve cosmological neural reconstructions [31] and to reduce the computational time of Bayesian inference routines. Therefore, we are confident that genetic algorithms are an excellent complementary element to the cosmological data analysis toolkit.