GAE: Graph-Augmented Evolution for Scientific Discovery via Reinforcement Optimization

# Graph-Augmented Evolution for Scientific Discovery via Reinforcement Optimization
Source: [https://arxiv.org/html/2607.10127](https://arxiv.org/html/2607.10127)
###### Abstract

Evolutionary program search guided by Large Language Models \(LLMs\) has emerged as a powerful paradigm for automated scientific discovery\. However, current approaches are fundamentally constrained by three bottlenecks: structurally blind parent selection, sparse whole\-program evaluation rewards, and static mutation operators that fail to adapt during search\. We presentGAE\(Graph\-AugmentedEvolution\), a framework that resolves these limitations through a tightly coupled, three\-pillar architecture\. First, arelational graph neural network \(GNN\)parses programs into typed computation graphs, producing structure\-aware embeddings\. Second, anRL\-optimized meta\-controllerleverages these embeddings to replace blind evolutionary sampling with a directed policy, dynamically selecting optimal parents and mutation directions based on reward history\. Third, anonline GRPO fine\-tuning loopcontinuously updates the LLM mutation operator at test\-time using group\-normalized evaluation rewards, directly aligning the model’s generation distribution with high\-fitness structural edits\. We evaluate GAE on a challenging scientific discovery task: symbolic regression for complex nonlinear oscillator systems\. By transforming stochastic search into a directed, self\-improving trajectory, GAE efficiently discovers closed\-form physical equations, consistently matching or outperforming static LLM\-driven baselines and achieving state\-of\-the\-art out\-of\-distribution performance\.

Machine Learning, ICML

## 1Introduction

The discovery of new algorithms and programs is a central frontier in AI\-assisted science\. Recent work follows a simple but powerful recipe: combine a large language model with programmatic evaluation and iterative search\. FunSearch\(romera2024mathematical\)first demonstrated that this loop can discover genuinely new mathematical objects, outperforming decades of hand\-crafted combinatorial search\. AlphaEvolve\(novikov2025alphaevolve\)generalized the principle to complete codebases, achieving state\-of\-the\-art results on hardware scheduling, matrix multiplication, and open problems in math and engineering\. Following this line of work, the overall procedure can be viewed as an iterative search framework in which an LLM serves as a program mutation operator that generates candidate code modifications, a task\-specific evaluator assigns fitness based on predefined scoring criteria, and a selection mechanism iteratively preserves and refines high\-performing programs within the search space\.

Despite these successes, three interrelated bottlenecks hinder effective evolutionary search, especially in reliably steering program mutations toward high\-quality regions of a large and complex program space\.\(1\) Reward sparsity\.Each fitness evaluation requires a complete program execution as a full training run in the case of program search, or a numerical integration in symbolic regression\. The archive is updated only when a child outperforms the current occupant of its MAP\-Elites cell, so the search process receives a single scalar signal indicating whether the candidate improves upon the incumbent in terms of the defined fitness score\. Programs that nearly beat an incumbent, for example achieving performance very close to the current best, or that contain partially correct or structurally meaningful sub\-expressions, are still discarded if they do not strictly improve the cell\. As a result, a large population of evaluated programs is reduced to a sparse sequence of binary improvement events, discarding fine\-grained relational information about proximity to improvement and the contribution of intermediate structural modifications\.\(2\) Uninformed parent selection\.AlphaEvolve samples parents by fitness rank with random diversity, treating the population as an unstructured bag of programs\. However, the program population may carry substantial structural information\. Effective program evolution depends on the existence of local structural similarity in the program space\. rather than merely high\-level fitness comparison\. Without such locality, small mutations induce unpredictable changes in program behavior, turning the search process into a largely stochastic exploration\. In contrast, when programs exhibit meaningful local similarity, neighboring programs tend to induce more consistent and gradual changes in functionality, allowing iterative mutations to progressively refine performance\.\(3\) Static mutation operator\.In AlphaEvolve, the LLM for mutation is fixed throughout the entire search process, which constrains the space of program edits to a static and pre\-trained distribution\. As a result, the system continues to sample from the same set of structural transformations regardless of the accumulated search experience\.

#### Our contribution\.

We presentGraph\-AugmetedEvolution, a framework that simultaneously resolves all three limitations within the OpenEvolve backbone\. Arelational GNNparses each program into a typed computation graph and is trained online to predict fitness \(§[3\.3](https://arxiv.org/html/2607.10127#S3.SS3)\)\. ADiscrete SAC meta\-controllertakes each program’s GNN embedding and learns which structural edit to apply next, balancing task\-score improvement against population novelty and complexity via automatic entropy regularization \(§[3\.4](https://arxiv.org/html/2607.10127#S3.SS4)\)\. AGRPO fine\-tuning loopuses group\-normalized evaluation rewards and the parent’s GNN embedding as a variance\-reduction baseline to update the LLM mutation operator online, progressively sharpening mutation quality throughout search \(§[3\.5](https://arxiv.org/html/2607.10127#S3.SS5)\)\. We evaluate GAE on an important scientific task—symbolic regression for nonlinear oscillator systems\. We utilize the benchmark datasets from LLM\-SR\(shojaee2024llm\)for our experiments\. In LLM\-SR \(§[4](https://arxiv.org/html/2607.10127#S4)\) experiment, GAE discovers novel programs that consistently match or surpass prior baselines on the benchmarks, achieving competitive or state\-of\-the\-art performance\. The remainder of the paper is organized as follows: §[2](https://arxiv.org/html/2607.10127#S2)surveys related work, §[3](https://arxiv.org/html/2607.10127#S3)describes the GAE framework in full, §[4](https://arxiv.org/html/2607.10127#S4)presents experimental results, §[5](https://arxiv.org/html/2607.10127#S5)discusses limitations and future directions\.

## 2Related Work

Prior work addresses at most one of the above\-mentioned bottlenecks in isolation\. EvoTune\(surina2025algorithm\)finetune the LLM via DPO and GRPO respectively, yet leave parent sampling uninformed\. Surrogate\-based optimization\(white2021bananas;falkner2018bohb\)improves evaluation efficiency but typically operates over fixed, hand\-encoded representations and does not support open\-ended code evolution\. Quality\-Diversity RL methods\(nilsson2021policy;faldor2023map;batra2023proximal\)replace random selection with policy\-gradient updates, but require dense per\-step rewards unavailable when fitness arrives only at program termination\. No existing method closes all three gaps simultaneously\.

GAE builds on LLM\-driven evolutionary search but departs from its static operator assumption\. FunSearch\(romera2024mathematical\), AlphaEvolve\(novikov2025alphaevolve\), and OpenEvolve\(openevolve2025\)all freeze model weights throughout search, relying on prompt diversity alone for exploration\. EvoPrompting\(chen2023evopromptinglanguagemodelscodelevel\)extends this idea through few\-shot prompting, yet still fixes the selection policy\. PACEvolve\(yan2026pacevolve\)addresses scaffold\-level failure modes, such as context pollution, mode collapse, and weak collaboration, through hand\-crafted rules, but leaves the mutation model and selection policy frozen\. ThetaEvolve\(wang2025thetaevolve\)introduces test\-time RL fine\-tuning to improve mutation quality, yet ties the weight update to a single task\-specific objective rather than a general reusable selector\. ShinkaEvolve\(lange2025shinkaevolve\)improves sample efficiency via a bandit\-based LLM ensemble and novelty\-based rejection filtering, yet the selection policy remains static and prompt\-driven with no learned component\. Our work adds a learned selector and fine\-tuning loop on top of this backbone, so both the*which\-to\-mutate*and*how\-to\-mutate*decisions improve over time\.

Maintaining population diversity while improving selection is precisely the goal of quality\-diversity RL, yet existing methods expose a fundamental mismatch with program evolution\. MAP\-Elites\(mouret2015illuminating\)archives diverse elites but samples parents uniformly\. Policy gradient extensions PGA\-MAP\-Elites\(nilsson2021policy\), DCG\-MAP\-Elites\(faldor2023map\), DCRL\-MAP\-Elites\(faldor2025synergizing\), and PPGA\(batra2023proximal\)improve this but all require dense per\-step rewards absent in whole\-program evaluation\. AURORA\(grillotti2022unsupervised\)relaxes hand\-specified feature axes by learning behavioural descriptors unsupervisedly\. GAE adopts this spirit for archive indexing while bridging the dense\-reward gap through graph\-structured signals propagated across the population\.

## 3Methodology

### 3\.1Background: MAP\-Elites and OpenEvolve

#### MAP\-Elites\.

Let𝒳\\mathcal\{X\}be a program space andℬ\\mathcal\{B\}akk\-dimensional behaviour descriptor space binned into cells𝒞\\mathcal\{C\}\. MAP\-Elites maintains an archive𝒜:𝒞→𝒳∪\{∅\}\\mathcal\{A\}:\\mathcal\{C\}\\to\\mathcal\{X\}\\cup\\\{\\emptyset\\\}mapping each cell to its current elite program\. The QD\-score aggregates quality and coverage:

QD​\(𝒜\)=∑c∈𝒞f​\(𝒜​\[c\]\)⋅𝟏​\[𝒜​\[c\]≠∅\]\.\\mathrm\{QD\}\(\\mathcal\{A\}\)=\\sum\_\{c\\in\\mathcal\{C\}\}f\(\\mathcal\{A\}\[c\]\)\\cdot\\mathbf\{1\}\[\\mathcal\{A\}\[c\]\\neq\\emptyset\]\.The deterministic update rule inserts a new programβ\\betainto𝒜\\mathcal\{A\}only if it strictly improves its cell’s incumbent, guaranteeingQD\\mathrm\{QD\}is non\-decreasing\.

#### OpenEvolve\.

OpenEvolve maintains a program archive𝒜\\mathcal\{A\}as a quality–diversity population, using a MAP\-Elites scheme in which programs are mapped onto a multi\-dimensional feature grid and the best performer is retained per cell\(openevolve\)\. At each generation, it samples a parent from𝒜\\mathcal\{A\}, assembles an LLM prompt from the parent code and top\-kkelites \(drawn separately from the parent, so inspiration programs differ from those shown to the LLM\), and queries the LLM for a mutated child via SEARCH/REPLACE diffs or full rewrites\. The child is scored by a user\-supplied evaluator under a cascade pattern of multi\-stage validation and admitted into𝒜\\mathcal\{A\}only if it dominates its cell’s occupant\. For diversity at scale,𝒜\\mathcal\{A\}is partitioned into multiple islands that evolve independently with periodic migration to prevent premature convergence, all orchestrated in parallel by a central controller\.

### 3\.2Graph\-Augmented Evolution

Graph\-Augmented Evolution extends the OpenEvolve loop with three co\-adaptive components, all underpinned by a shared relational GNN encoder \(§[3\.3](https://arxiv.org/html/2607.10127#S3.SS3)\) that parses each program’s abstract syntax tree into a typed computation graph and produces a structural embeddingz∈ℝdz\\in\\mathbb\{R\}^\{d\}\. Parents are drawn uniformly from a fixed\-size elite archive; the shared embedding then conditions a discrete Soft Actor–Critic meta\-controller \(§[3\.4](https://arxiv.org/html/2607.10127#S3.SS4)\) that maps the parent statezzto a mutation\-type hint, steering the LLM toward a structurally promising class of edit rather than an undirected one\. Its policy is trained on a three\-component reward \(task score, embedding\-space novelty, and a complexity penalty\) computed on the best child of each generation\. Whenever the LLM proposes a group ofGGchildren, online GRPO fine\-tuning \(§[3\.5](https://arxiv.org/html/2607.10127#S3.SS5)\) converts their group\-relative advantagesA^i∝Rβi−Rα\\hat\{A\}\_\{i\}\\propto R\_\{\\beta\_\{i\}\}\-R\_\{\\alpha\}\(child reward minus parent reward\) into an in\-place PPO\-clip update of the LLM weights, adapting the mutation policy without interrupting evolution while leaving the archive and the GNN unchanged\. The GNN encoder itself is updated separately everyKKgenerations via an MSE surrogate loss that regresses its embeddings onto observed rewards, after which all archive embeddings are refreshed\. The three components do not share gradients and operate at different timescales, keeping each module independently replaceable\. Figure[1](https://arxiv.org/html/2607.10127#S3.F1)illustrates the complete architecture, and Algorithm[1](https://arxiv.org/html/2607.10127#alg1)summarizes one generation of the GAE loop\.

![Refer to caption](https://arxiv.org/html/2607.10127v1/figures/GAEvolve_framework.png)Figure 1:GAE framework overview\.A relational GNN encoder parses each program’s abstract syntax tree into a typed computation graph and produces a structural embeddingzz, fine\-tuned online against observed rewards to track the evolving program distribution\. Parents are sampled from a fixed\-size elite archive; the parent embeddingzαz\_\{\\alpha\}then serves as the state of a Discrete Soft Actor–Critic meta\-controller \(policyπϕ\\pi\_\{\\phi\}with twin criticsQθ1,Qθ2Q\_\{\\theta\_\{1\}\},Q\_\{\\theta\_\{2\}\}\), which samples a mutation\-type hinta∼πϕ\(⋅∣zα\)a\\sim\\pi\_\{\\phi\}\(\\cdot\\mid z\_\{\\alpha\}\)that steers the mutation\. Each parent is mutated by a local LLM generatingG=8G=8candidates, which are scored by a numerical evaluator\. The SAC policy is trained on the best child of each generation with rewardr=clip​\(R​\(β\)\)\+λ1⋅novelty−λ2⋅complexityr=\\mathrm\{clip\}\(R\(\\beta\)\)\+\\lambda\_\{1\}\\cdot\\text\{novelty\}\-\\lambda\_\{2\}\\cdot\\text\{complexity\}, while the children’s group\-relative rewards adapt the LLM generation policy via GRPO\.Algorithm 1GAE — One Generation1:archive

𝒫\\mathcal\{P\}\(fixed size

NN\); GNN encoder

Φ:code→ℝd\\Phi:\\text\{code\}\\\!\\to\\\!\\mathbb\{R\}^\{d\}\(

d=128d\{=\}128\); SAC policy

πϕ\\pi\_\{\\phi\}over action space

𝒜\\mathcal\{A\}with twin critics

Qθ1,Qθ2Q\_\{\\theta\_\{1\}\},Q\_\{\\theta\_\{2\}\}; LLM operator

πθLLM\\pi\_\{\\theta\_\{\\mathrm\{LLM\}\}\}with reference

πref\\pi\_\{\\mathrm\{ref\}\}; replay buffer

ℬ\\mathcal\{B\}; GNN window

ℒ\\mathcal\{L\}; group size

GG; GNN interval

KK; warm\-up

WW; generation index

ss\.

2:// Parent selection and mutation hint

3:

α∼Uniform​\(𝒫\);zα←Φ​\(α\)\\alpha\\sim\\mathrm\{Uniform\}\(\\mathcal\{P\}\);\\quad z\_\{\\alpha\}\\leftarrow\\Phi\(\\alpha\)
4:

a←Uniform​\(𝒜\)a\\leftarrow\\mathrm\{Uniform\}\(\\mathcal\{A\}\)if

s<Ws<W, else

a∼πϕ\(⋅∣zα\)a\\sim\\pi\_\{\\phi\}\(\\cdot\\mid z\_\{\\alpha\}\)
5:// Generate and evaluate a group ofGGchildren

6:

\{β1,…,βG\}∼πθLLM\(⋅∣α,a\)\\\{\\beta\_\{1\},\\dots,\\beta\_\{G\}\\\}\\sim\\pi\_\{\\theta\_\{\\mathrm\{LLM\}\}\}\(\\cdot\\mid\\alpha,a\)
7:

Rβi←Eval​\(βi\)=−log10⁡NMSE​\(βi\),i=1,…,GR\_\{\\beta\_\{i\}\}\\leftarrow\\textsc\{Eval\}\(\\beta\_\{i\}\)=\-\\log\_\{10\}\\mathrm\{NMSE\}\(\\beta\_\{i\}\),\\;\\;i=1,\\dots,G
8:// GRPO update of the LLM operator \(§3\.5\)

9:

A^i←\(Rigrp−Rgrp¯\)/\(std​\(Rgrp\)\+ε\)\\hat\{A\}\_\{i\}\\leftarrow\\big\(R^\{\\mathrm\{grp\}\}\_\{i\}\-\\overline\{R^\{\\mathrm\{grp\}\}\}\\big\)/\(\\mathrm\{std\}\(R^\{\\mathrm\{grp\}\}\)\+\\varepsilon\),

Rigrp=Rβi−RαR^\{\\mathrm\{grp\}\}\_\{i\}=R\_\{\\beta\_\{i\}\}\-R\_\{\\alpha\}
10:Update

θLLM\\theta\_\{\\mathrm\{LLM\}\}via the PPO\-clip surrogate

ℒGRPO\\mathcal\{L\}\_\{\\mathrm\{GRPO\}\}\(Eq\. 8\)

11:// SAC meta\-controller update, best child only \(§3\.4\)

12:

β⋆←arg⁡maxi⁡Rβi;z′←Φ​\(β⋆\)\\beta^\{\\star\}\\leftarrow\\arg\\max\_\{i\}R\_\{\\beta\_\{i\}\};\\quad z^\{\\prime\}\\leftarrow\\Phi\(\\beta^\{\\star\}\)
13:

r←clip​\(Rβ⋆,rmin,rmax\)\+λ1​rnov​\(z′,𝒫\)−λ2​rcx​\(β⋆\)r\\leftarrow\\mathrm\{clip\}\(R\_\{\\beta^\{\\star\}\},r\_\{\\min\},r\_\{\\max\}\)\+\\lambda\_\{1\}r\_\{\\mathrm\{nov\}\}\(z^\{\\prime\},\\mathcal\{P\}\)\-\\lambda\_\{2\}r\_\{\\mathrm\{cx\}\}\(\\beta^\{\\star\}\)\(Eq\. 6\)

14:

ℬ\.push​\(zα,a,r,z′\)\\mathcal\{B\}\.\\mathrm\{push\}\(z\_\{\\alpha\},a,r,z^\{\\prime\}\); update

πϕ,Qθ1,Qθ2\\pi\_\{\\phi\},Q\_\{\\theta\_\{1\}\},Q\_\{\\theta\_\{2\}\}\(Eq\. 7\) if

\|ℬ\|≥32\|\\mathcal\{B\}\|\\geq 32
15:// Archive and GNN\-buffer update \(all children\)

16:for

i=1,…,Gi=1,\\dots,Gdo

17:

𝒫\.Insert​\(βi,Rβi,Φ​\(βi\)\)\\mathcal\{P\}\.\\mathrm\{Insert\}\(\\beta\_\{i\},R\_\{\\beta\_\{i\}\},\\Phi\(\\beta\_\{i\}\)\), keep top\-

NNby score

18:append

\(graph​\(βi\),Rβi\)\(\\mathrm\{graph\}\(\\beta\_\{i\}\),R\_\{\\beta\_\{i\}\}\)to

ℒ\\mathcal\{L\}if

Rβi\>−10R\_\{\\beta\_\{i\}\}\>\-10
19:endfor

20:// Periodic GNN fine\-tuning \(§3\.3\)

21:if

smodK=0s\\bmod K=0and

\|ℒ\|≥10\|\\mathcal\{L\}\|\\geq 10: fine\-tune

Φ\\Phi\(Eq\. 5\); refresh

zp​∀p∈𝒫z\_\{p\}\\;\\forall p\\in\\mathcal\{P\}

### 3\.3Program Graphs and GNN Encoder

#### Graph representation\.

We parse each candidate program into a typed computation graph𝒢=\(𝒱,ℰ\)\\mathcal\{G\}=\(\\mathcal\{V\},\\mathcal\{E\}\)using static Abstract Syntax Tree \(AST\) analysis\. Any parse failure produces a singleunknownnode, guaranteeing a valid graph is always returned\. Specifically for symbolic regression, nodes are mathematical operations and edges encode data\-flow\.

#### Relational GNN encoder\.

Each node carries anFF\-dimensional one\-hot encoding its node type\. We embed each program graph intoz∈ℝdz\\in\\mathbb\{R\}^\{d\}via a relational graph neural network\. Suppose we denote the archive size asNN, theinput projection is

h\(0\)=GELU​\(X​Win\+bin\),X∈ℝN×F,Win∈ℝF×h\.h^\{\(0\)\}=\\mathrm\{GELU\}\(XW\_\{\\mathrm\{in\}\}\+b\_\{\\mathrm\{in\}\}\),\\quad X\\in\\mathbb\{R\}^\{N\\times F\},\\;W\_\{\\mathrm\{in\}\}\\in\\mathbb\{R\}^\{F\\times h\}\.\(1\)with hidden widthhh\(the output dimension is produced only by the final pooling projection, where2​h=d2h=d\)\.

#### Relational message passing\.

We applyLLlayers of edge\-type\-conditioned convolution\. For each relationr∈ℰr\\in\\mathcal\{E\}, a dedicated weight matrixWr∈ℝh×hW\_\{r\}\\in\\mathbb\{R\}^\{h\\times h\}transforms the features of the source node\. Messages are mean\-aggregated at each destination and combined with a self\-loop:

hi\(ℓ\+1\)=LayerNorm\(GELU\(Wselfhi\(ℓ\)\+∑r=1\|ℰ\|1\|𝒩r​\(i\)\|∑j∈𝒩r​\(i\)Wrhj\(ℓ\)\)\),h\_\{i\}^\{\(\\ell\+1\)\}=\\mathrm\{LayerNorm\}\\\!\\left\(\\mathrm\{GELU\}\\\!\\left\(W\_\{\\mathrm\{self\}\}\\,h\_\{i\}^\{\(\\ell\)\}\\right\.\\right\.\\\\ \\left\.\\left\.\+\\;\\sum\_\{r=1\}^\{\|\\mathcal\{E\}\|\}\\frac\{1\}\{\|\\mathcal\{N\}\_\{r\}\(i\)\|\}\\sum\_\{j\\in\\mathcal\{N\}\_\{r\}\(i\)\}W\_\{r\}\\,h\_\{j\}^\{\(\\ell\)\}\\right\)\\right\),\(2\)
where𝒩r​\(i\)=\{j:\(j,i\)∈ℰr\}\\mathcal\{N\}\_\{r\}\(i\)=\\\{\\,j:\(j,i\)\\in\\mathcal\{E\}\_\{r\}\\,\\\}denotes the set of neighbors of nodeiiconnected by an edge of relation typerr, and\|𝒩r​\(i\)\|\|\\mathcal\{N\}\_\{r\}\(i\)\|is its cardinality\.

#### Global pooling\.

Mean and max pooling are concatenated and projected to the output representationzzas:

z\\displaystyle z=GELU​\(LayerNorm​\(\[h¯;h^\]​Wout\)\),\\displaystyle=\\mathrm\{GELU\}\\\!\\left\(\\mathrm\{LayerNorm\}\\\!\\left\(\[\\bar\{h\};\\,\\hat\{h\}\]\\,W\_\{\\mathrm\{out\}\}\\right\)\\right\),\(3\)h¯\\displaystyle\\bar\{h\}=1N​∑ihi\(L\),h^=maxi⁡hi\(L\)\.\\displaystyle=\\tfrac\{1\}\{N\}\\textstyle\\sum\_\{i\}h\_\{i\}^\{\(L\)\},\\;\\hat\{h\}=\\max\_\{i\}h\_\{i\}^\{\(L\)\}\.\(4\)

#### Pretraining and online fine\-tuning\.

Before evolution begins, the encoder is pretrained on synthetic programs with structural proxy scores \(feature coverage, operator diversity, complexity\), using a linear score\-prediction head trained jointly with the encoder\. During evolution, every 3 steps \(once≥10\\geq 10samples exist\), the encoder is fine\-tuned for 300 gradient steps on a sliding window of the 300 most recent\(zα,R​\(α\)\)\(z\_\{\\alpha\},R\(\\alpha\)\)pairs, wherezα∈ℝdz\_\{\\alpha\}\\in\\mathbb\{R\}^\{d\}is the GNN embedding of programα\\alphaandR​\(α\)∈ℝR\(\\alpha\)\\in\\mathbb\{R\}is its evaluator score\. Targets are standardized to zero mean and unit variance\. All population embeddings are then refreshed\. The training objective is:

ℒGNN=1\|ℒ\|​∑α∈ℒ\(μα−R​\(α\)\)2,\\mathcal\{L\}\_\{\\mathrm\{GNN\}\}=\\frac\{1\}\{\|\\mathcal\{L\}\|\}\\sum\_\{\\alpha\\in\\mathcal\{L\}\}\\\!\\left\(\\mu\_\{\\alpha\}\-R\(\\alpha\)\\right\)^\{2\},\(5\)whereμα=MLP​\(zα\)\\mu\_\{\\alpha\}=\\mathrm\{MLP\}\(z\_\{\\alpha\}\)is a scalar score prediction from a temporary head trained jointly with the encoder\.

### 3\.4RL\-Guided Mutation Selection

At each evolution step a parentα\\alphais drawn*uniformly at random*from the current island; the reinforcement\-learning agent does not select parents\. Instead, a Discrete Soft Actor\-Critic \(SAC\) agent chooses which*mutation*to apply toα\\alpha, conditioned on the parent’s graph embedding\.

#### State and action space\.

The state is the GNN embeddingzα∈ℝdz\_\{\\alpha\}\\in\\mathbb\{R\}^\{d\}of the parent program\. The action space is flat and fixed, with\|𝒜\|\|\\mathcal\{A\}\|discrete actions, each encoding a tuple \(mutation type, target node id, argument index\)\.

#### Networks\.

The policyπϕ​\(a∣z\)\\pi\_\{\\phi\}\(a\\mid z\)and the twin criticsQθ1,Qθ2:ℝd→ℝ\|𝒜\|Q\_\{\\theta\_\{1\}\},Q\_\{\\theta\_\{2\}\}:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}^\{\|\\mathcal\{A\}\|\}are three\-layer MLPs with LayerNorm and GELU activations that score every action from the state:

z→256LayerNorm→GELU→256GELU→ℝ\|𝒜\|\.z\\xrightarrow\{\\;256\\;\}\\mathrm\{LayerNorm\}\\to\\mathrm\{GELU\}\\xrightarrow\{\\;256\\;\}\\mathrm\{GELU\}\\to\\mathbb\{R\}^\{\|\\mathcal\{A\}\|\}\.Each critic is paired with a target network updated by Polyak averaging \(coefficientτ\\tau\)\.

#### Action selection\.

During a warm\-up ofWWsteps, actions are sampled uniformly at random to seed the replay buffer with diverse transitions\. Thereafter, the agent samples an action stochastically from the softmax policyπϕ​\(a∣zα\)\\pi\_\{\\phi\}\(a\\mid z\_\{\\alpha\}\)and no policy update is performed until the buffer holds at least3232transitions\.

#### Reward\.

The mutation producing the best childβ\\betareceives a three\-component reward:

r\\displaystyle r=clip​\(R​\(β\),rmin,rmax\)⏟rtask\+λ1​\(1−minp∈𝒫⁡cos⁡\(zβ,zp\)\)⏟rnovelty\\displaystyle=\\underbrace\{\\mathrm\{clip\}\\bigl\(R\(\\beta\),\\,r\_\{\\min\},\\,r\_\{\\max\}\\bigr\)\}\_\{r\_\{\\text\{task\}\}\}\+\\underbrace\{\\lambda\_\{1\}\\bigl\(1\-\\min\_\{p\\in\\mathcal\{P\}\}\\cos\(z\_\{\\beta\},z\_\{p\}\)\\bigr\)\}\_\{r\_\{\\text\{novelty\}\}\}−λ2​max⁡\(0,N​\(β\)−2020\)⏟rcomplexity,\\displaystyle\\quad\-\\underbrace\{\\lambda\_\{2\}\\,\\max\\\!\\Bigl\(0,\\,\\tfrac\{N\(\\beta\)\-20\}\{20\}\\Bigr\)\}\_\{r\_\{\\text\{complexity\}\}\},\(6\)whereR​\(β\)R\(\\beta\)is the evaluator score of the best child,𝒫\\mathcal\{P\}is the current population,N​\(β\)N\(\\beta\)is the AST node count ofβ\\beta’s return expression\. The task term uses the child’s absolute \(clipped\) score rather than a parent\-relative improvement\.

#### Online SAC training\.

Each transition\(zα,a,r,zβ\)\(z\_\{\\alpha\},a,r,z\_\{\\beta\}\)is stored in a prioritized replay buffer of capacity\. Every generation we takeKKgradient steps \(equal to the number of children per generation\) on priority\-sampled mini\-batches, updating the twin critics, the policy, and the entropy temperatureα\\alpha\(automatically tuned toward a target entropy of−12​log⁡\|𝒜\|\-\\tfrac\{1\}\{2\}\\log\|\\mathcal\{A\}\|\)\. The critics minimize an importance\-weighted TD loss

ℒ​\(θi\)\\displaystyle\\mathcal\{L\}\(\\theta\_\{i\}\)=𝔼\(z,a,r,z′\)∼ℬ​\[w​\(Qθi​\(z,a\)−y\)2\],\\displaystyle=\\mathbb\{E\}\_\{\(z,a,r,z^\{\\prime\}\)\\sim\\mathcal\{B\}\}\\\!\\left\[w\\,\\bigl\(Q\_\{\\theta\_\{i\}\}\(z,a\)\-y\\bigr\)^\{2\}\\right\],y\\displaystyle y=r\+γ​V​\(z′\),\\displaystyle=r\+\\gamma\\,V\(z^\{\\prime\}\),with bootstrapped value

V​\(z′\)=∑aπϕ​\(a∣z′\)​\(mini⁡Qθ¯i​\(z′,a\)−α​log⁡πϕ​\(a∣z′\)\)V\(z^\{\\prime\}\)=\\sum\_\{a\}\\pi\_\{\\phi\}\(a\\mid z^\{\\prime\}\)\\bigl\(\\min\_\{i\}Q\_\{\\bar\{\\theta\}\_\{i\}\}\(z^\{\\prime\},a\)\-\\alpha\\log\\pi\_\{\\phi\}\(a\\mid z^\{\\prime\}\)\\bigr\)\(7\)and importance weightww\. The policy maximizes the standard entropy\-regularized objective

J\(π\)=𝔼z∼ℬ\[∑aπϕ\(a∣z\)\(\\displaystyle J\(\\pi\)=\\mathbb\{E\}\_\{z\\sim\\mathcal\{B\}\}\\\!\\biggl\[\\sum\_\{a\}\\pi\_\{\\phi\}\(a\\mid z\)\\bigl\(mini⁡Qθi​\(z,a\)\\displaystyle\\min\_\{i\}Q\_\{\\theta\_\{i\}\}\(z,a\)−αlogπϕ\(a∣z\)\)\]\.\\displaystyle\-\\alpha\\log\\pi\_\{\\phi\}\(a\\mid z\)\\bigr\)\\biggr\]\.

### 3\.5GRPO LLM Fine\-Tuning

On every evolution, we construct a group ofGGchildren from the sampled parentα\\alpha,βg∼πθLLM\(⋅∣α,a\),g=1,…,G\\beta\_\{g\}\\sim\\pi\_\{\\theta\_\{\\mathrm\{LLM\}\}\}\(\\cdot\\mid\\alpha,a\),\\ g=1,\\ldots,G, withaadenoting the action selected by the SAC meta\-controllor\. A GRPO update is performed immediately after the group is evaluated\. Every valid child is truly evaluated, yielding rewardsRg\{R\_\{g\}\}and malformed completions are discarded before evaluation\.

The group reward is baseline\-corrected by the parent’s own \(true\) evaluated scoreRαR\_\{\\alpha\}:Rggrp=Rg−RαR^\{\\mathrm\{grp\}\}\_\{g\}=R\_\{g\}\-R\_\{\\alpha\}, i\.e\. each child’s reward is measured as its improvement over the parent\. Ifstd​\(Rgrp\)<10−8\\mathrm\{std\}\(R^\{\\mathrm\{grp\}\}\)<10^\{\-8\}\(all rewards identical\), the update is skipped\. The GRPO objective fine\-tunesθLLM\\theta\_\{\\mathrm\{LLM\}\}with a PPO\-clip surrogate and KL penaltyβkl=0\.1\\beta\_\{\\mathrm\{kl\}\}=0\.1to the reference model:

ℒGRPO​\(θLLM\)=\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{GRPO\}\}\(\\theta\_\{\\mathrm\{LLM\}\}\)=\{\}−1G​∑gmin⁡\(ρg​A^g,clip​\(ρg,1±ϵclip\)​A^g\)\\displaystyle\-\\frac\{1\}\{G\}\\sum\_\{g\}\\min\\\!\\bigl\(\\rho\_\{g\}\\hat\{A\}\_\{g\},\\;\\mathrm\{clip\}\(\\rho\_\{g\},1\\pm\\epsilon\_\{\\mathrm\{clip\}\}\)\\,\\hat\{A\}\_\{g\}\\bigr\)\(8\)\+βkl​\[DKL​\(πθLLM∥πθLLMref\)\],\\displaystyle\+\\beta\_\{\\mathrm\{kl\}\}\\,\\bigl\[D\_\{\\mathrm\{KL\}\}\(\\pi\_\{\\theta\_\{\\mathrm\{LLM\}\}\}\\,\\\|\\,\\pi\_\{\\theta\_\{\\mathrm\{LLM\}\}\}^\{\\mathrm\{ref\}\}\)\\bigr\],whereρg=exp⁡\(log⁡πθLLM​\(βg\|α\)−log⁡πθLLMref​\(βg\|α\)\)\\rho\_\{g\}=\\exp\\,\\\!\\bigl\(\\log\\pi\_\{\\theta\_\{\\mathrm\{LLM\}\}\}\(\\beta\_\{g\}\|\\alpha\)\-\\log\\pi\_\{\\theta\_\{\\mathrm\{LLM\}\}\}^\{\\mathrm\{ref\}\}\(\\beta\_\{g\}\|\\alpha\)\\bigr\),ϵclip\\epsilon\_\{\\mathrm\{clip\}\}is PPO clipping range, and the group\-normalized advantage of childgg

A^g=Rggrp−R¯grpstd​\(Rgrp\)\+ϵ\.\\hat\{A\}\_\{g\}=\\frac\{R^\{\\mathrm\{grp\}\}\_\{g\}\-\\bar\{R\}^\{\\mathrm\{grp\}\}\}\{\\mathrm\{std\}\(R^\{\\mathrm\{grp\}\}\)\+\\epsilon\}\.\(9\)The KL term is therefore evaluated for monitoring only and contributes no gradient\. The update ofθLLM\\theta\_\{\\mathrm\{LLM\}\}comes solely from the PPO\-clip surrogate, with clipping to the trust region\[1−ϵclip,1\+ϵclip\]\[1\-\\epsilon\_\{\\mathrm\{clip\}\},\\,1\+\\epsilon\_\{\\mathrm\{clip\}\}\]providing the effective regularization andβkl\\beta\_\{\\mathrm\{kl\}\}reported for diagnostics\.

Completions with reward below−10\-10\(evaluator crash sentinels\) are dropped before computing advantages\. After each group, all evaluated children are inserted into the fixed\-size elite population \(a top\-KKlist sorted by score, with the best child inserted first for elitism\), and the population is truncated to its maximum size\.

## 4Experiments

To evaluate GAEvolve’s capability of scientific discovery, we test its performance on the Nonlinear Oscillators task from LLM\-SR\. Nonlinear damped oscillators are described by differential equations that capture the complex interaction among an oscillator’s position, velocity, and the acting forces\. The goal is to discover a closed\-form expressionf^​\(x,t,v;θ\)\\hat\{f\}\(x,t,v;\\theta\)that predicts the accelerationx¨=d​v/d​t\\ddot\{x\}=dv/dtof a nonlinear harmonic oscillator, wherexxis position,ttis time,vvis velocity\. We evaluate GAEvolve on discovering the governing equation of a synthetically generated variant by minimizing the Normalized Mean Squared Error \(NMSE\) between predicted and observed trajectories\. The actual evolved program expresses a candidate symbolic equation as executable code\. The initial program is a naive linear modelv˙^=p0​x\+p1​t\+p2​v\+p3\\hat\{\\dot\{v\}\}=p\_\{0\}x\+p\_\{1\}t\+p\_\{2\}v\+p\_\{3\}, parsed into an symbolic regression expression tree via the graph builder, producing a 4\-node graph withadd,mul,input\_varandparamnodes\.

#### Dataset\.

We use thephys\_osc\(nonlinear harmonic oscillator\) domain of the LSR\-Synth subset of LLM\-SRBench\(shojaee2025llm\), predicting the accelerationv˙\\dot\{v\}from\(x,t,v\)\(x,t,v\)\. Each task supplies three disjoint sets: a training set, an in\-distribution \(ID\) test set, and an out\-of\-distribution \(OOD\) test set sampled from a different input time region\. For every candidate equation, the numeric constants are optimized by BFGS to minimize the mean squared error \(MSE\) on the training set\. The candidate is then scored by−log10⁡\(NMSE\)\-\\log\_\{10\}\(\\mathrm\{NMSE\}\)on the held\-out ID test set, which serves as the evolutionary fitness\. We report the median NMSE on the ID and OOD test sets, and the latter quantifying extrapolation\.

#### Experimental Setup\.

GAE maintains a fixed\-size elite population of 50 programs sorted by fitness\. Parents are sampled uniformly from this population, and a discrete Soft Actor\-Critic meta\-controller \(with 32 random warm\-up steps\) selects the mutation\-type hint passed to the LLM\. Each generation samples eight LLM\-generated candidates in parallel\. Before committing any LLM\-generated candidate to the expensive numerical optimization pipeline, we apply a lightweight four\-stage validation chain that filters out malformed programs early\. We \(i\) extract valid Python code from the raw LLM output, \(ii\) verify syntactic correctness, \(iii\) sandbox the candidate by stripping unsafe statements, and \(iv\) execute a fast smoke test on randomly generated dummy inputs to catch runtime errors\. Candidates failing any stage are discarded immediately, preventing invalid programs from entering the evaluator or corrupting the reinforcement learning replay buffer\.

For the OpenEvolve baseline, the MAP\-Elites archive is two\-dimensional \(expression tree depth, number of distinct AST node types as structural diversity\), with a population of 200 programs across 4 isolated islands and periodic migration every 15 generations\. Parent sampling uses fixed elite/exploration/exploitation, each at the ratio of 0\.3/0\.4/0\.3\. For the PACEvolve baseline, We use the official implementation with a single island\. The program database maintains a score\-sorted queue of at most 100 programs, and each parent is drawn by tournament selection of size 2 from the top\-4 programs\. Each iteration generates one candidate with the frozen LLM backbone with temperature at 0\.7 and top\-p at 0\.95\)\. Before code generation, a scratch\-pad step proposes a hypothesis, a repository capped at 5 ideas with periodic merging, and backtracking which restores an earlier repository state sampled from a power\-law distribution\.

We use Qwen3\.5\-35B\-A3B\(qwen3\.5\_2026\)as the LLM backbone of GAE\. We fix the LLM weights for PACEvolve\-Single and OpenEvolve and update it in full weights in GAE during the GRPO reinforcement learning period\. The GRPO has a group size ofG=8G=8, and uses a KL regularization coefficient ofβ=0\.1\\beta=0\.1together with a clipping range ofε=0\.2\\varepsilon=0\.2to stabilize policy updates\. Optimization is carried out using AdamW with a learning rate of1×10−61\\times 10^\{\-6\}\. Fine\-tuning is performed in an online manner, with updates triggered at every generation after candidate evaluation\.

MethodMean NMSEBest NMSENon\-LLM baselineuDSR\(landajuela2022unified\)−3\.95\-3\.95−4\.06\-4\.06LLM\-guided evolutionLLM\-SR−4\.06\-4\.06−4\.80\-4\.80OpenEvolve−5\.40\-5\.40−7\.11\-7\.11CodeEvolve−4\.97\-4\.97−7\.26\-7\.26ShinkaEvolve−5\.35\-5\.35−6\.35\-6\.35PACEvolve\-Sigle−5\.87\-5\.87−8\.23\-8\.23PACEvolve\-Multi−6\.11\-6\.11−8\.24\-8\.24Table 1:Benchmark results on the LLM\-SR Nonlinear Oscillators task\(yan2026pacevolve\), reported aslog10\\log\_\{10\}NMSE \(↓\\downarrow\)\. All LLM\-guided baselines run for 1000 iterations over 10 independent runs\. Both the mean and the bestlog10\\log\_\{10\}NMSE across runs are reported\.
MethodID NMSE↓\\downarrowOOD NMSE↓\\downarrowPACEvolve\-Single \(Gemini\)−5\.80\-5\.80−6\.10\-6\.10OpenEvolve \(Gemini\)−3\.14\-3\.14−3\.26\-3\.26OpenEvolve \(Qwen\)−3\.88\-3\.88−3\.94\-3\.94GAE \(ours\)−6\.87\\mathbf\{\-6\.87\}\(Mean\)−7\.24\\mathbf\{\-7\.24\}\(Best\)−7\.10\\mathbf\{\-7\.10\}\(Mean\)−7\.45\\mathbf\{\-7\.45\}\(Best\)Table 2:Comparison among PACEvolve, OpenEvolve \(both 1 run\) and our method GAE \(3 runs; mean and best reported\) on the Nonlinear Oscillators task under our evaluation setting \(150 iterations\), reported aslog10\\log\_\{10\}NMSE \(↓\\downarrow\) on in\-distribution \(ID\) and out\-of\-distribution \(OOD\) test sets\.

#### Main Results\.

First, we present results from PACEvolve\(yan2026pacevolve\)under their standard high\-budget evaluation protocol \(1000 iterations, 10 runs\), which compares against uDSR\(landajuela2022unified\), LLM\-SR, AlphaEvolve, ShinkaEvolve\(lange2025shinkaevolve\)\. Table[2](https://arxiv.org/html/2607.10127#S4.T2)reflects the performance of existing methods under their originally reported compute regime and provides a reference for overall competitiveness in the Non\-linear harmonic oscillator task\. Second, we conduct a controlled, low\-budget comparison under a fixed backbone setting, using either Qwen or Gemini 2\.5 Flash as the shared base model\. In this regime, we compare theOpenEvolvebaseline, the strong prior methodPaceEvolve, and our proposed methodGAE, all run for 150 iterations to ensure identical computational budget\. As shown in Table[2](https://arxiv.org/html/2607.10127#S4.T2),GAEachievesbestperformance compared to both baselines, with its ID NMSE at \-7\.24 and OOD NMSE at \-7\.45 across all methods\.

![Refer to caption](https://arxiv.org/html/2607.10127v1/figures/gae_metric_trajectory_final.png)Figure 2:Performance trajectory of the GAE method on the Non\-linear Harmonic Oscillator symbolic regression problem\. The curve shows the running\-bestlog10⁡\(NMSE\)\\log\_\{10\}\(\\text\{NMSE\}\)as a function of the number of evolution iterations\. Grey dots mark iterations at which a new best solution was discovered\. Lower values oflog10⁡\(NMSE\)\\log\_\{10\}\(\\text\{NMSE\}\)indicate better fit\.Figure[2](https://arxiv.org/html/2607.10127#S4.F2)shows the evolution trajectory of GAE over 150 iteration\. The method starts from an initial solution withlog10⁡\(NMSE\)=−1\.26\\log\_\{10\}\(\\text\{NMSE\}\)=\{\-1\.26\}and improves rapidly during the first 20 iterations, reachinglog10⁡\(NMSE\)≈−4\.64\\log\_\{10\}\(\\text\{NMSE\}\)\\approx\{\-4\.64\}by iteration 17\. Steady further gains accumulate through iterations 57–88, where three consecutive breakthroughs push the score tolog10⁡\(NMSE\)=−7\.19\\log\_\{10\}\(\\text\{NMSE\}\)=\{\-7\.19\}\. A final improvement at iteration 119 yields the best solution of the run, withlog10⁡\(NMSE\)=−7\.24\\log\_\{10\}\(\\text\{NMSE\}\)=\{\-7\.24\}\(NMSE≈5\.7×10−8\\approx 5\.7\\times 10^\{\-8\}\), after which the trajectory plateaus for the remaining 30 iterations\. Across all 150 iterations, 13 distinct improvements were recorded, demonstrating that GAE efficiently discovers a near\-exact symbolic form of the target equation through graph\-structured evolutionary search\.

Next, we demonstrate the best program discovered by each of all four pipelines \(Table[2](https://arxiv.org/html/2607.10127#S4.T2)\) and provide physical interpretations of each component in the resulting equation\.

#### GAE best program

f​\(x,t,v;𝜽\)=\\displaystyle f\(x,t,v;\\,\\boldsymbol\{\\theta\}\)=θ0​x\+θ1​x2​sgn⁡\(x\)\+θ2​x3\+θ3​x4​sgn⁡\(x\)\\displaystyle\\theta\_\{0\}\\,x\+\\theta\_\{1\}\\,x^\{2\}\\operatorname\{sgn\}\(x\)\+\\theta\_\{2\}\\,x^\{3\}\+\\theta\_\{3\}\\,x^\{4\}\\operatorname\{sgn\}\(x\)\(10\)\+\\displaystyle\+θ4​v\+θ5​v2​sgn⁡\(v\)\+θ6​v3\\displaystyle\\theta\_\{4\}\\,v\+\\theta\_\{5\}\\,v^\{2\}\\operatorname\{sgn\}\(v\)\+\\theta\_\{6\}\\,v^\{3\}\+\\displaystyle\+θ7​sin⁡\(t\)\+θ8​cos⁡\(t\)\\displaystyle\\theta\_\{7\}\\,\\sin\(t\)\+\\theta\_\{8\}\\,\\cos\(t\)\+\\displaystyle\+θ9​\(\|x\|​\|v\|\+x2​\|v\|\+\|t\|​\|v\|\+t2​v\+x​v2\)\\displaystyle\\theta\_\{9\}\\,\\bigl\(\|x\|\|v\|\+x^\{2\}\|v\|\+\|t\|\|v\|\+t^\{2\}v\+xv^\{2\}\\bigr\)
The termsθ0​x\\theta\_\{0\}x,θ1​x2​sgn⁡\(x\)\\theta\_\{1\}x^\{2\}\\operatorname\{sgn\}\(x\),θ2​x3\\theta\_\{2\}x^\{3\}, andθ3​x4​sgn⁡\(x\)\\theta\_\{3\}x^\{4\}\\operatorname\{sgn\}\(x\)correspond to the linear restoring force, quadratic and cubic nonlinearity, and quartic nonlinearity respectively\. The termsθ4​v\\theta\_\{4\}v,θ5​v2​sgn⁡\(v\)\\theta\_\{5\}v^\{2\}\\operatorname\{sgn\}\(v\), andθ6​v3\\theta\_\{6\}v^\{3\}represent velocity\-dependent damping\. External time\-dependent forcing is captured byθ7​sin⁡\(t\)\\theta\_\{7\}\\sin\(t\)andθ8​cos⁡\(t\)\\theta\_\{8\}\\cos\(t\)\. The five remaining terms\|x\|​\|v\|\|x\|\|v\|,x2​\|v\|x^\{2\}\|v\|,\|t\|​\|v\|\|t\|\|v\|,t2​vt^\{2\}v, andx​v2xv^\{2\}, all scaled by the shared parameterθ9\\theta\_\{9\}, are interaction terms\.

#### OpenEvolve best program \(Qwen\)

f​\(x,t,v\)=\\displaystyle f\(x,t,v\)=θ1​x\+θ2​x3\+θ3​x5\+θ4​v\+θ5​v3\\displaystyle\\;\\theta\_\{1\}x\+\\theta\_\{2\}x^\{3\}\+\\theta\_\{3\}x^\{5\}\+\\theta\_\{4\}v\+\\theta\_\{5\}v^\{3\}\+θ6​x​v\+θ7​x2​v\+θ8​x​v2\\displaystyle\+\\theta\_\{6\}xv\+\\theta\_\{7\}x^\{2\}v\+\\theta\_\{8\}xv^\{2\}\+θ9​sin⁡t\+θ10​sin⁡\(2​t\)\+θ11​cos⁡t\\displaystyle\+\\theta\_\{9\}\\sin t\+\\theta\_\{10\}\\sin\(2t\)\+\\theta\_\{11\}\\cos t\+θ12​x​sin⁡t\+θ13​v​sin⁡t\\displaystyle\+\\theta\_\{12\}x\\sin t\+\\theta\_\{13\}v\\sin t\+θ14​t\+θ15​t2\+θ16​t3\\displaystyle\+\\theta\_\{14\}t\+\\theta\_\{15\}t^\{2\}\+\\theta\_\{16\}t^\{3\}\+θ17​11\+x2\+θ18​11\+v2\\displaystyle\+\\theta\_\{17\}\\frac\{1\}\{1\+x^\{2\}\}\+\\theta\_\{18\}\\frac\{1\}\{1\+v^\{2\}\}\+θ19​e−x2\+θ20​e−v2\\displaystyle\+\\theta\_\{19\}e^\{\-x^\{2\}\}\+\\theta\_\{20\}e^\{\-v^\{2\}\}\+θ21​sinh⁡\(x\)\+θ22​cosh⁡\(x\)\\displaystyle\+\\theta\_\{21\}\\sinh\(x\)\+\\theta\_\{22\}\\cosh\(x\)\+θ23​x2​t\+θ24​v2​t\+θ25​x​v​t\.\\displaystyle\+\\theta\_\{23\}x^\{2\}t\+\\theta\_\{24\}v^\{2\}t\+\\theta\_\{25\}xvt\.
The first three terms \(θ0​x\\theta\_\{0\}x,θ1​x3\\theta\_\{1\}x^\{3\},θ2​x5\\theta\_\{2\}x^\{5\}\) represent a linear, cubic \(Duffing\-like\), and quintic restoring force\. The next two \(θ3​v\\theta\_\{3\}v,θ4​v3\\theta\_\{4\}v^\{3\}\) capture linear and cubic damping\. The mixed termsθ5​x​v\\theta\_\{5\}xv,θ6​x2​v\\theta\_\{6\}x^\{2\}v, andθ7​x​v2\\theta\_\{7\}xv^\{2\}are mixed linear damping with position, mixed cubic damping \(x2​vx^\{2\}v\), and mixed cubic damping \(x​v2xv^\{2\}\)\. Harmonic, higher harmonic, and cosine driving are handled byθ8​sin⁡\(t\)\\theta\_\{8\}\\sin\(t\),θ9​sin⁡\(2​t\)\\theta\_\{9\}\\sin\(2t\), andθ10​cos⁡\(t\)\\theta\_\{10\}\\cos\(t\), whileθ11​x​sin⁡\(t\)\\theta\_\{11\}x\\sin\(t\)andθ12​v​sin⁡\(t\)\\theta\_\{12\}v\\sin\(t\)introduce combined time\-position and time\-velocity interactions\. The termsθ13​t\\theta\_\{13\}t,θ14​t2\\theta\_\{14\}t^\{2\},θ15​t3\\theta\_\{15\}t^\{3\}are polynomial time terms\. The termsθ16/\(1\+x2\)\\theta\_\{16\}/\(1\+x^\{2\}\)andθ17/\(1\+v2\)\\theta\_\{17\}/\(1\+v^\{2\}\)are inverse terms for potential singularities\. The termsθ18​e−x2\\theta\_\{18\}e^\{\-x^\{2\}\}andθ19​e−v2\\theta\_\{19\}e^\{\-v^\{2\}\}model exponential decay or growth in position and velocity, andθ20​sinh⁡\(x\)\\theta\_\{20\}\\sinh\(x\)andθ21​cosh⁡\(x\)\\theta\_\{21\}\\cosh\(x\)provide smooth hyperbolic transitions\. Finally,θ22​x2​t\\theta\_\{22\}x^\{2\}t,θ23​v2​t\\theta\_\{23\}v^\{2\}t, andθ24​x​v​t\\theta\_\{24\}xvtare cross\-terms\.

#### OpenEvolve best program \(Gemini\)

f​\(x,t,v;𝜽\)=\\displaystyle f\(x,t,v;\\,\\boldsymbol\{\\theta\}\)=θ0​x\+θ1​sin⁡\(t\)\+θ2​v\+θ3​x2​v\\displaystyle\\theta\_\{0\}\\,x\+\\theta\_\{1\}\\,\\sin\(t\)\+\\theta\_\{2\}\\,v\+\\theta\_\{3\}\\,x^\{2\}v\(11\)\+\\displaystyle\+θ4​x2\+θ5​v​\|v\|\+θ6​x​v2\+θ7​x3\\displaystyle\\theta\_\{4\}\\,x^\{2\}\+\\theta\_\{5\}\\,v\|v\|\+\\theta\_\{6\}\\,xv^\{2\}\+\\theta\_\{7\}\\,x^\{3\}\+\\displaystyle\+θ8​cos⁡\(t\)\+θ9​v3\\displaystyle\\theta\_\{8\}\\,\\cos\(t\)\+\\theta\_\{9\}\\,v^\{3\}
The linear termθ0​x\\theta\_\{0\}xrepresents a linear restoring force, whileθ1​sin⁡\(t\)\\theta\_\{1\}\\sin\(t\)andθ8​cos⁡\(t\)\\theta\_\{8\}\\cos\(t\)together capture periodic driving or intrinsic oscillation via sine and cosine components\. The termθ2​v\\theta\_\{2\}vencodes linear damping\. The termθ3​x2​v\\theta\_\{3\}x^\{2\}vis a higher\-order interaction combining velocity\-dependent nonlinear stiffness with position\-dependent nonlinear damping\. The quadratic termθ4​x2\\theta\_\{4\}x^\{2\}introduces an asymmetric potential or nonlinear restoring force\. Rather thanv2v^\{2\}, the termθ5​v​\|v\|\\theta\_\{5\}v\|v\|is used to model quadratic drag, where the damping magnitude scales withv2v^\{2\}but the sign always opposes the direction of motion—a physically accurate form of turbulent or aerodynamic damping\. The termθ6​x​v2\\theta\_\{6\}xv^\{2\}is a higher\-order nonlinear interaction between position and velocity\. The cubic termθ7​x3\\theta\_\{7\}x^\{3\}introduces a Duffing\-oscillator\-like strong nonlinear restoring force, andθ9​v3\\theta\_\{9\}v^\{3\}adds cubic nonlinear damping\.

#### PACEvolve Best Program

f​\(x,t,v;𝜽\)=\\displaystyle f\(x,t,v;\\,\\boldsymbol\{\\theta\}\)=θ0​x\+θ1​v3\+θ2​v\+θ3​x​e−v2\+θ4​x3\\displaystyle\\theta\_\{0\}\\,x\+\\theta\_\{1\}\\,v^\{3\}\+\\theta\_\{2\}\\,v\+\\theta\_\{3\}\\,x\\,e^\{\-v^\{2\}\}\+\\theta\_\{4\}\\,x^\{3\}\(12\)\+\\displaystyle\+θ5​v​e−v2\+θ6​x​v2\+θ7​sin⁡\(t\)\\displaystyle\\theta\_\{5\}\\,v\\,e^\{\-v^\{2\}\}\+\\theta\_\{6\}\\,xv^\{2\}\+\\theta\_\{7\}\\,\\sin\(t\)\+\\displaystyle\+θ8​cos⁡\(t\)\+θ9​x​\|x\|\\displaystyle\\theta\_\{8\}\\,\\cos\(t\)\+\\theta\_\{9\}\\,x\|x\|
The linear termθ0​x\\theta\_\{0\}xrepresents a linear restoring force,θ2​v\\theta\_\{2\}vencodes linear damping, andθ1​v3\\theta\_\{1\}v^\{3\}introduces cubic nonlinear damping\. The cross\-termθ3​x​e−v2\\theta\_\{3\}x\\,e^\{\-v^\{2\}\}is a nonlinear coupling between position and velocity, where the Gaussian envelopee−v2e^\{\-v^\{2\}\}localizes the interaction to low\-velocity regimes\. The termθ4​x3\\theta\_\{4\}x^\{3\}adds a nonlinear restoring force of the Duffing\-oscillator type\. The termθ5​v​e−v2\\theta\_\{5\}v\\,e^\{\-v^\{2\}\}is a nonlinear damping contribution similarly localized by the Gaussian factor, attenuating its effect at large velocities\. The higher\-order cross\-termθ6​x​v2\\theta\_\{6\}xv^\{2\}captures nonlinear coupling between position and the square of velocity\. The termsθ7​sin⁡\(t\)\\theta\_\{7\}\\sin\(t\)andθ8​cos⁡\(t\)\\theta\_\{8\}\\cos\(t\)represent external periodic driving forces via sine and cosine components\. Finally,θ9​x​\|x\|\\theta\_\{9\}x\|x\|provides an additional nonlinear restoring force whose magnitude scales withx2x^\{2\}but whose sign tracks the direction of displacement, analogous to quadratic drag in the position domain\.

## 5Conclusion

We presentedGAE, a framework that augments LLM\-guided evolutionary search with three tightly coupled components: a relational GNN encoder that maps typed program graphs into structural embeddings, a contextual bandit selector that learns which program structures yield productive mutations, and a GRPO fine\-tuning loop that adapts the mutation LLM online\. Together, these components close three bottlenecks: reward sparsity, uninformed parent selection, and a static mutation operator\. Experiments on symbolic regression show that GAE achieves the lowest NMSE both in\-distribution and out\-of\-distribution among all baselines, attaining state\-of\-the\-art performance\.

Limitations\.Our method relies on a graph representation derived from Python ASTs, which encodes domain\-specific inductive biases\. Extending to other languages or domains requires redesigning node types and relations, limiting out\-of\-the\-box generalization and introducing non\-trivial adaptation overhead\. The benefit of graph\-guided policy learning depends on a sufficiently expensive evaluation regime, as the GNN and reinforcement optimization introduce additional overhead\.

## Impact Statement

This work provides a fundamental methodological advancement in LLM\-driven evolutionary algorithms\. By introducing a dynamic, structure\-aware paradigm—driven by relational graph encoders and online reinforcement optimization—our framework transforms evolutionary search from a largely stochastic process into a directed, self\-improving trajectory\. This algorithmic evolution is vital for scaling AI to tackle highly complex scientific problems, where the search spaces for viable mathematical or programmatic hypotheses are exceedingly vast and their evaluations computationally expensive\.

The immediate impact of this framework is demonstrated in its ability to automate the discovery of interpretable, closed\-form mathematical models from complex observational data\. This methodology lays the technical groundwork for more comprehensive scientific agent systems designed to solve multifaceted scientific problems across theoretical physics, fluid dynamics, and complex systems engineering\.

We foresee no significant harmful social consequences specific to this work\. The broader social consequences of our work are those that are well established when advancing the fields of Machine Learning and AI for Science, and we do not feel any must be specifically highlighted here\.

## References

## Appendix AGNN Architecture Details

#### Input projection\.

For taskτ\\tau, the one\-hot matrixX∈\{0,1\}N×\|𝒱τ\|X\\in\\\{0,1\\\}^\{N\\times\|\\mathcal\{V\}\_\{\\tau\}\|\}is projected to hidden dimensionhh:h\(0\)=GELU​\(X​Win\)h^\{\(0\)\}=\\mathrm\{GELU\}\(XW\_\{\\mathrm\{in\}\}\)\.

#### Relational convolution\.

For each relationr∈\{1,…,\|ℰ\|\}r\\in\\\{1,\\ldots,\|\\mathcal\{E\}\|\\\}, a dedicatedWr∈ℝh×hW\_\{r\}\\in\\mathbb\{R\}^\{h\\times h\}transforms source features; messages are mean\-aggregated and summed with a self\-loopWselfW\_\{\\mathrm\{self\}\}\. Each layer applies LayerNorm and GELU \(Eq\.[2](https://arxiv.org/html/2607.10127#S3.E2)\)\. This is a*relational GCN*, and no attention weights are computed\.

#### Fallback\.

Any program failing AST parsing produces a singleunknown/unknown\_layernode; the encoder still produces a validz∈ℝ128z\\in\\mathbb\{R\}^\{128\}\.

## Appendix BDecomposability of SAC, GRPO, and the GNN Encoder

The SAC meta\-controller, the GRPO fine\-tuning loop, and the GNN encoder optimize disjoint parameter sets \(\{ϕ,θ1,θ2\}\\\{\\phi,\\theta\_\{1\},\\theta\_\{2\}\\\},θLLM\\theta\_\{\\mathrm\{LLM\}\}, andθGNN\\theta\_\{\\mathrm\{GNN\}\}respectively\)\. No gradients flow between them, and each operates on its own timescale \(per generation, per group, and everyKKgenerations\)\. Disabling any one component degrades gracefully to a simpler baseline:

ComponentParametersTraining signalWithout itSAC meta\-controllerϕ,θ1,θ2\\phi,\\;\\theta\_\{1\},\\theta\_\{2\}r=clip​\(R​\(β⋆\)\)\+λ1​rnov−λ2​rcxr=\\mathrm\{clip\}\(R\(\\beta^\{\\star\}\)\)\+\\lambda\_\{1\}r\_\{\\mathrm\{nov\}\}\-\\lambda\_\{2\}r\_\{\\mathrm\{cx\}\}\(Eq\. 6\)Uniform mutation hintGRPOθLLM\\theta\_\{\\mathrm\{LLM\}\}A^g\\hat\{A\}\_\{g\}fromRggrp=Rg−RαR^\{\\mathrm\{grp\}\}\_\{g\}=R\_\{g\}\-R\_\{\\alpha\}\(Eqs\. 8–9\)Frozen LLMGNN encoderθGNN\\theta\_\{\\mathrm\{GNN\}\}ℒGNN\\mathcal\{L\}\_\{\\mathrm\{GNN\}\}: MSE onℒ\\mathcal\{L\}\(Eq\. 5\)Stale embeddings
Parent selection is uniform over the elite archive in all configurations and is unaffected by any of the three components\.

## Appendix CDiscrete Action Space — Mutation Operators

The Discrete SAC meta\-controller operates over a flat, fixed action space of\|𝒜\|=5×30×18=2700\|\\mathcal\{A\}\|=5\\times 30\\times 18=2700discrete actions, each encoding a tuple \(*mutation type*,*target node id*,*argument index*\)\. The target node id indexes the parent’s expression tree in BFS order \(capped at 30 nodes\), and the argument index selects an operator or constant from the vocabularies below\. The mutation\-type component of the selected action is additionally injected into the LLM prompt as a structural hint \(“prefer<type\>changes”\), steering generation toward the chosen edit class\. The full tuple is used when the mutation is applied directly at the AST level \(the fallback path when no valid LLM completion is produced\)\.

Table 3:Discrete mutation\-type vocabulary for the SAC meta\-controller \(symbolic regression\)\.ActionLabelDescription0add\_operatorWrap a node with a unary op, or insert a binary op above it1replace\_operatorReplace an existing operator node \(arity\-preserving\)2delete\_nodeRemove a unary wrapper, promoting its child3insert\_constantMultiply a sub\-expression by a constant4subtree\_rewriteReplace the node’s subtree with a fresh random treeThe operator vocabulary comprises 5 binary operators \(add,sub,mul,div,pow\) and 13 unary functions \(sin,cos,tan,exp,log,sqrt,abs,tanh,arcsin,arccos,arctan,sinh,cosh\); the constant vocabulary contains 9 values\{−2,−1,−0\.5,0\.5,1,2,3,π,e\}\\\{\-2,\-1,\-0\.5,0\.5,1,2,3,\\pi,e\\\}, giving an argument\-index space of sizemax⁡\(18,9\)=18\\max\(18,9\)=18\. All expressions are built from the three observable inputs: positionxx, velocityvv, and timett\. Numeric constants introduced asparams\[jj\]remain free parameters optimized by BFGS during evaluation\. Generated code applies numerical guards \(clipped exponents,ϵ\\epsilon\-protected division and logarithms\) so that every mutated program remains evaluable\.

## Appendix DAlgorithm

Algorithm 2Discrete SAC for Program Mutation1:GNN encoder

Φ\\Phi, evaluator

ℰ\\mathcal\{E\}, elite population

𝒫\\mathcal\{P\}, prioritized replay buffer

ℬ\\mathcal\{B\}, warm\-up

WW, total steps

SS
2:Initialise policy

πϕ\\pi\_\{\\phi\}, twin critics

Qθ1,Qθ2Q\_\{\\theta\_\{1\}\},Q\_\{\\theta\_\{2\}\}, targets

Qθ¯1,Qθ¯2Q\_\{\\bar\{\\theta\}\_\{1\}\},Q\_\{\\bar\{\\theta\}\_\{2\}\}, temperature

log⁡α\\log\\alpha
3:for

s=1,…,Ss=1,\\ldots,Sdo

4:

αprog∼Uniform​\(𝒫\)\\alpha\_\{\\mathrm\{prog\}\}\\sim\\mathrm\{Uniform\}\(\\mathcal\{P\}\)⊳\\trianglerightparent program, sampled uniformly

5:

z←Φ​\(build\_graph​\(αprog\)\)z\\leftarrow\\Phi\(\\texttt\{build\\\_graph\}\(\\alpha\_\{\\mathrm\{prog\}\}\)\)
6:

a∼\{Uniform​\(𝒜\)if​s<Wπϕ\(⋅∣z\)otherwisea\\sim\\begin\{cases\}\\mathrm\{Uniform\}\(\\mathcal\{A\}\)&\\text\{if \}s<W\\\\ \\pi\_\{\\phi\}\(\\cdot\\mid z\)&\\text\{otherwise\}\\end\{cases\}
7:

β←mutate​\(αprog,decode​\(a\)\)\\beta\\leftarrow\\texttt\{mutate\}\(\\alpha\_\{\\mathrm\{prog\}\},\\texttt\{decode\}\(a\)\)⊳\\trianglerightAST edit viaMutableSRTree

8:

R​\(β\)←ℰ​\(β\)R\(\\beta\)\\leftarrow\\mathcal\{E\}\(\\beta\)⊳\\trianglerightBFGS constant fit\+\+NMSE

9:

z′←Φ​\(build\_graph​\(β\)\)z^\{\\prime\}\\leftarrow\\Phi\(\\texttt\{build\\\_graph\}\(\\beta\)\)
10:

r←clip​\(R​\(β\),rmin,rmax\)\+λ1​rnov​\(z′,𝒫\)−λ2​rcx​\(β\)r\\leftarrow\\mathrm\{clip\}\\\!\\\!\\big\(R\(\\beta\),r\_\{\\min\},r\_\{\\max\}\\big\)\+\\lambda\_\{1\}\\,r\_\{\\mathrm\{nov\}\}\(z^\{\\prime\},\\mathcal\{P\}\)\-\\lambda\_\{2\}\\,r\_\{\\mathrm\{cx\}\}\(\\beta\)
11:Store

\(z,a,r,z′\)\(z,a,r,z^\{\\prime\}\)in

ℬ\\mathcal\{B\}
12:if

\|ℬ\|≥32\|\\mathcal\{B\}\|\\geq 32and

s≥Ws\\geq Wthen

13:Sample prioritized minibatch; update

Qθ1,Qθ2Q\_\{\\theta\_\{1\}\},Q\_\{\\theta\_\{2\}\},

πϕ\\pi\_\{\\phi\}, and

α\\alpha
14:Soft\-update targets

Qθ¯1,Qθ¯2Q\_\{\\bar\{\\theta\}\_\{1\}\},Q\_\{\\bar\{\\theta\}\_\{2\}\}with Polyak coefficient

τ\\tau
15:endif

16:

𝒫\.Insert​\(β,R​\(β\),z′\)\\mathcal\{P\}\.\\mathrm\{Insert\}\(\\beta,R\(\\beta\),z^\{\\prime\}\)⊳\\trianglerightkeep top\-NNby score

17:endfor

## Appendix EHyperparameters

Table[4](https://arxiv.org/html/2607.10127#A5.T4)collects all parameters introduced in the methodology \(§3\.2–§3\.5\), grouped by component\. Values marked†are set in the implementation but not stated explicitly in the main text\.

Table 4:Hyperparameters of GAE, grouped by component\.ParameterMeaningValue*Elite archive and evolution loop \(§3\.2, Alg\. 1\)*NNArchive size \(top\-NNelite population\)5050GGChildren generated per generation \(group size\)88KKGNN fine\-tuning interval \(generations\)33ddStructural embedding dimension128128*GNN encoder \(§3\.3\)*FFNode\-type vocabulary size \(one\-hot input\)1717\|ℰ\|\|\\mathcal\{E\}\|Number of typed edge relations66LLRelational message\-passing layers22hhHidden width \(pooled projection gives2​h=d2h=d\)6464—Gradient steps per fine\-tuning round300300—Sliding\-window size ofℒ\\mathcal\{L\}\(recent pairs\)300300—Minimum samples before fine\-tuning \(\|ℒ\|\|\\mathcal\{L\}\|\)1010*Discrete SAC meta\-controller \(§3\.4\)*\|𝒜\|\|\\mathcal\{A\}\|Action\-space size\(5×30×18\)\(5\\times 30\\times 18\)27002700—Policy/critic MLP hidden width256256τ\\tauPolyak target\-update coefficient0\.0050\.005WWWarm\-up steps \(uniform actions\)3232—Minimum buffer size before updates3232—Replay buffer capacity \(prioritized\)10,00010\{,\}000—Mini\-batch size6464γ\\gammaDiscount factor in the TD target0\.99†0\.99^\{\\dagger\}—SAC learning rate \(Adam, all networks\)3×10−4†3\\times 10^\{\-4\}\\,\{\}^\{\\dagger\}ℋ¯\\bar\{\\mathcal\{H\}\}Target entropy \(auto\-tunedα\\alpha\)−12​log⁡\|𝒜\|\-\\tfrac\{1\}\{2\}\\log\|\\mathcal\{A\}\|rmin,rmaxr\_\{\\min\},\\,r\_\{\\max\}Task\-score clipping bounds \(Eq\. 6\)−5,15\-5,\\;15λ1\\lambda\_\{1\}Novelty bonus weight \(Eq\. 6\)0\.10\.1λ2\\lambda\_\{2\}Complexity penalty weight \(Eq\. 6\)0\.050\.05—Complexity penalty onset \(AST nodes, Eq\. 6\)2020*GRPO LLM fine\-tuning \(§3\.5\)*ϵclip\\epsilon\_\{\\mathrm\{clip\}\}PPO clipping range \(Eq\. 8\)0\.20\.2βkl\\beta\_\{\\mathrm\{kl\}\}KL coefficient \(monitoring only, Eq\. 8\)0\.10\.1—Crash\-sentinel reward threshold \(dropped\)−10\-10—Advantage\-std skip threshold10−810^\{\-8\}