MadEvolve: Evolutionary Optimization of Cosmological Algorithms with Large Language Models

Tianyi Li; Shihui Zang; Moritz Münchmeyer

astro-ph.COFeb 17, 2026

MadEvolve: Evolutionary Optimization of Cosmological Algorithms with Large Language Models

Authors: Tianyi Li, Shihui Zang, Moritz Münchmeyer

Organizations: Department of Physics, University of Wisconsin-Madison

Abstract

We develop a general framework to discover scientific algorithms and apply it to three problems in computational cosmology. Our code, MadEvolve, is similar to Google's AlphaEvolve, but places a stronger emphasis on free parameters and their optimization. Our code starts with a baseline human algorithm implementation, and then optimizes its performance metrics by making iterative changes to its code. As a further convenient feature, MadEvolve automatically generates a report that compares the input algorithm with the evolved algorithm, describes the algorithmic innovations and lists the free parameters and their function. Our code supports both auto-differentiable, gradient-based parameter optimization and gradient-free optimization methods. We apply MadEvolve to the reconstruction of cosmological initial conditions, 21cm foreground contamination reconstruction and effective baryonic physics in N-body simulations. In all cases, we find substantial improvements over the base algorithm. We make MadEvolve and our three tasks publicly available at madevolve.org.

Figures & tables

Figure 1: Overview of the LLM-driven evolutionary pipeline in MadEvolve. The cycle proceeds as a closed loop: the Prompt Sampler retrieves parent programs from the Program DB to query the LLM Ensemble . The generated code diffs first undergo Parameter Optimization to refine continuous variables before being submitted to the Evaluator Pool . Evaluation scores update the database, completing the iteration. Finally, the best solution is transmitted to the central Report Generator to produce human-readable explanations of the discovered innovations.

Application	Baseline	Evolved	Improvement	Generations
BAO Reconstruction with StdRec ( $\bar{r}_{\text{BAO}}$ )	0.752	0.924	22.8%	1,165
BAO Reconstruction with IterRec ( $\bar{r}_{\text{BAO}}$ )	0.933	0.959	2.8%	386
21 cm Tidal Reconstruction ( $\bar{r}_{2D}$ )	0.743	0.971	30.7%	342
N-body baryonic displacements (test loss)	0.613	0.230	63%	233

Table 1: Summary of performance improvements across four evolutionary runs on three cosmological applications. BAO and tidal reconstruction metrics are cross-correlation coefficients (higher is better); tSZ uses L1 loss (lower is better). All evolved results represent the best algorithm discovered during the evolutionary run.

Metric	BAO	Iter. Recon.	Tidal	tSZ
Total generations	1,165	386	342	233
Successful programs	897	323	293	127
Success rate (%)	77.0	83.7	85.7	54.5
Improvement rate (%)	30.4	29.8	24.8	13.4

Table 2: Global statistics across the four evolutionary runs. Success rate indicates the fraction of LLM-generated programs that compile and execute without errors. Improvement rate indicates the fraction of all mutations that produce a valid program with higher fitness than its parent.

Model	$N$	Success (%)	Impr. (%)	Muts/Impr.
Gemini 3 Pro Preview	619	77.9	28.9	3.5
GPT-5.2	303	76.2	29.7	3.4
o4-mini	320	65.9	21.9	4.6
Gemini 2.5 Pro	291	82.1	28.9	3.5
Gemini 3 Flash Preview	228	89.5	32.9	3.0
GPT-5	186	84.4	20.4	4.9

Table 3: Per-model mutation statistics aggregated across all four applications.

N

: total mutations generated; Success: fraction producing valid, executable programs; Impr.: fraction of all mutations yielding a fitness score exceeding the parent’s; Muts/Impr.: average number of mutations needed to produce one fitness-improving program. Frontier models: Gemini 3 Pro, Gemini 2.5 Pro, GPT-5, GPT-5.2. Lightweight models: Gemini 3 Flash, Gemini 2.5 Flash, o4-mini. Not all models were used in every run.

Parameter	Value
Box size	$L=1000$ $h^{-1}$ Mpc
Particle count	$N_{p}=512^{3}$
Mass resolution	$m_{p}=6.56\times 10^{11}$ $h^{-1}M_{\odot}$
Grid resolution	$N_{\text{mesh}}=256^{3}$
Initial / final redshift	$z_{\text{IC}}=127$ , $z_{\text{final}}=0$
Cosmology	$\Omega_{m}=0.3175$ , $h=0.6711$ , $\sigma_{8}=0.834$

Table 4: Quijote simulation parameters for BAO reconstruction evaluation.

Figure 2: Evolution dynamics of BAO reconstruction algorithm discovery over 1,165 generations. The blue curve tracks the best training-set score

\bar{r}_{\text{BAO}}

achieved by any program in the population, while the green curve indicates the number of autodiff-optimizable parameters in the best-performing algorithm. The red dashed line marks the initial baseline performance (

\bar{r}_{\text{BAO}}=0.752

) achieved by the two-parameter Zel’dovich method. The evolutionary trajectory exhibits distinct phases: rapid initial improvement during the first 100 generations, a gradual increase through

\bar{r}_{\text{BAO}}\approx 0.85

–

0.88

as parameter count stabilizes near 10, a notable jump around generation 600–650, and continued refinement toward the final best score of

\bar{r}_{\text{BAO}}=0.924

found at generation 1,144. The +22.8% improvement demonstrates sustained evolutionary progress through the combination of LLM-proposed structural modifications and autodiff parameter optimization.

Algorithm	$\bar{r}_{\text{BAO}}$	Correlation	Runtime	Parameters
	( $k\in[0.01,0.5]$ )	(real space)	(seconds)
Baseline Zel’dovich (autodiff)	0.752	0.344	9.7	2
Evolved: Hybrid Augmented	0.924	0.581	76.1	10
Improvement	+22.8%	+68.9%	–	–

Table 5: Performance comparison of BAO reconstruction algorithms evaluated on 9 held-out test simulations.

Figure 3: A comparison of reconstruction performance for different algorithms. The results are averaged over nine simulations in the test set ( fiducial 1-9 ). The curves shows the cross-correlation function of dark matter field (black), standard reconstruction result (red), iterative reconstruction result (green), evolved algorithm result based on standard reconstruction (blue), evolved algorithm result based on iterative reconstruction (orange), and non-linear reconstruction result (purple) with respect to the initial condition. Left : Reconstruction is executed on a

256^{3}

mesh (which we evolved on). Right : Reconstruction is executed on a

512^{3}

mesh (which we did not evolve on). The shaded region corresponds to Fourier modes beyond the Nyquist frequency in the evolving configuration.

Figure 4: Evolution dynamics for the BAO reconstruction run starting from the iterative reconstruction baseline [ 37 ] . The blue curve tracks the best

\bar{r}_{\text{BAO}}

achieved by any program in the population across 386 generations. Starting from

\bar{r}_{\text{BAO}}=0.933

, the evolutionary search discovers differentiable refinements that push performance to

\bar{r}_{\text{BAO}}=0.959

(

+2.8\%

), demonstrating that the framework can improve upon an already well-optimized reconstruction algorithm.

Algorithm	$\bar{r}_{\text{BAO}}$	Correlation	Runtime
	( $k\in[0.01,0.5]$ )	(real space)	(seconds)
Iterative Reconstruction	$0.933$	$\mathbf{0.660}$	32.7
Evolved: Iterative Augmented	$\mathbf{0.959}$	$0.639$	37.1
Improvement	+2.8%	$-3.3\%$	+13.4%

Table 6: Performance comparison of evolved algorithm against iterative BAO reconstruction, evaluated on 9 held-out test simulations.

Figure 5: The cross-correlation coefficient $r(k_{\perp},k_{\parallel})$ of the intensity map with respect to the underlying matter field. The results are averaged over nine simulations in the test set ( fiducial 1-9 ).

Figure 6: Evolution dynamics of tidal reconstruction algorithm discovery over 342 generations. The blue curve tracks the best combined score

\bar{r}_{2D}

achieved by any program in the population, while the green curve indicates the number of tunable parameters in the best-performing algorithm. The red dashed line marks the initial baseline performance (

\bar{r}_{2D}=0.743

) achieved by the single-parameter tidal tensor method. The evolutionary trajectory exhibits rapid initial improvement as the system discovers Anisotropic Filtering , Generalized Anisotropic Potential , and Split-Weight Tensor Projection . By generation 52, the algorithm achieves

\bar{r}_{2D}=0.971

with only 9 parameters, capturing the essential physics of tidal reconstruction.

Algorithm	$\bar{r}_{2D}$	Runtime	Parameters
	(wedge region)	(seconds)
Baseline Tidal Tensor	0.743	38.4	1
Evolved: Split-Weight Anisotropic	0.971	81.0	9
Improvement	+30.7%	–	–

Table 7: Performance comparison of reconstruction algorithms.

Figure 7: A comparison of reconstruction performance for different algorithms. The results are averaged over nine simulations in the test set ( fiducial 1-9 ). Left : 1D cross-correlation coefficients

r(k)

. The red curve is the result from tidal reconstruction (baseline algorithm), the green curve is for the cluster fossil algorithm and the blue curve represents the performance of our new method. Right : Cross-Correlation coefficients binned into different directions.

\mu=0.83

is the direction closest to the line-of-sight, while

\mu=0.13

is the farthest.

Figure 8: Evolution dynamics of LDL tSZ prediction algorithm discovery over 233 generations. The plot shows progression of the cross-correlation coefficient

\bar{r}(k)

averaged over validation simulations. The evolutionary trajectory exhibits rapid initial improvement followed by a plateau and exploration phase, with a key breakthrough at generation 122 when the hybrid thermo-virial formulation is discovered. We discuss the physical interpretation of the discovered innovations in the next section.

Algorithm	Train Loss	Test Loss	$\bar{r}(k)$	Gen.
		(6-sim avg)	(all scales)	Found
Baseline LDL	0.371	0.613	0.943	0
Evolved: Physics-Informed	0.149	0.230	0.969	222
Improvement	60%	63%	+2.8%	–

Table 8: Performance comparison of tSZ prediction algorithms on CAMELS simulations (

64^{3}

mesh,

L=25\,h^{-1}\text{Mpc}

). Test metrics are averaged over six independent realizations (CV_5–CV_10) to assess generalization.

Figure 9: Scale-dependent performance comparison of baseline and evolved tSZ prediction algorithms. Top left : Cross-correlation $r(k)$ showing improved coherence across all scales. Top right : Transfer function $T(k)$ demonstrating better amplitude calibration. Bottom left : Power spectra comparison on log-log scale. Bottom right : Per-bin $r(k)$ improvement, with green indicating where the evolved model outperforms baseline. Results averaged over six test simulations.

Automated Algorithm Discovery

Automated algorithm discovery develops methods that generate and test executable algorithms for optimization, science, software, and machine learning.

See all 170 papers

Beyond AI as Assistants: Toward Autonomous Discovery in Cosmology
May 14, 2026Licong Xu, Thomas BorrettMulti-Agent LLM Systems AI Agents for Scientific Discovery
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK · Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK · Cavendish Astrophysics, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
MadEvolve: Evolutionary Optimization of Trading Systems with Large Language Models
May 21, 2026Yurii Kvasiuk, Tianyi Li, Owen Colegrove +1Quantitative Finance Evolutionary Optimization
Department of Physics, University of Wisconsin–Madison · Event Horizon Labs
ImprovEvolve: Basin-Hopping Meets LLM-Guided Evolutionary Search
Feb 10, 2026Alexey Kravatskiy, Valentin Khrulkov, Ivan OseledetsLarge Language Model-Guided Scientific Discovery Evolutionary Optimization
1MIRIAI, Russia · 2FusionBrain Lab, Russia · 3AXXX, Russia +1

May 14, 2026J/K move · Enter open · S save

Beyond AI as Assistants: Toward Autonomous Discovery in Cosmology

Licong Xu, Thomas Borrett

Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK · Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK · Cavendish Astrophysics, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

Multi-Agent LLM Systems AI Agents for Scientific Discovery Autonomous Scientific Discovery

Open paper AskPDF arXiv

Recent advances in artificial intelligence (AI) agents are pushing AI beyond tools toward autonomous scientific discovery. We discuss two complementary agentic systems for cosmology: CMBEvolve, which targets tasks with explicit quantitative objectives through LLM-guided code evolution and tree search, and CosmoEvolve, which targets open-ended scientific workflows through a virtual multi-agent research laboratory. As preliminary demonstrations, we apply CMBEvolve to out-of-distribution detection in weak-lensing maps, where it iteratively improves the benchmark score through code evolution, and CosmoEvolve to autonomous ACT DR6 data analysis, where it identifies non-trivial pair- and scale-dependent behaviour and produces analysis-grade diagnostics. These examples show how cosmology can provide both controlled benchmark tasks and realistic open-ended research problems for the development of AI scientist systems.

Abstract

Automated Algorithm Discovery

Explore similar work

Beyond AI as Assistants: Toward Autonomous Discovery in Cosmology

MadEvolve: Evolutionary Optimization of Trading Systems with Large Language Models

ImprovEvolve: Basin-Hopping Meets LLM-Guided Evolutionary Search