LLM-SR: Scientific Equation Discovery via Programming with Large Language Models
Symbolic regression becomes far more sample-efficient when an LLM proposes executable equation skeletons conditioned on a natural-language problem specification while numerical optimizers fit the placeholder parameters, with a FunSearch-style islands experience buffer maintaining diverse high-scoring hypotheses.

Mar 2025 · Paper
1. Motivation
Symbolic regression (SR) - recovering closed-form equations from data - is a core scientific-discovery task, but classic genetic-programming SR searches expression trees essentially blindly, requiring millions of candidate evaluations, while large pre-trained SR transformers generalize poorly outside their training distributions.1 LLMs carry substantial latent scientific prior knowledge, yet using them to generate full equations with numeric constants invites hallucinated parameters and physically meaningless expressions.1 LLM-SR’s key move is a decoupling: let the LLM do what it is good at (proposing plausible functional structure from domain priors) and let numerical optimizers do what they are good at (fitting constants to data). It also attacks a subtler problem - standard SR benchmarks such as the Feynman lectures are substantially memorized by LLMs - by building anti-recitation benchmarks.1
2. Methodology
2.1 Search Structure
Evolutionary program search over equation skeletons, organized as an islands model (independently evolving subpopulations, after PySR and FunSearch).12 Each iteration :
1. Sample an island uniformly
2. Sample k=2 in-context programs from that island (Boltzmann over clusters)
3. Prompt LLM (temperature 0.8) → batch of b=4 new skeletons
4. For each skeleton: optimize placeholder params → score (NMSE)
5. Update source island if a new best is found
LLM-SR runs ~2.5K iterations versus >2M candidate evaluations allowed for search-based baselines - the sample-efficiency claim is measured against that asymmetry.1
2.2 State Representation
A hypothesis is an executable Python program of the form def f(x, params): ... return y - the equation structure in code, with all numeric constants abstracted into a parameter vector params.1 The experience buffer stores (skeleton, score) pairs per island; within islands, programs are clustered by score “signature” to preserve diversity. Benchmarks ship with natural-language problem specifications (variables, constraints, objectives) that inject scientific priors into the prompt.1
2.3 Generation Mechanism
The prompt has four components: (i) an Instruction (complete the function body, respect physical meanings of variables); (ii) the Problem Specification; (iii) the Evaluation/Optimization function code; (iv) Experience Demonstrations - sorted skeleton programs showing an improvement trajectory, i.e., FunSearch-style best-shot prompting.1 The LLM performs mutation and crossover implicitly through in-context learning over these demonstrations; there is no explicit variation operator. The decoupled two-stage design is strict: the LLM proposes structure only, and parameters are then fitted by numpy + BFGS (scipy) or torch + Adam, with a 30-second optimization timeout per candidate.1
2.4 Evaluation & Grounding
Grounding is fully data-driven: score = normalized mean-squared error (NMSE) of the parameter-optimized program against observed data, reported on both in-domain (ID) and held-out out-of-domain (OOD) ranges - the OOD split detects overfit or physically meaningless expressions.1 Four purpose-built benchmarks (Oscillation 1, Oscillation 2 - nonlinear damped oscillators; E. coli growth; Stress–Strain) were designed to prevent LLM recitation, validated by showing higher LLM perplexity and slower discovery curves on them than on Feynman-Lecture problems, which the authors demonstrate are substantially memorized.1
2.5 Memory / Information Flow
The islands experience buffer is the search’s episodic memory:
new (skeleton, score) ──beats island best?──► admitted to source island
prompt construction: island ~ Uniform; program ~ Boltzmann over
score-signature clusters, P_i ∝ exp(s_i / τ_c)
Islands evolve independently from a trivial linear template (equation_v0); diversity is preserved by isolation plus cluster-based sampling rather than by migration (contrast FunSearch’s periodic culling of the worst islands).12
2.6 Selection & Propagation
Three nested selection rules: improvement-gated admission to the buffer (elitist per island); Boltzmann/softmax sampling favoring higher scores within islands for prompt construction; and return of the best-scoring program overall at termination.1 Selection thus operates on evaluated artifacts, not verbal reflection - placing LLM-SR on the artifact-memory side of the field’s memory divide, opposite ReEvo’s reflection texts.3
2.7 Anti-Slop / Anti-Hallucination Mechanisms
Four structural controls: (i) the skeleton+optimizer decoupling makes parameter hallucination impossible by construction - the ablation “w/o skeleton+optimizer,” where the LLM emits full equations with numbers, degrades catastrophically (ID NMSE $3.78 \times 10^{-1} vs \2.12 \times 10^{-7}$ for the full model on Oscillation 2); (ii) anti-recitation benchmark design with perplexity analysis documenting Feynman-benchmark memorization; (iii) execution plus numerical optimization grounds every hypothesis before it can enter memory; (iv) OOD evaluation exposes meaningless expressions.1
3. Results
NMSE (↓), Table 1, GPT-3.5 backbone, ~2.5K iterations vs 2M+ for baselines:1
| Model | Osc1 OOD | Osc2 ID | Osc2 OOD | E.coli OOD | Stress OOD |
|---|---|---|---|---|---|
| GPlearn | 0.5567 | 0.7551 | 3.188 | 1.039 | 0.4091 |
| E2E | 0.3722 | 0.1401 | 0.1911 | 1.4467 | 0.5867 |
| PySR | 0.3106 | 0.0002 | 0.0098 | 1.0141 | 0.1304 |
| uDSR | 0.0007 | 0.0032 | 0.0015 | 5.4584 | 0.1761 |
| LLM-SR (GPT-3.5) | 0.0005 | 2.12e-7 | 3.81e-5 | 0.0264 | 0.0516 |
| LLM-SR (Mixtral) | 0.0002 | 0.0030 | 0.0291 | 0.0037 | 0.0946 |
LLM-SR shows sharp error drops within its first iterations while baselines improve gradually, and the gap widens over time - evidence that prior-guided LLM structure proposal beats blind genetic-programming exploration per evaluation.1 Ablations on Oscillation 2: removing the natural-language prior, iterative refinement (ID $1.01 \times 10^{-1}), or the skeleton+optimizer decoupling (\3.78 \times 10^{-1}$) each substantially degrades performance; numpy+BFGS slightly beats torch+Adam, attributed to the LLM’s higher numpy proficiency.1
4. Limitations & Assumptions
Performance is “inherently tied to the quality and breadth of the LLM’s training data,” and iterative LLM queries plus optimization “could be prohibitive for large-scale problems.”1 The method assumes an informative natural-language problem specification exists (the w/o-Prior ablation shows real dependence). The 30s optimizer timeout caps equation complexity, and benchmarks use a single independent variable per problem. Notably, the SOTA claim has been complicated by the authors’ own follow-up benchmark LLM-SRBench (ICML 2025) and by a 2026 critique, “Set-Level Selection Beats LLM Evolution in Scientific Equation Discovery.”45
5. Critical Analysis
- Distinctive strength: the skeleton/optimizer decoupling is a clean, ablation-proven solution to LLM numeric hallucination, and the anti-recitation benchmark methodology (perplexity + discovery-curve analysis) is a transferable anti-slop instrument for any LLM-for-science claim.1
- Structural weakness / trade-off: island admission is gated on beating the island best, so memory keeps few stepping-stones; combined with single-variable, small-timeout benchmarks, it remains unclear whether the approach scales to high-dimensional, multi-physics problems where parameter fitting itself is the bottleneck.1
- Connections: LLM-SR explicitly adopts FunSearch’s islands model and best-shot demonstration prompts, adding scientific-prior prompting and parameter decoupling.12 Its score-signature clustering within islands is a lightweight quality-diversity mechanism echoing MAP-Elites’ archive-based illumination,6 and relative to ReEvo it marks the opposite memory design - implicit in-context crossover over evaluated artifacts versus explicit verbal reflection.3 AlphaEvolve generalizes the same program-search-with-evaluator architecture to full programs at datacenter scale.7
6. References
Footnotes
-
Shojaee et al., LLM-SR, ICLR 2025 (HTML v3): https://arxiv.org/html/2404.18400v3 · OpenReview: https://openreview.net/forum?id=m2nmp8P5in · code: https://github.com/deep-symbolic-mathematics/LLM-SR ↩ ↩2 ↩3 ↩4 ↩5 ↩6 ↩7 ↩8 ↩9 ↩10 ↩11 ↩12 ↩13 ↩14 ↩15 ↩16 ↩17 ↩18 ↩19 ↩20 ↩21
-
Romera-Paredes et al., FunSearch, Nature 625:468–475: https://www.nature.com/articles/s41586-023-06924-6 ↩ ↩2 ↩3
-
Ye et al., ReEvo, NeurIPS 2024: https://arxiv.org/abs/2402.01145 ↩ ↩2
-
Shojaee et al., LLM-SRBench, ICML 2025: https://arxiv.org/abs/2504.10415 ↩
-
Set-Level Selection Beats LLM Evolution in Scientific Equation Discovery (2026): https://arxiv.org/html/2607.04108v1 ↩
-
Mouret & Clune, MAP-Elites: https://arxiv.org/abs/1504.04909 ↩
-
Novikov et al., AlphaEvolve: https://arxiv.org/abs/2506.13131 ↩