LLM-SR: Scientific Equation Discovery via Programming with Large Language Models
Authors: Parshin Shojaee
Kazem Meidani
Shashank Gupta
Amir Barati Farimani
Chandan K. Reddy
Venue: ICLR 2025 Oral (2025)
Keywords: Symbolic Regression
Equation Discovery
Large Language Models
Evolutionary Search
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We introduce LLM-SR, an approach that harnesses Large Language Models (LLMs) to discover governing equations from data in an efficient, knowledge-guided manner.

Abstract: Mathematical equations have been unreasonably effective in describing complex natural phenomena across various scientific disciplines. However, discovering such insightful equations from data presents significant challenges due to the necessity of navigating extremely large combinatorial hypothesis spaces. Current methods of equation discovery, commonly known as symbolic regression techniques, largely focus on extracting equations from data alone, often neglecting the domain-specific prior knowledge that scientists typically depend on. They also employ limited representations such as expression trees, constraining the search space and expressiveness of equations. To bridge this gap, we introduce LLM-SR, a novel approach that leverages the extensive scientific knowledge and robust code generation capabilities of Large Language Models (LLMs) to discover scientific equations from data. Specifically, LLM-SR treats equations as programs with mathematical operators and combines LLMs' scientific priors with evolutionary search over equation programs. The LLM iteratively proposes new equation skeleton hypotheses, drawing from its domain knowledge, which are then optimized against data to estimate parameters. We evaluate LLM-SR on four benchmark problems across diverse scientific domains (e.g., physics, biology), which we carefully designed to simulate the discovery process and prevent LLM recitation. Our results demonstrate that LLM-SR discovers physically accurate equations that significantly outperform state-of-the-art symbolic regression baselines, particularly in out-of-domain test settings. We also show that LLM-SR's incorporation of scientific priors enables more efficient equation space exploration than the baselines.

Published as a conference paper at ICLR 2025
LLM-SR: SCIENTIFIC EQUATION DISCOVERY VIA PROGRAMMING WITH LARGE LANGUAGE MODELS
Parshin Shojaee'
Kazem M
Shashank Gupta²
Amir Barati Farimani
Chandan K. Reddy
Virginia Tech
Carnegie Mellon University
Allen Institute for Al
ABSTRACT
Mathematical equations have been unreasonably effective in describing complex natural phenomena across various scientific disciplines. However, discovering such insightful equations from data presents significant challenges due the necessity of navigating extremely large combinatorial hypothesis spaces. Current methods of equation discovery, commonly known as symbolic regression techniques, largely focus extracting equations from data alone, often neglecting the nomain-specifft prior knowledge that scientists typically depend on. They also employ limited representations such expression trees, constraining the search space and expres- siveness of equations. To bridge this gap, we introduce LLM-SR, a novel approach that leverages the extensive scientífic knowledge and robust code generation capa- bilities of Large Language Models (LLMs) to discover scientific equations from data. Specifically, LLM-SR treats equations as programs with mathematical opera- tors and combines LLMs' scientific priors with evolutionary search over equation programs. The LLM iteratively proposes new equation skeleton hypotheses, draw- ing from its domain knowledge, which are then optimized against data to estimate parameters. We evaluate L.LM-SR on four benchmark problems across diverse scientific domains (e.g., physics, biology), which we carefully designed to simulate the discovery process and prevent LLM recitation. Our results demonstrate that LLM-SR discovers physically accurate equations that significantly outperform state-of-ther -frt--r-tt---ttttee symbolic regression baselines, particularly in out-of-domain test settings. We also show that LLM-SR's incorporation of scientific priors enables more efficient equation space exploration than the baselines [I]
INTRODUCTION
The emergence of Large Language Models (LLMs) has marked a significant milestone in artificial intelligence, showcasing remarkable capabilities across various domains (Achiam et al. 2023). As LLMs continue to evolve, researchers are exploring innovative ways harness their potential for solving complex problems such as scientific discovery (Wang et al. 2023a A14Science & Quantum 2023). Their ability to process and comprehend vast amounts of scientific literature, extract relevant information, and generate coherent hypotheses has recently opened up new avenues for accelerating scientific progress (Zheng et al.) 2023b; Ji et al. 2024). Additionally, by leveraging their ability to understand and reason with the help of programming and execution, L.Mss have shown the potential to enhance automated reasoning and problem-solving capabilities for general natural language and mathematics optimization tasks, e.g., prompt optimization and heuristic discovery (Meyerson et al. 2023 Yang et al.) 2023al Madaan et al.) 2024 Romera-Paredes et al. 2024). Motivated by these strengths, LLMs could be particularly helpful for the task of equation discovery, a fundamental task in science and scientific discovery.
Discovering accurate symbolic mathematical models from data is an important task in various scientific and engineering disciplines. The task of data-driven equation discovery (also commonly known as Symbolic Regression (SR)), aims to find abstract mathematical equations from data
"Equal contribution. Contact: parshinshojace@wr.chh mmcidani ammeidani@@anddew.cm....h
Code and data available: htttp:://// thub. wtwpbcoooooo...cccccccocc...ooot deep-
htmpm//eeppppemlll11mmmmm symbo1 ic-mathematttt..... C8/1.IMSSK

---

Published as a conference paper at ICLR 2025
Equation Skeleton Hypothesis as Pogram
Insinslum
Hypothesis Generation
Speciitaatio
breleetion:
Prompt
LLM
Eqersonne Somestrratin:
pa wwwppaa/pppaaaaaaaaaaa (1P
waleiin paran[[]]
+b
Fetch Experience
Experience Management
Data-driven
Dalaset
E vvllatioo
plimize Parameters (Torchnddos 15RP)
Experience Manager
-6.43
Evaluate
paronc 2.1
burre
Experience Buffer
parroee
LLM-SR framework, consisting three main steps: (a) Hypothesis Generation, where LLM generates equation program skeletons hased on a structured prompt; (b) Data-driven Evaluation, which optimizes the parameters of cach equation skeleton hypothesis and assesses its fit to the data; and (c) Experience Management, which maintains diverse buffer of high-scoring hypotheses to provide informative in-context examples into LLM's prompt for effective iterative refinement.
Figure
observations such that these equations predictive of the underlying data, are interpretable, and generalize to unseen data from the same physical phenomena. Finding such equations offers several advantages over simply estimating predictive model, the resulting mathematical functions provide insights into the underlying physical processes, enable extrapolation beyond the observed data, and facilitate knowledge transfer across related problems (Langley 1981 Schmidt & Lipson 2009). However, while evaluating the fit of a proposed equation is relatively straightforward, the inverse process of obtaining these mathematical equations from data challenging problem, known to be NP-hard (Virgolin & Pissis) 2022). Current equation discovery methods encompass a wide variety approaches from evolutionary search algorithms (Cranmer 2023 Mundhenk al. 2021
La Cava et al., 2021) to advanced deep learning methods using Transformers (Biggio et al. 202
Kamienny et al.] 2022). Most of the traditional symbolic regression techniques are built on top of Genetic Programming (GP) (Koza 1994) evolutionary methods, representing mathematical equations as expression trees and searching the combinatorial space of possible equations through iterative mutation and recombination. However, these methods often struggle with the complexity of the vast optimization space and do not incorporate prior scientific knowledge, which leads to suboptimal solutions and inefficient exploration of the equation search space. Similarly, the design of current general LLa-b optimization frameworks (Meyerson et al.) 2023] Rooera-PPaedee et al.. 2024 Yang et al. 2023a) also have several key limitations terms of domain knowledge integration and diverse exploration which are critical for equation discovery. Thus, there is need for specialized equation discovery methods that effectively integrate prior scientific knowledge into the navigation of vast equation search space, strategy akin to a scientist's 's reliance on foundational scientific knowledge when formulating hypotheses for scientific discovery.
To address all these limitations, we introduce LLM-SR (shown in Fig. T. novel framework that combines the strengths of L.LMs, reliable optimizers, and evolutionary search for data-driven equation discovery. At its core, LLM-SR is an iterative hypotheses refinement method that generates, evaluates, and refines equation hypotheses based on data-driven feedback. Specifically, LLM-SR first prompts the LLM to propose new equation hypotheses (Fig. [a)), then evaluates their fit on the observed data using off- the-hhh--heff optimizers (Fig. [b)), and uses this data-driven feedback and a carefully maintained dynamic memory of previous equations (Fig. [(c) iteratively guide the search towards better equations. LLM-SR leverages the scientific knowledge embedded in LLMs using
short descriptions of the problem and the variables involved in given system to generate educated
2

---

Published as a conference paper at ICLR 2025
hypotheses for equation skeletons (i.e., mathematical structures with placeholder parameters for numeric coefficients and constants). The LLM's in-context learning and crossover capabilities (Meyerson et al., 2023) are then employed to refine the suggested equation skeletonss in iterative process. By representing equations as Python programs, we take advantage of I LLM's ability to generate structured and executable code (Li et al. 2023 Shhoaaee et al. (023) while providing flexible and effective way to represent general mathematical relations. The program representation also f'ciilitates direct and differentiable parameter optimization to better optimize the coefficients or constants in the generated equations.
To leverage LLM's scientific prior knowledge yet prevent the risk of LLM recitation (Wu al. 2023) in equation discovery (observed for common benchmarks like Feynman (Udrescu & Tegmark 2020), we designed four custom benchmark problems across physies, biology, and materials science for the evaluation of LLM-SR. By incorporating synthetic modifications to physical models and experimental datasets, these problems aim to simulate the real discovery processes (see App. C and for details). We evaluated LLM-SR using GPT-3.5-turbo (Brown 2020) and Mixtral-8x7BB (Jiang et al. 2024) as backbone LLMs. Results demonstrate that LLM-SR consistently outperforms state-of- the-art art symbolic regression methods, discovering physically accurate equations with better fit and generalization in both in-domain (ID) and out-of- (OOD) test settings. By leveraging the scientific prior knowledge, LLM-SR explores the equation search space more efficiently, requiring fewer iterations to find accurate equations. Our ablation analysis also highlights the crucial role of data-driven feedback, iterative refinement, and program representation LLM-SR's performance. The major contributions of this work are as follows:
We introduce LLM-SR, a novel framework that leverages domain-specific prior knowledge and code generation capabilities of LLMs combined with off-the-shelf optimizers and evolutionary search for data-driven scientific equation discovery.
We create four benchmark problems spanning physics, biology, and materials science, designed to simulate real-world discovery and prevent LLM recitation risks for evaluation of LLM-SR.
We show that LLM-SR outperforms state-of-the-ar symbolic regression methods by navigating the equation search space more efficiently and discovering more accurate equations with better out-of-domain generalization.
We demonstrate through a comprehensive ablation study that natural language problem descriptions. program representation, data-driven feedback, and iterative hypothesis refinement are all essential components for LLM-SR's success.
2
LLM-SR METHODOLOGY
2.1 PROBLEM FORMULATION
In the task of data-driven equation discovery, also known as symbolic regression (SR), the goal is to find concise symbolic expression j aproxximaannnggg an unknown function / Rd + R. Given dataset D= = {(x,, H/ SR methods seek to uncover the hidden mathematical relationship such that f(x,) RJ yı. Vi. The discovered equation should not only accurately fit the observed data points but also exhibit strong generalization capabilities to unseen data while maintaining interpretability.
Current SR methods typically represent equations using techniques such as expression trees (Cranmer) 2023), prefix sequences (Petersen et al., 2021 Biggio et al. (022]). or context-free free grammars (Brence al. 2021). These representations provide structured and constrained search spaces, enabling evolu- tionary algorithms like genetic programming to explore and find candidate expressions. In contrast, we employ program functions to directly map inputs æ targets y: def f (x) This approach offers greater expressiveness in mathematical relations but expands the search space significantly. To navigate this vast program space effectively, we leverage LLMs for their scientific knowledge and code generation capabilities. Let wa denote a pre-trained LLM with parameters 0. We iteratively sample equation program skeletons J {f: f ]e). aiming to maximize the reward Scorey (f. D) for a given scientific problem T and dataset D: f = arg maxy EdeD [Scorey (f. D)). Our approach, I.LM-SR, prompts the I.LM to propose hypotheses based on problem specifications and demonstrations of previously discovered promising equations. The I.I.M generates equation program skeletons with placeholder parameters, which are then optimized using robust. Python optimizers. Promising hypotheses are added to dynamic experience buffer, guiding subsequent in-context example updates and equation refinement. Below we explain the key components of this framework, shown in Fig. I in more detail.
return y.
3

---

Published as a conference paper at ICLR 2025
Instruction
Vou are helpfut acst stant tacked with stsrovening sutheratiral fanctios structures scientáfic systerg. complete the "squation" below, considaring the physical completing
furction. letts think
Problem Specificttion
Fied nathomatical function akelaton that represents ecceteratios donped nonlinear oset llttor syaten with driving forcs, given date on poakison, seloesty, and siee.
Evalustion and Optimization
def (aata equation):
Purtrler Optnialonn NamPy)
Ballsae equatson dats ehservatsime.
Eititnia partat Hind dis Deiyy
datar dietsinory containing "ingets" aud "ostguta' errage
sqnatsin: A Jsncison representssy the squatios to be evelosted.
Acture:
befsces prodicted outysts.
Load dais observeliioo
inpats, cutputa datal isppts")) ), data ('oepputpt
imputs
Opttntze peroneters bared data
lepert serd
Return mvalaction score retura score
def witth.. (Nll))
Equation Program Example
ait tormetlani,, it:
def egintiio_sriie (t, parsex): Rachesstiesi fuactson for acceleration in dasped nosisecar asositaton
qptasiue bedt optca Marland pussstacsl), leto
Arga:
1: Pine
spticium giall lamm tedvere))
Position
Velocity
panses: Runorsc parsecters to be optsnized
(nme))
doselerasson are
instial. return paruaa [0]+t pazaza [1]*x parana [2] paruna [3]
/saction structure
Function Complete
def squatioo, vi (s, R, Va paraza): "e>Zeproudd verssion of equation_a@*s
Figure 2: Example of initial input prompt for the nonlinear oscillator discovery task, including problem specification, evaluation and optimization function, and the initial input equation example.
2.2 HYPOTHESIS GENERATION
The hypothesis generation step (Fig. aa) utilizes pre-trained LLM to propose diverse and promising equation program skeletons. Our prompt structure, shown in Fig. consists of the following compo- nents: Instruction: A clear directive for completing the function body, emphasizing consideration of physical meanings and relationships among input variables; Problem Specification: concise description of the scientific problem, including key variables, constraints, and objectives; Evaluation and Optimization Function: The function used to assess the data-driven quality and fitness of pro- posed equations after parameter optimization; and Experience Demonstration: In-context examples of equation skeleton programs and their improvement trajectory.
At each iteration t, we sample a batch of equation skeletons =;= (f from the LLM meद fo Te(-)p:) where Pt is the constructed prompt. We employ stochastic cemperature-basee sampling to balance exploration (creativity) and exploitation (prior knowledge) in the hypothesis space. Sampled equation programs are executed, and those failing to execute or exceeding a maximum execution time threshold are discarded to ensure validity and computational efficiency.
2.3 HYPOTHESIS OPTIMIZATION AND ASSESSMENT
After generating equation skeleton hypotheses, we evaluate and score them using observed data (Fig. [b)). This process involves optimizing the parameters of each hypothesis and then assessing its fitness. We decouple equation discovery into two steps: (1) discovering the equation program structures (skeletons) using the I.M,, and (ii) optimizing the skeleton parameters/cooffiiients based on data. The LLM is responsible for generating equation skeletons and the core logic the program, while the numeric values of the parameters are represented as placeholders in the form of parameter vector par ams (as shown in Fig. 2). These placeholders are subsequently optimized to fit the data. Each equation program skeleton f E F, is function of the form: "def E (x, parans):
return y". We employ two optimization approaches: numpy+BFGS: nonlinear optimization method using sci py library (Fletcher 1987), and torch+Adam: A stochastic gradient-based
4

---

Published as a conference paper at ICLR 2025
optimization algorithm using PyTorch (Kingma & Ba. 2014). The choice between these methods depends on the problem characteristics and equation skeleton complexity. The numpy+BFGS is preferred for problems with fewer parameters, while Adam is more suitable for larger-scale problems benefiting from efficient gradient computation through differentiable programming.
After optimizing the skeleton parameters (params"), we assess the fitness of equation program hypotheses by measuring its ability to capture underlying patterns in the data. We compute predicted target values as: y [x, params"). The fitness evaluation score 8 is then calculated as the negative Mean Squared Error (MSE) between predicted and true target values: & = Sc,oer([1,2) -MSE y).
2.4 EXPERIENCE MANAGEMENT
To efficiently navigate the search landscape and avoid local minima, LLM-SR employs an experience management step (Fig ()). This process maintains diverse population of high-quality equation programs in a dynamic experience buffer and samples from this population to construct informative prompts for subsequent LLM iterations. Let Pa denote the experience buffer at iteration t, storing pairs of equation skeleton hypotheses and their corresponding scores s). We adopt an islands model (Cranmer, 2023 Romera- Paredes et al.] 2024) with independently evolving islands, initialized with a copy of the equation program example from the initial prompt (equat t on_vv in Fig. (2). At each iteration I, new hypotheses F, and their scores are added to the source island (from which the in-context examples of prompts were sampled) if they improve upon the current best; Pj + Pi U {(f f E Fr. s S D), shea] where Pi the i-th island and Sheat is its current best score. Within each island, equation programs are clustered based on their signature (defined by their score) to further preserve diversity.
To construct informative prompts, we then sample equation programs from the experience buffer using a two-stage method. First, uniformly select a random island from the m available. Second, sample & equation programs from the selected island using (a) Cluster selection via Boltzmann sampling, favoring higher scores: P: Σ. (x/=)) where 8; is the mean score of the i-th cluster
exp(s./T.)
and Te is a temperature parameter. (b) Individual program sampling, favoring shorter programs: P(f:) x exp(-1./Tp) where I is the normalized program length and Tp is a temperature parameter. The sampled programs are then included in the prompt as in-context experience demonstrations, guiding the LLM generating new equation program hypotheses. Detailed sampling procedures are provided in Appendix B
Algorithm | presents the simplified pseudo-code of the LLM-SR framework. The experience buffer Po is initial- ized with initial prompt. using a simple linear equation skeleton a template (c.g.. Fig. for the nonlinear oscil- lator problem). This initial structure serves as a baseline for the LLM to modify operators and structures based its domain knowledge. Each iteration ł involves: (i) sam- pling & in-context examples from P1-1. (60) updating the prompt, (iii) generating b equation program skeletons from the LLM, and (tv) evaluating and potentially adding these to P. if they improve upon the best score s". This process leverages the [.M'' generative capabilities to re- fine equation structures guided by the evolving experience buffer. The algorithm returns the best-scoring program f" and its score s* as the optimal solution, iteratively ex- ploring the equation space while balancing exploitation of promising structures with exploration of new possibilities.
Algorithm 1: LLM-SR
Input :LLM Ta dataset D. problem T, T itteraaions, & in-context
examples,
samples per prompt
# Initialize population Po InitPop()
f + null,-oo -o0 for t + to 7 do
# Sample examples from buffer
< {c,1ー=1-
SampieExp(P:r-a)
# Prompt with new examples
p + MakkeeesShooPPrmpp((hc (B)
# Sample from LLM
Fi (( fj ~ TT(-Ip)
# Evaluation and population update
for f F,
(corer(f.D) if + > 8 then
8" f,8 P
3
EXPERIMENTS
end
3.1 BENCHMARKS AND DATASETS
end
The Feynman benchmark (Udrescu & TTemarkk 2020), comprising 120 fundamental physics problems from Feyn- man Lectures on database
end Output :/" S
database s is the current
5

---

Published as a conference paper at ICLR 2025
standard benchmark for evaluating symbolic regression techniques in scientific equation discov- ery. However, our investigation reveals that LLMs have significant memorization issues with these well-known physics equations, potentially undermining their effectiveness in assessing LLM-based equation discovery approaches. For instance, L.LMSRR rapidly achieves low data-driven errors within few iterations (< 20) on Feynman problems, suggesting a recitation of memorized information rather than a process of discovery (full results in App. C). To address these limitations and provide a more robust evaluation, we introduce novel benchmark problems across three scientific domains. Our benchmark design focuses on: (i) Custom modifications to physical models to prevent trivial memorization; (tt) Complex synthetic non-linee oon-linerr relationships to challenge creative exploration; and (ii) Realistic scenarios with experimental data to reflect real modeling processes. These benchmarks are designed to challenge the model's ability to uncover complex mathematical relations while leveraging its scientific prior knowledge, simulating conditions for scientific discovery. To validate our new benchmarks' effectiveness, we compared LLM response perplexity (using open-source Mixtral-8x 8x7B) and equation discovery error curves (using GPT-3.5) between Feynman problems and our new benchmarks. Results show lower perplexity (Fig. 9 in App. C and sharper discovery curves (Fig. in App. C for Feynman problems, suggesting that both LLM backbones have more likely memorized common Feynman equations, while our benchmarks present novel challenges requiring reasoning and exploration. We next discuss these new benchmark problems in detail:
Nonlinear Oscillators Nonlinear damped oscillators, ubiquitous in physics and engineering, governed by differential equations describing the complex interplay between an oscillator's position, velocity, and acting forces. The general form of these equations is typically expressed as: H f(1, .) = 0 where / time, x is position, and f(4, x i)) represents nonlinear forces. To challenge LLM-based equation discovery methods beyond common oscillator systems (e.g.. Van der Pol, Duffing), we introduce two custom nonlinear designs: ! Poooillattt
1: ú = F sin((wr) ov3 3 yrv - x cos(z); and ② Osci socillatioo 2: é = F sin(wit) ax.3 fFrr 8x exp()x)) where v = è represents velocity, and w, a. 3, Y, 5 are constants. These two forms, serving as a proof of concept, are carefully designed to incorporate a combination of challenging yet solvable nonlinear structures (including trigonometric, polynomial, and exponential) that are distinct from well-known oscillator systems. More details on the design rationale and data generation are provided in App. D.I
The growth of Escherichia coli (E. coli) bacteria has been widely studied in microbiology due to its importance in various applications, such as biotechnology. and food safety. Discovering equations governing E. coli growth rate under different conditions is crucial for predicting and optimizing bacterial growth. The bacterial population growth rate has been modeled using a differential equation with the effects of population density (B), substrate concentration (S), temperature (T), and pH level, which is commonly formulated with multiplicative structure: dB (B, T, pH) = f& [B) fS( fT( (pu( To create a challenging benchmark that leverages LLMs' prior knowledge while preventing trivial memorization, we introduce novel nonlinear formulations for Sr(T) and (pH) These custom functions maintain key characteristics of established models while introducing complexities that require exploration and discovery rather than recall. The complete mathematical formulations, along with the data generation process and parameter ranges, are detailed in App. D.2]
Bacterial Growth
Material Stress Behavior The stress-strain relationship of materials under varying conditions, particularly as function of temperature and material type, is fundamental to structural design and analysis across engineering disciplines. This benchmark problem leverages a real-worid experimental dataset from (Aakash et al. 2019), comprising tensile tests on Aluminum 6061-T651 across range of temperatures. The inclusion of this benchmark serves multiple purposes: (i) It challenges LLsed equation discovery methods with experimental data, moving beyond synthetic or idealized problems. (M) Unlike the previous benchmarks, there is no predetermined theoretical model structure for this problem, necessitating creative modeling approaches from LLMs. In other words, modeling for this type of task is mostly empirical and the stress-straine strain-tempprature relations may vary significantly based the specific material and experimental condition, preventing trivial memorization. More details on this problem and experimental data are provided in App. D33
3.2 EXPERIMENTAL SETUP
We compare LLM-SR against state-of-the-an symbolic regression (SR) methods, including evolutionary-based approaches like GPlearr (Genetic Programming) and PySR (multi-island
https://githubs http://wwww..ooommioaaccaaaaeeee

---

Published as a conference paper at ICLR 2025
Oscillation
Oscillation 2
E. coli growth
Stress-Strain
Model
IDI
0OD
IDI
OOD1
IDL
OODL
IDJ
0ODJ
GPlearn
0.0155
0.5567
0.7551
3.188
1,081
1.039
0.1063
0.4091
NeSymReS (Biggio al. 2021
0.0047
0.5377
0.2488
0.6472
N/A (d
>
0.7928
0.6377
E2E
(Kamiehny et 2024
0.0082
0.3722
0.1401
0.1911
0.6321
1.4467
0.2262
0,5867
DSR (Ptterse et al. 2021
0.0087
0.2454
0.0580
0.1945
0.9451
2.4291
0.3326
1,108
uDSK
Landajuela
2022
0.0003
0.0007
0.0032
0.0015
0.3322
5.4584
0.0502
0.1761
PySR
Crammer
0.0009
0.3106
0.0002
0.0098
0.0376
1.0141
0.0331
0.1304
LLM-SR (Mixiral)
7.89e-8
0.0002
0.0030
0.0291
0.0026
0.0037
0,0162
0.0946
LLM-SR (GPT-3.5)
4.65e-7
0.0005
2.12c-7
3.81e-5
0.0214
0.0264
0.0210
0.0516
Table 1: Quantitative performance comparison of LLM-SR (with GPT-3.5 and Mixtral backbones), and SR baseline models on different scientific benchmark problems measured by Normalized Mcan Squared Error. "N/A" refers to incompatibility of E. coli Growth dataset for the NeSymReS baseline (limited to < 3d data).
asynchronous evolution) (Crammer 2023). and deep taariing-based methods such as DSR (rein-
forcement learning for expression generation) (Petersen ct al. (202I) and uDSR (extending DSR with Genetic Programming search at decoding) (Landajuela et al. (022). We also evaluate against pre-trained Transformer SR models: NeSymReS (Biggio et al. 2021 and E2E (Kamienny et al. 2022). This selection provides a comprehensive evaluation across different SR paradigms. We allow all search-based baselines to run for over 2M iterations until convergence to their best performance. In LLM-SR experiments, each iteration samples b = 4 equation skeletons per prompt with temperature 0.8, optimizes parameters via numpy+BFGS or tormh+Adam (with 30 seconds timeout), and uses k = 2 in-context examples from the experience buffer for refinement. We run LLM-SR variants for around 2.5K iterations in all experiments. More details on the implementation and parameter settings of each baseline as well as implementation specifics of LLM-SR, including experience buffer structure, prompt refinement strategy, and parallel evaluation are provided in App. A
3.3 QUANTITATIVE RESULTS
Accuracy Table 1 compares the performance of LLM-SR (using GPT-3.5 and Mixtral backbones) against state-of of-e--of-the-art symbolic regression methods across various scientific benchmarks. Perfor- mance is measured using Normalized Mean Squared Error (NMSE), with lower values indicating better performance. LLM-SR with both backbones consistently outperform baselines, despite running for fewer iterations (2.5K vs. 2M+ for baselines). To assess generalization capability, we evaluate per- formance on both in- domamn (ID) and out-of-domain (OOD) test sets. The performance gap between LLM-SR and baselines is more pronounced in the OOD setting, suggesting superior generalization of L.M--RRss discovered equations. For instance, on the F. coli growth problem, LLMMSR achieves an OOD NMSE of 0.0037, significantly outperforming other methods (with OOD NMSE 1).
Among baselines, PySR and uDSR show the best performance, while Transformer SR models (NeSymReS, E2E) perform poorly, likely due to limited general- ization from their pretraining on com- mon benchmark distributions to our novel datasets. These results demon- strate LLM-SR's effectiveness in discov- ering accurate and generalizable equa- tions across diverse scientific domains.
Osclation
Derlatinn
home VIE
Itertion
Rerition
Efficiency Fig. 3 shows the perfor- mance trajectories of LLM-SR variants and symbolic regression baselines across different scientific benchmark problems, depicting the best fitting scores achieved over search iterations. By leveraging scientific prior knowledge, LLM-SR ex- plores considerably lower number of equation candidates in the vast optimiza- tion space compared to symbolic regres- sion baselines that lack this knowledge. This is evident from the sharp drops in
Nomalter
Harattoo
Etaration
Figure 3: Best score trajectories of LLM--SR with GPT-3.5 and Mixtral against SR baselines across different benchmark prob- lems. LLM-SR discovers accurate equations more efficiently, requiring fewer iterations. Baselines fail to match LLM-SR even after 2M iterations.
7

---

Published as a conference paper at ICLR 2025
Oscilation
LLMSR (SFT
Felllsee deping fonce
Briing qpsisguraaatl (0)
NeSymRe5
-1i5 Fieee 7. Day7
Belller nettrrig 2014
pare 60
wellieer pereni72)
014 o00 4330003014 E3G771)
wnd -iscog-free depong soreriag.foe p)
Jone
DSR
-' www.opori (h5d) 1110--1222 12009
-3 2rG23) 100 B5p)
uDSR
LLM-SR (Rieba))
decping "p pecaofil ate(s)
Dosping Lare
F- 13 05c
beng
de.sing pucueoti ap sisppecou Eil x paecas (43)
dupping arosdes jenowr
PySR
datring
-E120eperee2" + 632" sange+ car)
-05ti 0500 0.09 X022250)
lef egpntinneenaae, Y, V)
MeSymita5 -0677 2u
Oscilation
LLM-RR
towooar fuve seseruction paeesa iSl rependent an pestiten and velacity diaging Ligh
draving purazalblee sin(t) Persodic esternet [1]
EZE
Jossi aeceterction retern drisiag centoring Interaction dappay parsae [i]
0ze61 01607- E (1 652) (012attat 331))nn5535 oube)
03 +S8LL el) ESP
DSR mot-ariitt
UDSR
LLM-SR (Micral)
Jcover dascig parkcalllve a Proporticnal to tetecits setee aduie paress aucidleeng acmiisl Pasttii interaaiiss Thtst
LIt -aaze 2a -aa taa -ectr---h Uhe ane (9e mea
EC abatt endlring
eeiiserr Juping
PySR
-5BK+ 0.00s - 033 knd 9.1) 1227 9:77l
65 nhre a 2.50000
(a) Ground Truth
(b) LLM-SR Equations
(c) Baseline SR Models
Figure 4: Discovered equations for Oscillation 1 (top) and Oscillation 2 (bottom) problems: (a) True equations and their phase diagram; (b) Equation program skeletons identified by LLM-SR, with simplified forms obtained after parameter optimization; and (c) Equations found using SR baselines. Shaded green terms denote recovered symbolic terms from true equations.
the error curves for LLM-SR variants, indicating they efficiently navigate the search space by ex- ploiting domain knowledge to identify promising candidates more quickly. In contrast, the symbolic regression baselines show much more gradual improvements and fail to match LLM-SR's perfor- mance even after 2M+ iterations. The performance gap between LLM-SR and baselines also mostly widens over iterations, highlighting the effectiveness of L.MSs acting as mutation (or crossover) agents in LLM-SR's ' iterative refinement process.
3.4 QUALITATIVE ANALYSIS
Discovered Equations Fig. presents the final discovered equations for both Oscillation problems using LLM-SR and other symbolic regression baselines. A notable observation is that equations discovered by LLM-SR have better recovered the symbolic terms of the true equations compared to baselines. Also, LLM-SR provides explanations and reasoning steps based on scientific knowledge about the problem. leading to more interpretable terms combined as the final function. For example, in both problems, LLM-SR identifies the equation structure as a combination of driving force, damping force, and restoring force terms, relating them to the problem's 's physical characteristies. In contrast, baselines generate equations lacking interpretability and understanding of the physical meanings of variables and the relations between them. These equations appear operations and variables without a clear connection to the problem's 's underlying physical principles. App. G provides more detailed qualitative analysis of the final discovered equations for other benchmark problems (Figs. 23 and as well as the equations discovered over the performance trajectory of LLM-SR's iterations (Figs. 22p.
combination of mathematical
Generalization Fig. 5 compares predicted distributions obtained from LLM-SR, and competing baselines (PySR and uDSR) with the ground truth distribution of E. coli growth problem. The shaded region and black points indicate in-domain (ID) data, while the rest represent out-of-domain (OOD). Results show that distributions obtained from LLM-SR align well with the ground truth, not only for ID data but also for OOD regions. This alignment demonstrates the better generalizability of equations discovered by LLM-SR to unseen data, likely due to the integration of scientific prior knowledge in the equation discovery process. In contrast, PySR and uDSR tend to overfit the observed data, with significant deviations in OOD regions.
8

---

Published as a conference paper at ICLR 2025
crowe Pass
This overfitting behavior highlights their limited ability to generalize beyond the training data and capture the true physical underlying patterns. Detailed analyses for other benchmark problems provided in App.
Ground Trush
Caeee mar
LLM-SR
Coows Drent Mare
3.5 ABLATION STUDY
We conducted am ablation study on the Oscil- lation 2 problem using GPT-3.5 the LLM backhone to investigate the impact of LLM-SR's key components (Fig- 6 more detailed results in App. Our findings reveal the crucial role of each component in the model's performance. The "w/o Prior" variant, which removes the nat- ural language description of the scientific prob-
Population Density
Subcrate Concentration
Tamperature
pH
Figure 5: Comparison of E. coli growth rate distributions from LLM- SR. PySR, and uDSR.
lem and its variables, led to a considerable performance drop. This highlights the importance of incorporating prior domain knowledge equation discovery. The "w/o Program" variant, which restricts LLM hypothesis generation to single-line mathematical expressions, also had negative but less severe impact on performance, denoting the importance of programming flexibility in this task. The "w/o Iterative Refinement" variant, equivalent to the LLM sampling without the optimization loop, led to substantial performance drops (NMSE: 1.0le-1 in-domain, 1.81e-1 OOD), emphasizing the importance of the evolutionary search and optimization process in LLM-SR's success. The "w/o skeleton + optimizer" variant, which requires end-to-end equation generation without separate parameter optimization step (i.e., generating hypotheses full equations along with their numeric parameters), also significantly worsened results (NMSE: 3.78e-1 in-domain, 3.75e-1 OOD). This highlights the effectiveness of our two-stage approach- opproach-generating equation skeletons followed by data-driven parameter optimization-inn in navigating complex combinatorial optimization space of discrete equation structures and continuous parameters.
We compared two optimization frameworks: numpy+BFGS and torch+Adam. In our experiments, the numpy+BFGS variant performed slightly better compared to torch+Adam. This difference is most likely attributed the LLM's higher proficiency in generating numpy code rather than inherent superiority of the optimization method for this task. LLM-SR relies on direct and differentiable parameter optimization, capability not present in current symbolic regression methods. Combining LLM-SR with LLM backbones that are better in generating Py Torch code could potentially enhance equation discovery by leveraging differentiable parameter optimization in future.
In-drrain
10
Noraalzze MSE
10
10-4
zond
Loces
wib Prior
Wie Program
N/o Iteralive Hefinmmntt
A/o Stoleton Optimiarr
Torch Adam
Wi N Numay+BFGs Ogt mptir
Optimizer
Figure 6: Ablation results on the Oscillation 2 problem, showing the impact of problem specification, program representation, iterative refinement, parameter optimization, and optimization frameworks on LLM-SR's performance.
4 RELATED WORK
LLMs and Optimization While LLMs have shown remarkable capabilities in various domains, their performance often falls short tasks requiring high precision and complex reasoning. To address this. researchers have explored combining LLMs with feedback mechanisms (Madaan et al.,
2024| Yang et al. 20236 Haluptzok et al. 2022) and integrating them into iterative optimization loops
(Lehman al. 2023 Liu et al. 2023
Wu et al. 2024
Lange al. 2024). Recently, LLMs have

---

Published as a conference paper at ICLR 2025
been successfully applied in prompt optimization Yang et al. 2023a; Guo et al. 2024). data-driven analysis (Majumder et al. 2024, Zheng et al., 2023b). and neural architecture search (Chen et al.
2023 Zheng al., (2023a). Most related to work is FunSearch (Romera-Parededr et al. 2024)
that combines L.Ms with systematic evaluators to search for heuristics that push the boundaries in solving some established open mathematical problems. Building upon these ideas, L.MMSRR framework employs LLM as an optimizer, leveraging its scientific prior knowledge and data-driven evaluators discover mathematical equations underlying scientific observations.
LLMS for Scientific Discovery The integration of LLMs into scientific tasks has recently garnered significant attention, offering transformative potential across various fields such as drug discovery, biology, and materials science (Wang et al. 2023a A14Science & Quantum 2023). Specifically, recent studies have demonstrated the capacity of LsMs to propose scientifically plausible and poten- tially novel hypotheses by leveraging their extensive domain knowledge and reasoning capabilities (Majumder et al., 2024 ZZenn et al. 2023b Qi et al.
2023 Ji et al. 2024). Also, when equipped with external tools and scientific simulators, LLM agents have shown promise in automated statistical discovery and reasoning (Li et al. 2024 Wang et al. 2023b; Ma ct al. 2024). Despite the increasing exploration of LLMs in scientific contexts and question answering, their potential for tasks such as equation discovery and symbolic regression remains largely unexplored. Our work extends this line of research by introducing a novel approach for equation discovery that combines LLMs' scientific prior knowledge and code generation with data-driven evaluation.
Symbolic Regression Symbolic regression (SR) methods can be broadly categorized into search- based approaches, learning-based models, as well as hybrid learning and search methods. Search- based approaches mainly explore the space of equation structures and parameters using evolutionary algorithms or reinforcement learning (Schmidt & Lipson, 2009, Cranmer 2023: Petersen et al. [021] Sun al. 2023). They offer interpretable results but often struggle with scalability and efficinccy Learning- based models, on the other hand, leverage large-scale synthetic data and Transformer models to learn the mapping between numeric input observations and output mathematical expressions (Biggio et al. [202 Kamienny et al. 2022). Hybrid methods aim to combine the strengths of both approaches, guiding the search by employing neural priors to improve the expressiveness and efficiency of the discovery process (Landajuela et al., 2022 Shojaee et al., 2024 Mundhenk ct al. 2021 Meidani et al.] (023). Despite the progress made by these approaches, they often face limitations such as the lack of scientific prior knowledge incorporation and the restricted expressiveness of traditional equation representations like expression trees. While there have been some works incorporating prior knowledge by using declarative bias and structures with pre-defined grammars (Todorovski & Dzeroski 1997; Todorovski & Džeroski 2007). these methods do not leverage the power of LLMs for this task. Our work advances this research direction by utilizing LLMs to efficiently search the combinatorial optimization space of equation discovery and generate meaningful equation structures based on the embedded scientific prior knowledge.
CONCLUSION AND FUTURE WORK
In this work, we introduced LLM-SR, a novel approach to equation discovery that leverages the scientific knowledge and code generation capabilities of Large Language Models (LLMs). By treating equations as programs and combining LLM LLM--genrated educated hypotheses with evolutionary search, our method demonstrates superior performance on benchmark problems across diverse scientific domains, particularly in out-of-domain test settings. Despite its promising results, L.MMSRR has limitations. The method's 's performance is inherently tied the quality and breadth of the LLM's training data, which may lead to biases or gaps in certain scientific domains. Additionally, the computational cost of iterative LLM queries and parameter optimization could be prohibitive for large-scale scale problems. Future work could focus on integrating domain-spectie LMs and retrieval- augmented learning techniques to enhance the relevance and accuracy of generated equations; and incorporating human domain experts in the pipeline to improve the scientific plausibility. The creation of more comprehensive benchmarks, designed to simulate true discovery processes and prevent L.M recitation, is also crucial for rigorous evaluation of LLM-based equation discovery methods.
ACKNOWLEDGMENTS
This research was partially supported by the U.S. National Science Foundation (NSF) under Grant No. 2416728.
10

---

Published as a conference paper at ICLR 2025
REFERENCES
B.S. Aakash, JohnPatrick Connors, and Michael D. Shields. Stress-strain data for aluminum 6061- 1651 from 9 lots 6 temperatures under uniaxial and plane strain tension. Data in Brief, 25: 104085, 2019, ISSN 2352-3409. doi: https://doi.or org/ (ttp:/.ib..211.14400000.
Josh Achiam, Steven Adler, Sandhini Agarwal, Lama Ahmad, Ilge Akkaya, Florencia Leoni Aleman, Diogo Almeida, Janko Altenschmidt, Sam Altman, Shyamal Anadkat, et al. Gpt-4 technical report. arXiv preprint arXiv:2303.08774, 2023.
Microsoft Research Al4Science and Microsoft Azure Quantum. The impact of large language models scientific discovery: preliminary study using gpt-4. arXiv preprint arXiv:2311.07.36.1 2023.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. In Marina Meila and Tong Zhang (eds.). Proceedings of the 38th International Conference Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 936-945. PMLR, 18-24 Jul 2021.
Jure Brence, Ljupëo Todorovski, and Sašo Džeroski. Probabilistic grammars for equation discovery. Knowledge-Based Systems, 224: 1070777, 2021. ISSN 0950-7051, doi: https://doi.org/10... 1016//
knosys.2021.107077
Tom B Brown. Language models are few-shot learners. arXiv preprint arXiv:2005.14165 2020.
Angelica Chen, David Dohan, and David So. Evoprompting: Language models for code-level neural architecture search. In A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in Neural Information Processing Systems, volume 36, pp. 7787-7817. Curran Associates, Inc., 2023.
Miles Cranmer. Interpretable machine learning for science with pysr and symbolicregression. jl. arXiv preprint arXiv:2305.01582. 2023.
Roger Fletcher. Practical Methods of Optimization. John Wiley Sons, New York, NY, USA, second edition, 1987.
Qingyan Guo, Rui Wang. Junliang Guo, Bei Li, Kaitao Song, Xu Tan, Guoqing Liu, Jiang Bian, and Yujiu Yang. Connecting large language models with evolutionary algorithms yields powerful prompt optimizers. In The Twelfth International Conference on Learning Representations, 2024.
Patrick Haluptzok, Matthew Bowers, and Adam Tauman Kalai. Language models can teach them- selves program better. arXiv preprint arXiv:2207. 14502, 2022.
Heng Ji, Qingyun Wang, Doug Downey, and Tom Hope. Scimon: Scientific inspiration machines optimized for novelty. In ACL Anthology: Proceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 279-299, University of 2024.
Illinois Urbana- CattannaCCammpaignnCABBBB
Albert Q Jiang, Alexandre Sablayrolles, Antoine Roux, Arthur Mensch, Blanche Savary, Chris Bamford, Devendra Singh Chaplot, Diego de las Casas, Emma Bou Hanna, Florian Bressand, et al. Mixtral of experts. arXiv preprint arXiv:2401 44400804 2024.
Pierre-Alexandre Kamienny, Stéphane d' A'scol,,, Guillaume Lample, and Francois Charton. End- to-end symbolic regression with transformers. In Advances in Neural Information Processing Systems, 2022.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint 142.6690, 2014.
John R. Koza. Genetic programming as means for programming computers by natural selection. Statistics and Computing, 4(2):87-112, Jun 1994, ISSN 1573-1375. doi: 10. 10.1007//F011555 75355..
William La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabricio de Franca, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason Moore. Contemporary symbolic regression methods and their relative performance. In J. Vanschoren and S. Yeung (eds.), Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks, volume 2021.

---

Published as a conference paper at ICLR 2025
Mikel Landajuela, Chak Lee, Jiachen Yang, Ruben Glatt, Claudio P. Santiago, Ignacio Aravena, Terrell N. Mundhenk, Garrett Mulcahy, and Brenden K. Petersen. A unified framework for deep symbolic regression. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems, 2022.
Robert Tjarko Lange, Yingtao Tian, and Yujin Tang. Large language models as evolution strategies. arXiv preprint urXiv:2402. 200::2222..881,, 18381, 2024.
Pat Langley. Data-driven discovery of physical laws. Cognitive Science, 5(1):31-54, 1981.
Joel Lehman, Jonathan Gordon, Shawn Jain, Kamal Ndousse, Cathy Yeh, and Kenneth O Stanley. Evolution through large models. In Handbook of Evolutionary Machine Learning, pp. 331-366. Springer, 2023.
Michael Y Li, Emily B Fox, and Noah D Goodman, Automated statistical model discovery with language models. arXiv preprint arXiv:2402. 10.2..77777 /7879, 2024.
Raymond Li, Loubna Ben Allal, Yangtian Zi. Niklas Muennighoff, Denis Kocetkov, Chenghao Mou, Mare Marone, Christopher Akiki, Jia Li, Jenny Chim, al. Starcoder: may the source be with you! arXiv preprint arXivv2305.06/661 2023.
Tennison Liu, Nicolás Astorga, Nabeel Seedat, and Mihaela der Schaar. Large language models to enhance bayesian optimization. The Twelfth International Conference on Learning Representations, 2023.
Pingchuan Ma, Ssun-Hsuan Wang, Minghao Guo, Zhiqing Sun, Joshua B Tenenbaum, Daniela Rus, Chuang Gan, and Wojeiech Matusik. Lim and simulation as bilevel optimizers: A new paradigm to advance physical scientific discovery. arXiv preprint 27Xiv:240500978: 2024.
Aman Madaan, Niket Tandon, Prakhar Gupta, Skyler Hallinan, Luyu Gao, Sarah Wiegreffe, Uri Alon, Nouha Dziri, Shrimai Prabhumoye, Yiming Yang, et al. Self-refine: Iterative refinement with self-feedback. Advances in Neural Information Processing Systems, 36, 2024.
Bodhisattwa Prasad Majumder, Harshit Surana, Dhruv Agarwal, Sanchaita Hazra, Ashish Sab- harwal, and Peter Clark. Data-driven discovery with large generative models. arXiv preprint arXiv:2402. 13610, 2024.
Michael de la Maza and Bruce Tidor. An analysis of selection procedures with particular attention paid to proportional and boltzmann selection. In Proceedings of the 5th International Conference on Genetic Algorithms, pp. 124-131, San Francisco, CA, USA, 1993. Morgan Kaufmann Publishers Inc. ISBN 1558602992.
Kazem Meidani, Parshin Shojaee, Chandan K Reddy, and Amir Barati Farimani. Snip: Bridging mathematical symbolic and numeric realms with unified pre-training. In The Twelfth International Conference Learning Representations, 2023.
Elliot Meyerson, Mark Nelson, Herbie Bradley, Adam Gaier, Arash Moradi, Amy K Hoover, and Joel Lehman. Language model crossover: Variation through few-shot prompting. arXiv preprint 1rX1v:2302.12170, 2023.
Terrell N. Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio Santiago, Daniel faissol, and Brenden K. Petersen. Symbolic regression via deep reinforcement learning enhanced genetic programming seeding. In A. Beygelzimer, Y. Dauphin, Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, 2021.
Brenden K Petersen, Mikel Landajuela Larma, Terrell N. Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathemat- ical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, 2021.
Biqing Qi, Kaiyan Zhang, Haoxiang Li, Kai Tian, Sihang Zeng, Zhang- Ren Chen, and Bowen Zhou. Large language models are zero shot hypothesis proposers. arXiv preprint 2.l::2311.596, 2023.
12

---

Published as a conference paper at ICLR 2025
Bernardino Romera-Paredes. Mohammadamin Barekatain, Alexander Novikov, Matej Balog, M. Pawan Kumar, Emilien Dupont, Francisco R. Ruiz, Jordan S. Ellenberg, Pengming Wang. Omar Fawzi, Pushmeet Kohli, and Alhussein Fawzi. Mathematical discoveries from program search with large language models. Nature, 625(7995):468-4777 Jan 2024. ISSN 1476-4687. doi:
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. Science Advance, 324(5923):81-335 2009. ISSN 1036-8075. doi: 10.1 126//cience..... 165iieenn.e.1619..... 165893.
Parshin Shojaee, Aneesh Jain, Sindhu Tipimeni, and Chandan K Reddy. Execution based code generation using deep reinforcement learning. arXiv preprint 1rXiv:2301.13816 2023.
Parshin Shojace, Kazem Meidani, Amir Barati Farimani, and Chandan Reddy. Transformer- based planning for symbolic regression. Advances in Neural Information Processing Systems, 36, 2024.
Fangzheng Sun, Yang Liu, Jian-Xun Wang, and Hao Sun. Symbolic physies learner: Discovering governing equations via monte carlo tree search. In The Eleventh International Conference on Learning Representations, 2023.
Ljupco Todorovski and Saso Dzeroski. Declarative bias in equation discovery. In Proceedings of the Fourteenth International Conference on Machine Learning, ICML '97, pp. 376-384, San Francisco, CA,, USA, 1997. Morgan Kaufmann Publishers Inc. ISBN 1558604863.
Ljupčo Todorovski and Saão Džeroski. Integrating Domain Knowledge in Equation Discovery, pp. 69-97. Springer Berlin Heidelberg, Berlin, Heidelberg, 2007. ISBN 978-3-540-73920-3, doi: 01.1007/978-3-54--33222
4.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: A physics-inspired method for symbolic regression. Science Advances, 6(16):eaay2631, 2020. doi: 10. 10.11266scadyv.aay263..
Marco Virgolin and Solon P Pissis. Symbolic regression is NP-hard. Transactions Machine Learning Research, 2022. ISSN 2835-8856.
Hanchen Wang. Tianfan Fu, Yuanqi Du, Wenhao Gao, Kexin Huang, Ziming Liu, Payal Chandak, Shengchao Liu, Peter Van Katwyk, Andreea Deac, Anima Anandkumar, Karianne Bergen, Carla P. Gomes, Shirley Ho. Pushmeet Kohli, Joan Lasenby, Jure Leskovec, Tie- Yan Liu, Arjun Manrai, Debora Marks, Bharath Ramsundar, Le Song. Jimeng Sun, Jian Tang, Petar Veličković. Max Welling, Linfeng Zhang, Connor W. Coley, Yoshua Bengio, and Marinka Zitnik. Scientific discovery in the age of artificial intelligence. Nature, 620(7972):47-60, Aug 2023a. ISSN 1476-4687. doi: 10. (0.1038//41586--23-002222
Ruocheng Wang, Eric Zelikman, Gabriel Poesia, Yewen Pu, Nick Haber, and Noah D Goodman. Hypothesis search: Inductive reasoning with language models. arXiv preprint (300:22009006600) 2023b.
Xingyu Wu, Sheng-hao Wu, Jibin Wu, Liang Feng, and Kay Chen Tan. Evolutionary computation in the era of large language model: Survey and roadmap. arXiv preprint arXiv:2401.1003, 2024.
Zhaofeng Wu, Linlu Qiu, Alexis Ross, Ekin Akyürek, Boyuan Chen, Bailin Wang, Najoung Kim, Jacob Andreas, and Yoon Kim. Reasoning reciting? exploring the capabilities and limitations of language models through counterfactual tasks. arXiv preprint arXiv:2307.02477, 2023.
Chengrun Yang. Xuezhi Wang, Yifeng Lu, Hanxiao Liu, Quoc Le. Denny Zhou, and Xinyun Chen. Large language models optimizers. arXiv preprint arXiv:2309. 230.0030999 2023a.
Kaiyu Yang, Aidan Swope, Alex Gu, Rahul Chalamala, Peiyang Song, Shixing Yu, Saad Godil, Ryan Prenger, and Anima Anandkumar. LeanDojo: Theorem proving with reureeedl-augmeened language models. In Neural Information Processing Systems (NeurIPS), 2023b.
Mingkai Zheng, Xiu Su. Shan You, Fei Wang, Chen Qian, Chang Xu, and Samuel Albanie. Can gpt-4 perform neural architecture search? arXiv preprint arXiv:2304. 10970, 2023a.
Yizhen Zheng, Huan Yee Koh, Jiaxin Ju, Anh TN Nguyen. Lauren T May, Geoffrey I Webb, and Shirui Pan. Large language models for scientific synthesis, inference and explanation. arXiv preprint arXiv:2310.07996. 2023b.
13

---

Published as a conference paper at ICLR 2025
APPENDIX
A BASELINE IMPLEMENTATION DETAILS
A.1 MODELS
We compare LLM-SR against several state-of ofute-of-the-ee art Symbollc Regression (SR) baselines, encom- passing diverse range of methodologies from traditional evolutionary approaches modern deep learning-based techniques. The baselines include:
GPlearn GPlearn is a pioneering and standard genetic programming (GP) SR approach. We use the open-source gplearn" package with the following parameters: Population size: 500, Tournament size: 20, Maximum generations: 2 million. Most of the hyperparameters are set from default setting.
PySR PySR (Cranmer 2023) is an advanced SR method that employs asynchronous multi-island GP-based evolutions. We implement PySR using the open-source pys package with the following settings: Number of populations: 15, Population size: 33, Maximum iterations: 2 million. Except for the number of iterations, other parameters are the same as the default setting in PySR. This configuration leverages the power of parallel evolution over a long time allowing for a diverse and robust search of the equation space.
Deep Symbolic Regression (DSR) DSR (Petersen et al. 2021) employs an RNN-based rein- forcement learning search over symbolic expressions. We implement DSR using the open-source deep-symboll wtt-ssyymllicopptt.. impptimmiin zat ion (DSO) "d package with standard default parameters: Learn- ing rate: 0.0005, Batch size: 512, and Maximum iterations: 2 million. This approach allows for guided search through the space of symbolic expressions, leveraging the power of deep learning to inform the exploration process.
Unified Deep Symbolic Regression (uDSR) uDSR (Landajuela et al.) 2022) extends DSR by incorporating additional linear token and GP search at the decoding stage. We also implement uDSR using the DSC package with the same default parameters as DSR. This unified approach aims to combine the strengths of deep learning and traditional GP methods.
Neural Symbolic Regression that Scales (NeSymReS) NeSymReS (Biggio et al. 2021 is the pioncering pre-trained Transformer SR model for expression skeleton generation. We implement using the Neura a Symbol cRegress Tonoott les repository with the following default parameters: Number of datapoints passed Transformer: 500, and Expression sampling size: 32. It is important to note that this model is limited to pre-training with < 3 variables. Consequently, we only apply the model datasets with dmaz 3, excluding the Bacterial Growth problem (which has 4 variables) for evaluation of this model.
End-to-Ena Symbolic Regression (E2E) E2E (amiieny et al. 2022) is a more recent end-to-end pre-trained Transformer SR approach. We implement it using the symbol cregressior Face- book repository with the following default parameters: Number of datapoints passed to Transformer: 200, and Expression sampling size: 10. This model is also pre-trained for problems with 10 variables.
A.2 DATA PREPROCESSING AND MODEL EXECUTION
For the pre-trained Transformer SR models (NeSymReS and E2E), data normalization is crucial. We apply standard normalization to the input data before feeding it to these Transformer models to ensure optimal performance. For the search-based methods, we allow all baselines (GPleam, PySR, DSR, and uDSR) to run for over 2 million iterations until convergence to their best performance. In the experiments, each baseline undergoes 5 replications. The best results obtained were then documented and reported. This extensive evaluation process, with a large number of iterations and search evaluations, ensures robust assessment of each model's capability to converge towards optimal solutions and effectively explore the vast equation space of each problem.
ckegression
14

---

Published as a conference paper at ICLR 2025
Pind the sathenatinal fanction akcieion that regresenta E, Cols bactariel groath nate, given dsta on population decesty of bactersai apectes, ashstrate concontration, terpenature, and pll leses.
def (data equation):
one
Busluste equation date observations.
one votB, pR, parsas):
fanction bactersal grouth
Arga:
Popslation density
Suestrste concestration
Terpersture
pi: pil Level
paruva: Ilaeeric paranelers to be optinised
Reiura:
Grouth vate
4 xsap te rastial enseple of foectton stmunture
return paraza (0) pa]ana paraas[2] parara [3] parans
ded aguation_vi (B, pl, parans): "n<Isprosed
(a) Bacterial Growth Rate
Find the nathenatical fusction skeiston that represencs strese, gives date on atrusn and Leuparature in Alsnimsur red for both elaxise and plaarii regions nvn
vallute(data,, equation):
Euclucte the aquatioe on data cbservations.
equation vO paranz):
Nathematica: function for strese in Alluninisa rod
Arga:
Stroin level
Tesparaiare
Wuneric parancters optseized
Return:
aseple insiial exarple of faxcison straciore reture params [0] parans[]] [I]+T <T parans [2]
def uquation_vice vi parama): en Ieproved version of equatios s0
Material Behavior Analysis
Figure 7: Example of input prompts program body for (a) E. Coli Growth and (b) Stress-Strain problems, with problem specification, and the initial equation program example (set as simple linear equation skeleton). better readability, the details of evaluation function are not included this figure. Check Fig. for details.
B DETAILS OF LLM-SR METHOD AND IMPLEMENTATION
Hypothesis Generation and Data-driven Evaluation Fig. provided an example of specification for Nonlinear Oscillator problem. Here, Fig. 7 showcases illustrative examples of prompts and specifications tailored for the Bacterial Growth and Stress-Strain problems. These prompts contain descriptions of the problem and relevant variables, expressed in natural language. By providing this context, the language model can leverage its existing domain knowledge about the physical meaning and relations of variables to generate scientifically plausible hypotheses for new equation programs. Fig. also shows a more detailed example of prompt and specification for L.LM-SR that prompts the model generate differentiable equation programs in PyTorch using tensor operations. The prompt suggests using differentiable operators and replacing non-differentiable components (e.g.. if-else conditions) with smooth differentiable approximations.
Our experiments employ either Mixtral-8xx -xrBB (using NVIDIA RTX 8000 GPUs with 48GB memory each) GPT-3.5-turbo (via OpenAl API) the language model backbone. During each prompting step, the language model generates distinct equation program skeletons using a generation temperature of 0.8. This temperature setting is chosen based on preliminary experiments to
15

---

Published as a conference paper at ICLR 2025
Find authenstical fencison form that fitz dala. This feaction form contats and unery aathemstical that differentiable, .xx
inport torch
inport torch.
import torch. optin optin
class Kodel (torch, no Nodule):
def ...nit.... (selt))
super (Nodel, self) .i.t)....ntt_...0
# Isitialsze sodel parsacters
self for torch. (araa Paraneter (torch. tenzor(1. (1.0))))))))))))) (P)])
def forvard(self,, return
torch. torch.Tenser:
Fouard pasa model
equation(*, self params)
def evaluate (data: dict) float:
Rusluate equation seput and ostyut obseruations. aan
# Losd (rue data obzerustions
outputs data('inpute)) data data('outputa')
Defise optceeaer
optinizer optin A 1r=0.001)
optinize equation akeletos paransters
model. train()
nue sue_iterations 10000
(mun for 1 rango stt_tttttt
# Zero the gradvents
optimizer zero.grad()
FFruunn
predicted outputs inputs modei
y-pred model (idgut))
Corpate
torch Dean((y_pzed outpute)
Backaard pasa: compute tosa gradient with nespect to paraneters Loss backward()
# Ujdatt paranclers using op sptinizer
optimizer atep()
Allesrs evaination acore
return -loss. 15ea(() if not (torch. isman(1oss) Borch.isinn (1oss)) any() else
Mone
def equation(x: torch.
paruza: ParaneterList)
Tensor::
ana
Arga:
(tarch, Tenaor): Teput data.
parsns (Sorch. nn, Persnederrisst) List of model suneric constant paraneters,
Returs:
torch. Tenaor: The resuld of applying the authenstical functiion to the
return parama [0] paraaaa [+]
Figure An example of prompt structure, containing problem specification, evaluation and optimization function, and equation program with pytorch tensor operations.
balance creativity (exploration) and adherence to the problem constraints and reliance the prior knowledge (exploitation). To control the length and the complexity of the generated equations and prevent overparameeterization,
set the maximum number of parameters (length of pa pamm vector) as 10 in all experiments. The generated equation skeleton programs are then evaluated to gather feedback. In this framework, deploy e 4 evaluators operate concurrently. This parallelization allows for rapid and efficient assessment of the generated programs per prompt. Evaluation is constrained by time and memory limits set programs that exceed these limits are disqualified and considered as discarded hypotheses by returning None scores. This constraint ensures timely progress and resource efficiency in the search process.
seconds and M 2GB, respectively. Equation
Experience Buffer Management The system stores equation hypotheses and their corresponding data-driven scores in an experience buffer. It uses an islands model with multiple populations (m 10 islands) to maintain diversity. Each island initialized with simple equation, which can be customized for domain-specitf problems. At each iteration, new hypotheses and their fitness scores are added to their originating island if they improve upon the island's 's best score. To maintain the quality and diversity of the experience buffer, we follow [Romera-Paredes, periodically reset the worst- porss-pprforming islands. Every Treset iterations (every 4 hrs), we identify the m islands whose best equation programs have the lowest fitness scores. All the equation programs in these islands are discarded, and each island is reinitialized with a single high- high-performing equation program, obtained by randomly selecting one of the surviving m/2 islands and copying its
al. (0244) and
16

---

Published as a conference paper at ICLR 2025
highest-scoring equation program (favoring older programs in case of ties). This reset mechanism allows the framework to discard stagmant or unproductive regions of the equation program space and focus promising areas. Within each island, further cluster the equation programs based on their signature, which is defined as the equation program score. Equation programs with identical signatures grouped together, forming clusters within each island. This clustering approach helps preserve diversity by ensuring that equation programs with different performance characteristics are maintained in each population.
Experience Sampling To construct informative prompts for the LLM, we sample equation programs from the experience buffer and update the prompt to include new experience demonstration in-context examples. Similar to (Romera-Paredee we randomly select an island from the m available islands. Then, within the selected island, we sample & equation programs (lypically, k 2) to be included as in-context examples in the prompt. When sampling equation programs within an island, we employ a two-step approach. First, we sample cluster based on its evaluation score, favoring clusters with higher scores (i.e., higher-quality equation programs). Let denote the score of the i-th cluster, defined as an aggregation (e.g.. mean) of all the scores in the signature that characterizes that cluster. The probability P. of choosing cluster i is given by:
et al. 2024), here we use two-stage sampling process. First,
exp
= To (1 21 mod N
P.
N
. exp
where Ta is the temperature parameter, is the current number of equation programs in the island, and To = 0.1 and N 10, 000 are hyperparameters. This selection approach known as the Boltzmann selection procedure (Maza & Tidor] 1993). Once cluster is selected, we sample an equation program within that cluster, favoring shorter programs. Let 4 denote the negative length of the i-th program 4,-min, b, within the chosen cluster (measured as the number of characters), and let max,7 1p 10-8- We set
the probability of selecting each equation program proportional to exp (-1,/Tp). where Tp is temperature hyperparameter. The sampled programs are then included in the prompt as in-context experience demonstration, providing the LLM with relevant and diverse examples to guide the generation of new equation programs. By maintaining diverse and high-quality population in the experience buffer and employing a strategic sampling approach, the experience management enables the LLM-SR framework to effectively explore the space of equation programs and iteratively refine its search based on the most promising candidates.
C LIMITATION OF FEYNMAN BENCHMARK PROBLEMS
The Feynman benchmark (Udrescu Tegmark 2020), consisting of 120 fundamental physics problems from the Feynman Lectures on Physics, is widely used to evaluate symbolic regression techniques in scientific equation discovery. However, our investigation indicates that L.LMs have likely memorized many of these well-known physics equations. This memorization poses challenge when using the Feynman benchmark to assess LLM-based equation discovery methods, as may not accurately reflect the models' true discovery capabilities. This section elaborates these limitations and provides evidence supporting the necessity of our newly designed benchmark problems.
To quantify the potential mem- orization of Feynman problems by LLMs, we first conducted comparative perplexity analysis. Fig. iil
Perplexity Analysis
PPL on Moatre
v 9.
lustrates the median perplexity of Feynman problems against our new benchmarks (Oscillation 1, Oscilla- tion 2, and E. Coli Growth) using the Mixtral-8x-7B model as the LLM backbone. Perplexity, in this con- text, is calculated only for the generation of equations given scientific context: p( ples of input prompts and outputs used for perplex- ity computation across different benchmarks are pro-
Feyrman Problems
Ciclation
Decilation?
E.col Growth
Exam-
Figure Perplexity (Mixtral) comparison of Feynman benchmark and our new designed bench- mark problems
vided in Fig. 10 Mathematically, perplexity is defined as: PPL exp
IN Ler,loogee logp(zz/ze))
where N is the number of tokens in the generated equation, and p( (x<s)) is the probability of token
17

---

Published as a conference paper at ICLR 2025
Feynman Problem Examples
New Benchmark Problems
Sceentoe Dontee
Equarion
Scirtiis Connett
Fryustion
Reuatiin of grasttational fores (P) hetuens Twn taxees et) reperater by diatanss grasitationa)
Spartiee of outhee doatlaso eccitlatin for acoddrretice (eoie)) aoeen podiiise (e),
dusee
(32-93322 (22-23)72)
and
P("quation Contoot) 9 PPL- 1.24
Teynman L.. Problem
JCFqquation Cooneett PPL 6.24
Oscllation Problem
Scardite Contic
Equation
Sclentifi: Cornot
Equation
Equation wownten pertickle excs (a), uslactry (O,, dietuees
aze(thrta)
Eqartise ronuimean excazlatar for accelorarics (drid)) goaen (1), poctitionn tol,
betwein and
EEEuatioon Contex)) PPL:2.44
Peyrmmnt 10.1d Prellen
#(Equation Contert) PPL: 6.27
Cacilaation 2 Frntlen
Scaentie Contee
Equation
Sejertits: Connest
Frgustion
feuation relating evngy (8), mass (m), speed (cl,
<2)
(12 sereC vi
single pepsiation E. taus proulacize rice (18/012 popelaaiio (8), sabetrate concetratiin 153, tenganature
oettt side
toent astr ttecttra (acagghn) ie alagtel lpe wocees alnjer
M(quuation (Comlot)) 2.0-
Foynman 148.2 2 Probiom
((liquation [Conbeet) PPL- 6.50
Coll Growth Frobtem
Figure 10: Examples of input context and output equations for Feynman equations and our new benchmark problems in Perplexity experiments with Mixtral LLM backbone.
given the preceding tokens, derived from the logits of the Mixtral model. The significantly lower perplexity observed for Feynman problems indicates a higher certainty the LLM's predictions for these equations, and higher chance of LLM recitation rather than reasoning and discovery. This suggests a high likelihood of memorization of well-known Feynman equations by the LLM, potentially due their prevalence in scientific training data. It is worth noting that have excluded the Stress-Strain problem from this analysis due to its experimental nature and lack of a predetermined theoretical model structure, precluding the calculation of perplexity.
Discovery Error Curve Analysis To further validate the need for new benchmarks, we com- pared equation discovery error curves between Feynman problems and our new benchmark problems. Fig. II presents the performance of LLM-SR with GPT-3.5 backbone across var- ious problems, showing the best score trajectory of Normalized Mean Squared Error (NMSE) against the number of iterations. For Feynman benchmark problems, LLM-SR achieves low NMSE scores within very few iterations, often in single pass. This rapid convergence further supports the hypothesis that I.LMS have likely memorized these fundamental physics equations due to their ubiquity in training data. Qualitative examples in Figs. 12 and 13 also provide addi- tional evidence. The LLM's one-pass responses to several Feynman problems not only demon- strate functional accuracy but also often recite the exact form of the corresponding physics pressions, suggesting direct recall rather than discovery process. In contrast, our newly designed benchmark problems present novel challenges requiring reasoning and exploration.
10a
Normalized MSE
10t
Oscillation
Oscillation 2
Growth
10-*
Stress Stras-
Feynman
200
400
600
800
1000
1200
1400
Iterations
Figure 11: Trajectory of Normalized MSE score over iterations for LLM-SR (GPT-3. 3.5) 5) Feynman bench- mark problems versus new benchmarks
D ADDITIONAL DETAILS ON NEW BENCHMARK PROBLEMS
The datasets used this study include both publicly available and newly generated data. The material stress behavior analysis dataset (stress-strain) is publicly available under the CC BY 4.0 license and can be accessed at
The remaining datasets (Oscillation 1, Oscillation 2, and E. coli Growth) generated for this work and released under the MIT License as part of the LLM-SR SR GitHub repository: https::/gituub.cm//eep.s.ym mmahemmatie//IIMMS5P
18

---

Published as a conference paper at ICLR 2025
Input Prompt
Find the natheratical function skeleton that represents the ressldant snfenssty due isterference cokenest light
mns
def (data equation):
Pvaluate equation observations.
vo 12, delta, parans):
Nathonstical function for the ressitant intenesty due interference of 500 woves
Arge: Intenaslyy firet uave
Interesty of second uene
delta: Phase difference betuees tue aases
paraas: Paneric parsneters be opiswized
Aeturn:
Resuitant intesesty
steple initial exoeple of function structure return paranr [0]+11 parame[l3+12 parans[2]sdalts parase [i]
def equation
parans):
siooopooeed verasos of equation
LLM Response
parans parana
(21+(()) cos(delta)
parans [i]-I1 [2]
Feynman 1.37.4 Problem
1 + 12 + 2, 12 co)
Figure 12: example of LLM response Feynman 1.37.4 problem, demonstrating LLM recitation without iterative search. For better readability, the details evaluation function are not included this figure. Check Fig. 2 for details.
D.1 NONLINEAR OSCILLATOR EQUATIONS
In this work, we introduce two novel nonlinear oscillator systems as part of our effort to create more challenging and robust benchmarks for LL.M-based equation discovery methods. These custom- designed oscillators extend beyond commonly studied systems like Van der Pol, Rayleigh, or Duffing oscillators, presenting unique challenges that test the reasoning and discovery capabilities of LLM- based approaches.
The general form of nonlinear damped oscillator equations is typically expressed as: ž + f(t, 21 2) O, where represents time, x represents position, and ftt, x, æ represents nonlinear forces. Our custom designs expand upon this framework, incorporating rich combination of nonlinear terms to create systems that are challenging yet solvable. We simulate two nonlinear oscillators using the sol volve_ie ve_ivp function from the scipy library generate data.
Dataset
Time range
initial values
F
3
7
Oscillator
(x=0.5, v=0.5)
0.8
0.2
0.5
Oscillator 2
(0, 50)
[x=0.5, v=0.5)
0.3
0.5
1.0
5.0
0.5
1.0
Table 2: Parameter values Oscillator datasets.
The parameters and initial values for these simulations provided in Table equations for our oscillator systems are as follows:
The governing
Oscillator 1:
5 F Fsi((wx) ox0 - Bxx - -xx - x cos(2r)
Oscillator 2:
0 Fii((wt) an Br v dx exp(71c)
19

---

Published as a conference paper at ICLR 2025
Input Prompt
Find outheracical Junction skeletos representa relatsonsksp angle of Kncidence asgle of refraction.
def evaluste equation)::
von
Evslunts the equalion on daís obasroutions.
def equation, vō (n, theta, parane): #atheratical /unction for angls of refraction
Arge:
Rafractive inder
thete: Asgle of sneidenes
panaxs; Muterie paraxeters to be optsnized
Return:
Refraction ongle v.e steple
exapple of fancéson stracture
return paraen[0] paramali] +theta params [2]
def equalion, (u, Limin, puzsae):: ennTogrooud versies equatson
LLM Response
return parame [o]
np. arcsia
up appeii(theta)) parane
Feynman 1.26.2 Problem
0, = aresin( (0))
Figure 13: An example of LLM response Feynman 1.26.2 problem, demonstrating LLM recitation without iterative search. For better readability, the details evaluation function are not included this figure. Check Fig. for for details.
where v = å represents velocity, and F, wi, ar, 33 are constants specific to each oscillator system. These equations carefully deisgned to incorporate diverse set of nonlinear structures, including trigonometric, polynomial, and exponential terms. This design choice serves multiple purposes:
Challenging Complexity: The combination of various nonlinear terms creates rich dynamical system that is more complex than common oscillator systems, making the equation discovery task non-trivial.
While complex, these equations still solvable and still represent physically plausible systems, incorporating recognizable clements such as nonlinear damping and position- dependent restoring forces.
Realistic Physics:
Novelty: By deviating from well-known oscillator forms, we reduce the likelihood of LLMs simply reciting memorized equations, thus testing their reasoning and discovery capabilities in the context of data-driven scientific equation discovery.
Fig. 14 Illustrates the phase plane diagrams of these nonlinear damped oscillators, visually demonstrating the complex dynamics arising from the interplay of nonlinear driving forces, restoring forces, and damp- ing forces. These diagrams highlight the rich behav- ior that makes these systems challenging for equa- tion identification tasks. To effectively evaluate the generalization capability of predicted equations, we employ strategic data partitioning scheme. The simulation data is divided into three sets based on the trajectory time: (1) Training set, (2) In-domain validation set, and (3) Oft-oo--dmaan validation set. Specifically, utilize the time interval T [0, 20) to evaluate the out-of--omain generalization of the discovered equations.
Diphoormet BC
Figure 14: Phase diagrams of trajectories corre- sponding to custom oscillators: (a) Oscillator 1 and Oscillator
20

---

Published as a conference paper at ICLR 2025
D.2 E. COLI GROWTH RATE EQUATIONS
In the domain of microbiology, understanding and modeling the growth dynamics of Escherichia coli (E. coli) is of paramount importance due its wide-ranging applications in biotechnology, food safety, and fundamental biological research. To advance the LLM-based equation discovery approaches in this field, we have developed a novel benchmark problem centered around E. coli growth rate modeling. The growth rate of bacterial populations, including E. coli, is typically modeled by a differential equation that incorporates multiple environmental factors. This mathematical model commonly takes a multiplicative form:
dB dt = (T. fz()) fs(S) fT( (pH(
where represents bacterial population density, is substrate concentration, i' temperature, and pH represents the acidity or alkalinity of the growth medium. To create a benchmark that is grounded in biological prior knowledge yet is challenging, we have extended this framework with custom differential equation:
dB MmmonBB
S
tanh (k(T zul)
(pH pH min pHmmar pHmmin
I c(T' Edewayy ) exp (-IpHI-
-pl
sin
dt
K, K.+S7
This equation incorporates several key components, each designed to test different aspects of equation discovery systems:
Population Density (fB): We maintain a linear relationship with B, reflecting a simple single- population scenario. This choice allows the focus to remain on the more complex environmental dependencies.
Substrate Concentration (fs): We employ the well-established Monod equation, K,+S
which
has been comnerstone of bacterial growth modeling since its introduction by Jacques Monod in 1949. This inclusion serves as prior knowledge, allowing us to evaluate how well discovery methods can identify known relationships within more complex overall structure.
Temperature Dependency (fr): We introduce novel formulation, tanh(k(TT=oo)
-ddeayy which cap- tures the non-monotonic response of bacterial growth to temperature changes. This function combines a hyperbolic tangent term, representing the initial growth acceleration with temperature, and quartic decay term, modeling the rapid decline in growth rate at high temperatures. This formulation presents new and challenging form for LLLM-based equation discovery methods, introduces operators and structures not commonly seen in the literature of this scientific context.
pH Dependency (fa)): Our custom pH function, exp (- |pH pHopel) sin (pIl-pll..aa )m pil -pl combines exponential and trigonometric terms to model the complex relationship between bacterial growth and pH levels. This formulation captures both the optimal pH range for growth and the symmetric decline in growth rate as pH deviates from the optimum. It also poses challenging setting for LLM-based equation discovery methods with structures uncommon in the relevant scientific context.
Fig. 15 illustrates the behavior of our custom-designed (T') and fpt (pH) functions in comparison established models from the literature for temperature and pH impact in bacterial growth. As observed, our custom models maintain key characteristics of bacterial growth responses while introducing complexities that challenge equation discovery methods. The temperature dependency model shows a sharper optimal peak and more rapid decline at high temperatures compared to traditional models, while the pH dependency model exhibits a narrower optimal range with steeper declines outside this range.
This carefully constructed benchmark serves multiple purposes: (i) It leverages LLMs' prior knowl- edge of bacterial growth patterns and common mathematical functions used in biological modeling; (ii) It prevents trivial LLM recitation by introducing novel combinations of functions and opera- tors that go beyond standard models; and (ii) It challenges equation discovery systems to identify complex, biologically plausible relationships from data, simulating the process of scientific discovery.
21

---

Published as a conference paper at ICLR 2025
Cardivl Ierperstum Model kattouetly Model
Caidioll Mecel
Zaietring
Gaustion
Seponentia Cuttor Hodel
Cutto Mecel
10
00
Temperature (TC)
2H
Figure 15: Scheme of some established models from literature for temperature and pH impact in bacterial growth compared to custom-destt dessgndd model behavior.
Temperature
20
250
100
150ºC
Stress (MPa)
200
250
300
150
100
0
0.00
0.05
0.10
0.15
Strain
Figure 16: Stress-strain curves of Aluminium 6061-T651 under various temperatures (data from (Aakash et al. 2019
D.3 MATERIAL STRESS BEHAVIOR ANALYSIS
The analysis of material stress behavior with regard to temperature variations is critical area of study in materials science and engineering. Our focus experimental data of Aluminium 6061-T651 (Aakash et al. 2019) provides a new experimental case study for LLM-based equation discovery methods. Fig. 16 presents stress-strain curves for this alloy across range of temperatures, offering rich insights into its mechanical behavior. For this problem, the experimental data represents the tensile behavior of material under uniaxial tension for 6 different temperatures from 20 to 300 We allocate the data corresponding to T C00 for use as the out-of-domain validation set.
The stress-strain curves in Fig. 16 reveal several key features: Temperature Dependence: As temperature increases from 20 to 300CC, we observe a significant decrease in both yield strength and ultimate tensile strength. Elastic Region: The initial linear portion of each curve represents the clastic region, where deformation is reversible. The slope of this region, known as Young's modulus, appears decrease with increasing temperature, indicating reduced stiffness higher temperatures. Plastic Region: Beyond the yield point, the curves exhibit non-linear behavior characteristic of plastic deformation. The shape of this region varies with temperature, suggesting changes work hardening behavior. Failure Region: The endpoints of the curves indicate material failure. Notably, the strain at failure generally increases with temperature, implying enhanced ductility at higher temperatures. Complex Non-linearity: The overall shape of the curves exhibit piece-wise form, particularly in the plastic region, displays complex non-linear behavior that varies significantly with temperature and cannot be simply modeled with closed-form mathematical expressions.
The complex, nonlinear behavior of materials under varying experimental conditions poses significant challenges in developing comprehensive theoretical models. Moreover, unlike many physics problems where equation forms might be known or suspected, stress-strain- tures--strin-temppeauuuu relationships for specific materials often lack universally accepted theoretical model. This absence necessitates
22

---

Published as a conference paper at ICLR 2025
12
12-1
04
443-00
3 Rt
3 Pl
10
Py5R
LUM SA
LUM 5A
LU
PySR
LLMSSR
LLM SA Iw/
SM
(a) Oscillation 1
(b) Osciliation
domaii
000
100
12
Narmaland
PySn
LUX
Pysk
LLM--A tnto Program)
LLWAA
(wto Prior)
(wi Program)
Priar)
E. Coli Growth
(d) Streest-traat
Figure 17: Ablation comparison of LLM-SR variants (GPT- 3.5 backhone) and PySR on four henchmarks, showing in-domain and out-of-domatte errors for (a) Oscillation 1, (b) Oscillation 2, (c) E.coli Growth, and (d) Stress-Strain problems. LLM-SR consistently outperforms PySR and ablated versions, highlighting the importance of both prior knowledge and equation program representation components.
more exploratory, data-driven approach empirical equation discovery. By incorporating this real-world experimental materials science problem into our benchmark problems, we aim to evaluate the ability of LLMM-RR to discover physically meaningful and interpretable equations in a domain where empirical modeling often dominates
E ADDITIONAL ABLATION STUDY
To evaluate the key components of our LLM-SR framework, we conducted additional ablation studies across all benchmark datasets. These experiments specifically target two main motivations behind LLM-SR: (1) leveraging scientific prior knowledge and (2) utilizing code generation capabilities of LLMs the context of scientific equation discovery. We present the results of these ablation studies in Fig. 17 which compares the performance of three LLM-SR variant models: (i) LLM-SR (w/o Prior): This variant removes the incorporation of woolem--specife prior knowledge from the LLMs' input prompt; (ii) LLM-SR (w/o Program): This variant eliminates the use of program representation for hypotheses; and (iti) LLM-SR: This variant is the final version of LLM-SK model including all components. We evaluated these models on all benchmark problems. For each problem, we assessed performance under both in-domain and out-of-domain (OOD) test settings, providing a robust evaluation of generalization capabilities in discovered equations.
Fig. 17 shows that both ablated components are critical to LLM-SR's success. Across all datasets and in both test settings, w/o Prior and w/o Program consistently underperformed compared the full LLM-SR model. In both Oscillation 1 and Oscillation 2 problems, prior knowledge and code generation capabilities demonstrated comparable impacts on model performance. This held true for both in-domain and OOD settings, suggesting that both components contribute similarly to the model's understanding and generalization oscillatory systems. For the E. coli Growth problem, we observed a more pronounced effect of prior knowledge. This indicates that domain-specific knowledge plays a particularly crucial role in modeling for this problem. Interestingly, while both components had similar impacts in in-domain performance in stress-stras strain problem, prior knowledge appeared to have a more substantial effect on OOD performance. This aligns with our intuition that OOD generalization should correlate more strongly with prior domain knowledge of the system.
23

---

Published as a conference paper at ICLR 2025
Stress Strain
Oscillation 2
Model
ID4
00D,
IDL
OODJ
LLM-SR (wlo multi-island & sampling)
0.0257
0.1010
6.23e-6
0.0008
L.M.SRR
0.0210
0.0516
2.12e-7
3.81e-5
Table 3: Impact of multi-island design on LLM-SR performance (with GPT-3.5 backbone) across two benchmark
problems measured by Normalized Mean Squared Error.
We further conducted ablation experiments to analyze the impact of the multi-island buffer design and corresponding sampling strategy. Table compares the performance of LLM-SR with and without the multi-island components on two benchmark problems. Specifically, we evaluated variant that uses only one island in the buffer and employs simple deterministic top-k selection approach for in-context example selection. The results demonstrate that the multi-island design positively impacts LLM-SR performance in both in-domain and OOD settings. Qualitatively, we observe that the number islands plays critical role in balancing exploitation and exploration. With fewer islands, the framework exhibits reduced exploration capabilities, generating less diverse equation hypotheses and converging prematurely to equation structures produced in early iterations.
F ADDITIONAL EXPERIMENTS
F. NOISE ROBUSTNESS ANALYSIS
To evaluate the robustness of LLM-SR imperfect data conditions commonly encountered in real- world settings, we conducted a systematic analysis model performance under varying levels of noise. We specifically focused on the Oseillation benchmark problem, introducing controlled Gaussian noise with different standard deviations (a {0, 0.01, 0.05, 0.1)) to the training data. This analysis provides insights into how the incorporation of prior knowledge through LLM-SR affects equation discovery performance when dealing with noisy measurements.
(ID)
1.4
LLM-SR (COD)
PySR (ID)
1.2
(OOD)
1.0
NMSE
0.2
0.0
0.01 Noise Level (a)
0.05
Figure 18: Noise robustness analysis of LLM-SR compared PySR on the Oscillation benchmark. The plot shows normalized mcan squared error (NMSE) for both in-domain and oft-of-domain (OOD) predictions under different levels of Gaussian noise (a = 0, 0.01, 0.05, 0.1). While performance degrades with increasing noise for both methods, LLM-SR maintains better robustness, particularly for OOD predictions, demonstrating the value of incorporating domain knowledge in noisy real-world settings.
Fig. 18 presents a comparative analysis between LLM-SR and PySR under different noise conditions. The results demonstrate that while increasing noise levels generally degrade performance across all methods, LLM-SR exhibits notably better resilience to noise compared to traditional approaches. Specifically, moderate noise levels (a 0.01, 0.05), LLM-SR maintains significantly lower NMSE compared to PySR, particularly in out-of -o-maan predictions.
This enhanced robustness can be attributed the incorporation of domain knowledge through LLM prompting, which helps constrain the search space to physically plausible solutions even in the presence of noise. These findings highlight an important aspect of symbolic regression: data
24

---

Published as a conference paper at ICLR 2025
quality decreases, the value of incorporating prior knowledge becomes increasingly significant. This observation aligns with the broader principle that when evidence (data) becomes less reliable, the role of priors in inference becomes more crucial.
F.2 LLM-BASED OPTIMIZATION BASELINES
Several frameworks have emerged recently exploring the integration of LLMs into optimization tasks (Meyerson al. 2023 Yang et al. 2023a Romera-Paredes al. 2024). Among these works, LMX (Meyerson et al. 2023) has included symbolic regression as one of its experimental tasks, though primarily as a proof-of-concept rather than aiming achieve state-of-the-ar performance. Their implementation does not incorporate domain-specifer prior knowledge, generates complete equations with L.M instead of optimizable skeletons, uses older LLM models (Galactica, Pythia), and lacks several design elements present in LLM-SR such equation-as-pre ps-prgram representation and multi-island dynamic memory management for diverse exploration. We conducted experiments to evaluate LLM-SR alongside LMX and FunSearch Paorees-Paredes et al. [024). For LMX, we utilized their open- implementation directly, while for FunSearch, we adapted their prompt and feedback design to suit the equation-as as-prgram task. All methods were run for the same number of iterations as LLM-SR, using GPT-3.5 as the LLM backbone.
Oscillation 1
Oscillation 2
Model
IDI
OOD1
IDI
OOD|
DSR
0.0087
0.2454
0.0580
0.1945
uDSR
0.0003
0.0007
0.0032
0.0015
PySR
0.0009
0.3106
0.0002
0.0098
LMX
0,5031
0.9371
FunSearch
0,4840
8.059
0.7234
0.5861
LLM-SR (wla skeletoo--pttimzee))
0,1371
0.6764
0.3780
0.3751
LLM-SR (w/o Prior)
0.0001
0.0029
4.65c-5
0.0071
L.M-SRR (w/o Program)
0,0001
0.0035
3.22e-5
0.0043
L.LM-SR
4.65e-7
0.0005
2.12e-7
3.81e-5
Table 4: Additional experimental results comparing different optimization methods on oscillator problems measured by Normalized Mean Squared Error.
Results in Table show that LLM-SR and its variants achieve better performance than LMX and FunSearch in both in-domain and OOD test settings. Even when the skeleton+optimizer design is ablated (allowing LLM to generate complete equations without placeholder parameters), performance drops significantly but still outperforms both LMX and FunSearch. Further ablation of prior knowl- edge or program representation components shows less dramatic impact, though the full LLM-SR incorporating all components achieves the best performance overall. Notably, traditional SR baselines (DSR, uDSR, PySR) also outperform the existing LLM-based approaches on both datasets, supporting our focus on comparing against ttate-of-the-an SR methods in the main experiments.
F.3 LLM BACKBONES
We have also conducted experiments with larger LLMs as the backbone for LLM-SR (than current Mixtral and GPT-3.5 backbones). Due the computational budget limitations, we only conducted experiments on the oscillation 1 and oscillation 2 datasets. For fair comparison, we ran this GPT-4o backbone model for the same number of iterations (2500) as other LLM-SR runs. Table 5 shows that on both datasets, LLM-SR (with GPT-4o) provides slightly better performance than LLM-SR with GPT-3.5 and Mixtral backbones, particularly in the OOD settings (reported numbers NMSE in Table D These findings align with expectations that improved LLMs with better knowledge, reasoning, and programming capabilities have potential to enhance performance in LLM-SR framework.
Oscillation
Oscillation 2
Model
IDI
OOD
IDI
0OD]
LLM-SR (Mixtral)
7.89c-8
0.0002
0.0030
0.0291
LLM-SR GG.T3355
4.65e-7
0.0005
2.12--7
3.81e-5
LLM-SR (GPT 40)
7.29c-6
6.75c-8
4,276-10
1.29-66
Table 5: Additional experimental results comparing different LLM backbones oscillator problems measured by Normalized Mean Squared Error.
25

---

Published as a conference paper at ICLR 2025
Equation: 0.3 sin()) -0.5x* -x 5x.exp(0.55
(13 5.0 e²)
001/5
0.3 sin(t) 5.0x 0.03v 0.06'' 2.3x²
0.1100
Normalized MSE
0.0125
0.3 5.0x 0.03v 0.07v 2.3x 0.25v/vl
(t3 + 5.0(1 - e²) + 0.02v 0.23v/o|
|
(.3 5.0(1
0.0255
0.0000
Iterations
Figure 19: Performance Trajectory of LLM-SR (GPT-3. along with the best-scoring simplified equations (after parameter optimization) over iterations on the Oscillation problem. Green highlights indicate recovered symbolic terms from the true equation.
G ADDITIONAL QUALITATIVE RESULTS
Discovery Trajectory In this section, we evaluate the progress of generated equations using LLM-SR over iterations. This analysis can illustrate the qualitative evolutionary refinement of the discovered equations.
Fig. shows the NMSE values for the Oscillation dataset. For simplicity, we have provided the simplified equation versions of programs with their optimized parameters. We observe that some of the common nonlinear terms such as sinisoidal driving force term are found early in the search, while more complicated nonlinear terms are found later in the search. An interesting observation here is that LLM-SR + 5(1 exp(x)) at the end. By evaluating
while ground truth equation for this dataset is = 0.3 sin(t) 0.5c x-U - 5x dar-exp(0.522
has discovered the equation ú 0.3 sin(t) - 0.5c the different terms in these two forms, we observe that in fact 5(1 exp(x)) )) R -5x exp(0.5x) for (-2, 2), which is the approximate range of displacement in this dataset.
Fig. 20 illustrates the evolution of equation program skeletons for the Oscillation problem. It can be observed that the model attempts incorporate various nonlinear terms corresponding to driving, restoring, and damping forces, as evidenced by comments or variable names within the code, aiming enhance accuracy.
Similarly, Fig. 21 presents an annotated performance curve illustrating LLM-SR's performance on the E. coli growth rate equation discovery benchmark problem. It becomes apparent that the model recognizes the potential presence of optimal values for temperature and pH from the early iterations which comes from model prior knowledge about the common bell-shaped effect of these variables on the growth rate (check Fig. 15) To enhance accuracy. the model necessitates exploration and incorporation various nonlinear forms. Notably, LLM-SR directs its evolutionary changes towards the more critical and variable aspects of the problem, specifically the pH and temperature effects, as opposed to other components such as substrate concentration represented by the Monod equation K+S. Additionally, the figure demonstrates LLM-SR' comprehension that different components of the function should be multiplied together in the final step, underscoring how prior domain knowledge of the problem structure can guide LLM-SR's evolutionary steps.
Fig. 22 displays three distinct equation skeleton programs discovered by LLM-SR for the stress-strain problem over search iterations. As in previous cases, we notice the model's enhancement through exploration and incorporation of various nonlinear terms into the equations, particularly in the plastic region terms. An additional significant observation for this problem that stress-strain relationships often exhibit piece-wise behavior (as it can also be observed in Fig. 16), which fosed-form expression
26

---

Published as a conference paper at ICLR 2025
def equation (paraeo v9:
def eguation (parses, X, 93: Deswing force
leattorng ferce
restaring, "paraee[]]
driving. force paramelo] dinpppranee x)
dhaping, forze [parae([))
Peteckty tene linezr relocity paraze (21
driving_foree paraas [2] C ge+2) sp.siga(o)
chrviing joret
Daeping Joree
nealinear daxping parans (pppaass)))
nonitneer restoring
parama ges3
Popitson-eeleefp interactson
interactios tera "patemu[5]
Bue up all feneo lo zet lotal force tosel, forse raatoring.freee fesping,foeee friving,.orea mentttas..eetoring
8ax all ferms to the totel ferce
forel driving force linear, velocity zonlinoor, damping interaction
retura totui_tores
raturn total, ForcH
equation (garans, Posston-usiocsty enteraction demping nonlinear, damping -pareas vesporans(])
Veloeilp darpiny
Normalizd MIEE
aolocitydduppin "paasua verparamm
Drining fonce driving.forc parame up ap.ein(parann [6] x parans(e))
8ox vp all terns to get total force
30"4
tocal, force monlineax dasping relocisy.duageee driving_force
zuburn cotal_ Force
Iterations
Duaginn daxping "parans
parans
zp. abp.abs(v)
8 Brsuseg Sern driving paranas E20 ap man (paraan [S3 parans
Sun all terms ged total, force driving
Figure 20: Performance Trajectory of LLM-SR (Mixtral) along with the best-scoring equation program skeletons (before parameter optimization) iterations the Oscillation I problem.
structures traditional symbolic regression models struggle to model. However, LLM-SR represents equation skeletons as programs, thus, it can employ conditional rules (If-Else) or their continuous relaxations, utilizing step-wise nonlinear differentiable functions such as the s si sigmoid d function model smooth piece-wise relations. This diffcrentiability and smooth approximation of if-else conditions are particularly helpful for the parameter optimization step, providing smooth functions for the optimizer to navigate.
Discovered Equations Fig. 23 and Fig. 24 depict the equation programs identified by LLM-SR and other leading symbolic regression baselines (DSR, uDSR, and PySR) for the E. coli growth and the stress-strain problems, respectively. The diverse range of equation forms identified by different symbolic regression methods reflects the challenges posed by these datasets. Notably, in both datasets, the SR methods yield either lengthy highly nonlinear equations that are not aligned with the prior knowledge of the systems, as evidenced by their poor out-of-domain (OOD) performance scores in Table I In contrast, LLM-SR finds flexible equation programs that are more interpretable and aligned with the domain -smaiii scientific priors of the systems.
Behavior Discovered Models Fig. 25 and Fig. 26 offer qualitative comparison by visually presenting the outputs of the equations obtained using LLM-SR, PySR, and uDSR. Upon examination. it becomes evident that the predictions generated by LLM-SR SR exhibit a notable degree alignment with the ground truth data. This alignment suggests that LLMSSR effectively captures the underlying patterns and dynamics the data.
27

---

Published as a conference paper at ICLR 2025
(nf aqurtinngarass, 3, 1, 24)
opactio./peres #0: parsaafol Sptiest Eaguratuure Gatatas
1,aps pareza Opptenal teperstore for Oxtetes 28.9po
Popslatton desatty affect
ww..gnwhhurete De
Souttrste onsesiration effect (beee)
Jopslation offect (practh unitection high deurttia)) peplletio_eefees up powett, peress DJ)
Treprrstune degradenc
effect parmaat)) T.eps)
dskscracs somrstration affest (Rowdd) wdstreee.etfees
pecuzn [4]
pul degredmnn
edfect sp.erpt-pee paccaaall istt eot) Nuriese grast refe
Zeperaterr depmiiascy (iasssii))
Tapi)rs))
sar.powtsteaee paraan]))
sffect parmm ipt pkospc)=e))
Gabiise all offrete io greath vafe sempecetoreaattes gl.effect pecmma(7))
Gayiste currylng esperityy sflees II -1-- B peccas (1)
tone)
Coavioe ali offerts get the griwch wato gront (ees growc cote pepelesiie eddect rebotente edfect currving stfect)
def reuation B, fa, Ra, M. paraas
Hil:
Normalized
resperttrre
estctation fester
ap.etpigenno,
ops3))
Pepuistion desesly offect pepilasionaatoon up. powee algha)
Subuerate consentration offect (losed) subtteute eifees
on5
teparasuraaafoor sp-arpl-Rt up-aqarast Taap03
Pespersture dependency
det scoatize/ppereta, 3, T. p00:
Z.ggi. purecall Cpiemut coperstture Getarto
dapendrazy
PRaefect ap.espl-ts to ap.aquaze(pt plLopt))
iibutroie coneeirutioo offect Sood (orstet)]
R03
droth.rots 86
repprrttree desendercy (arrhentus) ap. 30 (T
pepulation, affacs sfserata.affeett
sesperature affect pl,affect desth tecpffactor pat Loclor
effect fespeemncca decas frem apties! pe)
recure grte
Papulaties dessity Loniaatien (lagiates grouch
Cubiss offects get groath
retura poots_race
Figure 21: Performance Trajectory of LLM-SR (GFT-3.5) along with the best-scoring equation program
skeletons (before parameter optimization) iterations the E. coli growth problem.
28

---

Published as a conference paper at ICLR 2025
9d tuesinnlpermm, 29:
6et (eaatiospparacs, 6, Ti:
depesdeace
Sepack the porseetare
initial, paraza (0)
alph.. tesa, e0, nige.y.. alga.0, q. parazo
Prasattion fector betoees clasctc and plaatic regsons
Lineer clastsc *espsnar retpor alpan
Alpto 1 (Eeporaca Stpp (ppanee (5)-0-ppccoo00)))
Cinear chaiis slaesic_tes term [T prraa][]+*
plaatic +iges bets / [t wp.oddrdeite (Gs 60)
(ap.capl 1 21100>20)
Amiisonn plasinc lew plassic cerm paraza [9]
reptons
tranestion (
+))
deciiso alastee and plactto recpesios attres sotesad attere slpha) (elartic, plantac.vaaahhiieee
gootted cranstcion 0.5 (1 mp tanani ia ()))
respenees
sctal alartic.rerpate moothed tranissii) plasst_.reepoor
return Jiowas
caocthe trancition
del uquaLicnsqummss, 1): Depact the parsatters
stphs, beta, delta, +0, sigra.y, eignnLo, 9. Ts À, paraas
Zeppenehreddepeneet
sp.esp(D
# Bertic-plaftic trensiiion term with teperaaaar dependence sp.ncp(-delte ((a atrsinn
zipha
Field tore with lecx2 betu (1 pepldelle (fe wo s00l))
dependence
Smotth stop function for tranctcios. moooth_stop
tp sp-orpt-q
Ainisneor Aersiesin tare sare3 - 0.5 (1+ np.tanh(10 (a a0) moothh meothustep) Co <0032
Rerations
Serad
taras)
i Shocth aproeieatio for stress calcslation mmoohh, apprex
sp.eepp exp (smcott_atep 0..)))
Crebise termo find
stease tare: (turaz Serwi) (temns smcoth, (Sormt semes)) Gecoth, appees terad snooth, sppros) (ters) sermal)
Figure 22: Performance Trajectory LLM-SR (Mixtral) along the best-scoring equation program skeletons (before parameter optimization) over iterations on the Stress-Strain problem.
29

---

Published as a conference paper at ICLR 2025
def equaect (mprrrmm, dhguck Sh,
LLM-SR (Mixtral)
T.ept, pil.opt parame
Zergeratore factor Semperatura
ag.wxp(-Xt op.abatt T_opt) (CT T_opt) 2)
12
23)
pdf Jactor
Fll factor 1/ (1 Kh. (np abolpa p(_opt) F))
Inashetson factor
Inhibitioo_facter
wtppelpppppaheeeee pll_ope) wepp-np.pppeaeee
pl.opt)),
Taggeratore-depentert groath factor
h_factor T_opt)
grouch_factor factor up. Lanh fufactor "p.arciaa((_
fal_facoor tobibdtoon,fottor
Rinatte factor (Nonod)
factorr / (Ss
Grouth
grouth_rate
Kinetie_factor where(
kinetic_factor up opp (kiaetic_factee
kinetic. factor)
Coobine affectx ged final grouth rate en texperature foctor foctor T.opb) / 10,
frcttt up an.iniiam ap.abo(T
return rote
dof auuationpparama B, S, I, p(): Uspack the parancters
LLM-SR (GPT 3.5)
Na, alpba, Ka, Kp, Kd, gomna_t, gazmo_P, T_opt, pil_opt parana
Temperafure rodulalion factor tanp_ractor / (1 ep. appeeappeen_tt (T T_opt)))
all modulation factor pilfactor (1 ep.exp(gamm.p
plt_opt))
Poputation denasty effect population offect porppoeett <lphe)
enmenttraionn affect (Honad)
substrate, sabetrate_effecc effect = (8 Ka)
Tesperature dependency
tenperaturo_offec. ap.exp(-Kt ap.aqurroT T_apt))
plt dependency
pil_etfect mp. ap.exp(-Kp np opprropoppp pll_opt))
Racterial death death_rata
Conbena all offacta to get the growth nets
gate population
teaperabure_cfft pll_effect death_rate) tesp_factor factor
grosth_rate
(log(lo)(HH)
DSR
pl (b s.pll (log(h)
(-93257 x to 329004 10 for- 18890 10-10000-02000 50005 10 181ptt
* 642027 10 4.070 10 fast 678274 x 10 *asr 1.5007 10 48spat
78s 548619 "ay" 127722
Ppr
*881503 10 ABpas" 1.4007 x "Bpht s.mss x 20--23881111 x
Aphy 131312
1.768 10 "smpel 441799 x Agy 200799 x.10 Aspat" 2.5544 15phl
10-"s 439662
10-"rpt 420817 10-Prptt
=1.1249 x 10-2y- 122334 for" 4.8089 10-"pa"== -0- -22310)
xeep (tn p1t++
pHf
p+
sin 721 sintcos (in.0005))) 0.2))*
PySR
Figure 23: Final discovered equations from LLM-SR and other leading SR baseline methods (DSR, uDSR, PySR for E. coli bacterial growth rate problem.
30

---

Published as a conference paper at ICLR 2025
def equation (parane, 6, 17:
LLMSS (Mitral)
# lhpaak the parans Lers
alpha, bota, delta, e0, eigat.s, sigaa.0, 4, parane
factor
atrnin_scale
up. arR
Blastic plactic trassition term with temperature dependence
taral alpha a / (1 ng.exp(-dalta (fe ed atrain__cale)==2))
# Field stress term with bespersfare dependence aignly
up. sp.asp(-dalta
ob atroin_eall))=23)
Smooth step function for tranaition smoohh
up. exp (e x)))
# Hostinear hardening tens tarad 0.5
(10 (a e0) smooth. <ao0th_rtep)
gield term
tormo oigeo_0 terns)
i Smooth spprozinatios for stress caiculation smooth.
ap. 10 (ezcoth.stey 0.533)
Combise
fied strese
abreax tern (bare) term)) (saraS sppen (tarsd (tems)) (cmooth_appenn approx tarrd ceooth_spprax) (termz berxi))
return atrecc
daf equation (parane,, o, T): # Uspach the paraneterd
LLMSR (GPT 3.5)
22, 24, X6, 29, parsra
# Strass the clastic region alastic_etress
k2
# Stress in the plastic region alastic.atresa
2)
gStreza in the facture regson failure.atee k6 up. exppxx
(T
# Smooth transstson from clastic to plantsc regton tlattic.to.plattss
up aep(-8 (a k9)))
# Smooth tronsstsos froe plastic to failare region ap.axp(-k10
plattic,tt_ffilus
ka)))
# Stresses in atl regions
total_stzeas
olastic
elastic, toeeii_tt___oo plaasas))))
taattic.to_fallle
plastic_streza x clastic_to, paaiit__to_ppaaee plastic * (1 sosstic_to_fallrre
failure_strees plastic_to_feille
sin (exp e Txp (T²) +e +e- +e-r))
DSR
- 3.51ee + 7 3.812 0.92e - 0.30c 7 -409
+ 1.11 2017 087T log log (2²) 2.71
sin 1.136"
PySR
cos(e)
0.004
Figure 24: Final discovered equations from LLM-SR and leading SR baseline methods (DSR, uDSR, PySR) for Stress-Strain Strain problem.
31

---

Published as a conference paper at ICLR 2025
mnon
Time: (0,
Displacement
Velocity (V]
Figure 25: Qualitative evaluation of the performance of LLM-SR on Oscillation problem compared to uDSR
and PySR baselines, Plots show the target acceleration with respect time, displacement, and velocity.
(b)
Temperature: 300
Dinnal Pyen
Strss paPn)
5500
6.00
0211
0100
Figure 26: Qualitative evaluation of L.M--SR performance for Stress-Strain problem compared to uDSR and PySR haselines. Plots show the target stress with respect to strain and temperature.
32


Reviews:
[{"review_id": "jlNTNzasbs", "reviewer": "Reviewer_4Qbb", "summary": "This paper introduces LLM-SR, a new approach to discovering mathematical equations from data using LLMs' scientific knowledge and code generation. LLM-SR treats equations as program structures. LLM-SR achieves significant improvements over baseline symbolic regression methods.", "questions": "1. A more thorough analysis of the computational cost, especially comparing LLM-SR to traditional symbolic regression techniques, would help gauge the practical feasibility of this approach.\n\n2. Providing a real-world example of how LLM-SR assists in scientific problem solving could make the results more compelling.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1. Innovative approach: It is innovative to use LLMs for scientific equation discovery and symbolic regression. LLM-SR integrates LLMs' prior scientific knowledge and and code generation ability to achieve good performance.\n\n2. Robust evaluation: The authors have carefully designed benchmark problems that align well with the real-world applications.\n\n3. Comprehensive analysis: The qualitative analysis effectively evaluates the model's OOD generalizability  and the ablation study provides insight into the contribution of each of model's component.\n\n4. Well written: The paper is well written with clear data presentation.", "weaknesses": "The paper does not provide detailed insight into how LLM-SR handles edge cases, such as noisy data, or highly complex functions. This could be important for real-world applications where data may be imperfect or equations are inherently intricate.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mZrXklOXz9", "reviewer": "Reviewer_89Tg", "summary": "The paper introduces **LLM-SR**, a novel framework for discovering scientific equations using **Large Language Models (LLMs)**. The approach integrates LLMs' scientific knowledge and code generation capabilities with evolutionary search and optimization techniques to generate, evaluate, and refine mathematical equation hypotheses iteratively. The key steps involve generating equation skeletons, optimizing parameters using Python-based tools, and managing an experience buffer to iteratively improve the search process.\n\n**Main Contributions:**\n1. **Novel Framework**: LLM-SR combines LLMs' strengths with evolutionary algorithms to navigate complex equation discovery.\n2. **Integration of Scientific Priors**: The method leverages LLMs' embedded scientific knowledge for more efficient hypothesis generation.\n3. **Benchmark Creation**: The authors design new benchmark problems in physics, biology, and materials science to test the method's capabilities and avoid memorization of well-known equations.\n4. **Performance Advantage**: LLM-SR reportedly outperforms traditional symbolic regression methods, showing better accuracy and generalization, especially in out-of-domain tests.\n5. **Ablation Study**: Demonstrates the importance of various components, such as problem specification and iterative refinement, in the overall performance of LLM-SR.", "questions": "1. **Clarification on the Necessity of Complex Mechanisms**:\n   - *Question*: Could the authors elaborate on the rationale behind incorporating complex mechanisms such as the sampling strategy and island model? How do these elements contribute to the overall performance, especially compared to simpler alternatives?\n   - *Suggestion*: Including comparative results or ablation studies that illustrate the impact of these mechanisms on the method’s performance would be helpful to better understand their necessity and specific contributions.\n\n2. **Verification of LLM Prior Knowledge Utilization**:\n   - *Question*: How do the authors verify that the LLM’s embedded scientific prior knowledge is effectively utilized during the equation generation process? Are there examples or analyses demonstrating cases where this prior knowledge played a significant role?\n   - *Suggestion*: We suggest adding controlled experiments that isolate and assess the contribution of prior knowledge. Comparing models with and without domain-specific pretraining might further strengthen this aspect.\n\n3. **Evaluation on Real-World Datasets with Noise and Outliers**:\n   - *Question*: Do the authors plan to evaluate LLM-SR on more challenging, real-world datasets that contain noise, outliers, or measurement errors? How do they foresee the method performing under such conditions?\n   - *Suggestion*: Expanding the experiments to include data with inherent noise and outliers, or discussing the method’s anticipated robustness in these cases, could help demonstrate its applicability to complex, real-world scenarios.\n\n4. **Impact of LLM Programming Performance on Results**:\n   - *Question*: Have the authors considered examining the extent to which the LLM’s programming capabilities, rather than scientific reasoning, might influence the results? How do they ensure that the outcomes reflect genuine equation discovery rather than sophisticated code generation?\n   - *Suggestion*: Providing an analysis or qualitative examples that differentiate between the model’s programming skills and its scientific reasoning in generating equations would be valuable to clarify the contributions and limitations of the approach.\n\n5. **Interpretability of Results**:\n   - *Question*: The authors state that the generated equations are more interpretable and scientifically meaningful. Could they provide examples or a detailed comparison showing how LLM-SR’s output stands out in terms of interpretability relative to baseline methods?\n   - *Suggestion*: Presenting case studies or examples where the equations discovered by LLM-SR provide clearer or more scientifically insightful interpretations would reinforce this claim.\n\n6. **Broader Comparison with Hybrid Methods**:\n   - *Question*: How does LLM-SR compare with recent hybrid approaches that incorporate domain knowledge or leverage deep learning for symbolic regression? Are there particular strengths or limitations of LLM-SR compared to these methods?\n   - *Suggestion*: Including comparisons with other state-of-the-art hybrid methods or discussing potential improvements and challenges in adapting LLM-SR could provide more context on its positioning and competitive advantages.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. **Innovative Application of LLMs for Scientific Equation Discovery**:\n   - The authors introduce an innovative approach by leveraging the scientific knowledge and code generation capabilities of large language models (LLMs) for equation discovery. Representing equations as code and integrating data-driven feedback and optimization provides a novel perspective for scientific modeling.\n\n2. **Experience Buffer and Iterative Optimization**:\n   - The use of an experience buffer to maintain high-quality equation examples and facilitate iterative optimization is a unique feature. This mechanism helps the model learn and recall effective generation patterns, improving the quality and efficiency of equation exploration.\n\n3. **Comprehensive Multi-Domain Benchmarking**:\n   - The paper demonstrates the method’s applicability across multiple scientific domains, including physics, biology, and materials science. The newly designed benchmarks, which aim to prevent memorization of known equations, reflect the authors' commitment to scientific rigor.\n\n4. **Promising Experimental Results**:\n   - Despite certain challenges and areas for improvement, LLM-SR shows promising performance compared to existing symbolic regression methods, particularly in out-of-domain tests. This indicates potential for tackling complex scientific problems and suggests a valuable direction for future research.", "weaknesses": "**Complexity and Necessity of the Method**:\nThe authors introduce complex mechanisms such as sampling strategies and an island model, showcasing innovative efforts. However, clearer justification is needed to explain the necessity of these elements and their specific contributions to enhancing the method’s performance. Simplifying and elucidating these mechanisms' actual benefits would improve the clarity and persuasiveness of the method.\n\n**Utilization and Verification of Prior Knowledge**:\nAlthough the authors assume that the LLM leverages its scientific prior knowledge to generate reasonable equations, the paper lacks detailed analysis and experimental validation of how this mechanism functions and its tangible effects. Further exploration and evidence would more convincingly support this claim and make the method’s core advantage clearer to readers.\n\n**Breadth and Realism of Experimental Design**:\nThe current experiments are based on data generated from known equations, providing a basic evaluation. However, data in real-world scientific research often contains noise and outliers, and the method’s performance under such conditions is a critical indicator of its practical applicability. Including datasets with measurement errors and outliers would help demonstrate the method’s robustness and generalizability in complex, realistic scenarios.\n\n**Impact of Programming Capabilities on Experimental Results**:\nThe paper emphasizes the LLM's ability to generate programming code but lacks analysis of whether these programming capabilities interfere with the experimental results. The boundary between the LLM’s programming skills and scientific reasoning remains unclear, which might mean that experimental outcomes reflect coding proficiency rather than pure scientific inference. Clarifying the role of these aspects in equation generation would aid in better understanding the method’s true strengths and limitations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "RlYB0Qp5Lj", "reviewer": "Reviewer_BtD7", "summary": "This paper proposes a method for symbolic regression (SR) that combines LLM-based evolutionary search for program synthesis with numerical optimizers (such as BFGS and Adam) to tune the parameters of the generated programs. The authors, additionally, present four new SR problems that were created with synthetic modifications to address concerns regarding LLM memorization. Finally, ablations are provided to analyze the benefits of different design choices.", "questions": "1. L243 briefly described how population initialization is performed. Is the linear equation skeleton manually written? How many skeletons are provided at initialization?  \n2. In L149, should the expectation be over the data points in $\\mathcal{D}$?  \n3. In Appendix E, could you clarify what LLM-SR (w/o Prior) refers to? From the text described in the main paper, I could not find mention of the described \"problem-specific prior knowledge\" (L1220).", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. To the best of my knowledge, framing SR as equation _template_ discovery seems novel and a useful contribution.\n2. The analyses presented to demonstrate LLM memorization of the Feynman dataset are quite convincing.\n3. The new benchmark problems are well-motivated and look like sound contributions.\n4. The paper is, overall, clearly written and easy to follow, and visualizations are used well in various places.", "weaknesses": "My main concerns with the work are twofold:\n1. The paper reads in a manner that moderately overclaims novelty. The main idea appears to be the combination of numerical optimizers with LLM-based optimization over templates, besides the contribution of new SR problems. However, the motivation, currently, is primarily devoted to comparisons with non-LLM methods, which seems unfair. Since the framing of SR problems within the context of LLM-based optimization is not new (Table 1 from [1] is a useful reference), it would be more appropriate to motivate the work by describing existing LLM-based SR methods, clearly stating their problems, and then showing how the new method addresses them.  \n2. My second concern is regarding the choice of baselines. While I appreciate the inclusion of several standard SR methods, I think the comparisons do not precisely show the benefits of the paper's contributions.  \n  (a) Given prior works such as [1, 2, 3] in LLM-based SR/optimization, it would have provided utility to see the performance of the proposed method in comparison to its closest counterparts. Alternatively, if the argument is that some of the ablations subsume these methods, it would then be useful to see those ablations presented as the main baselines on all the datasets instead of only on a partial set.  \n  (b) The ablation results from Figure 6 show that a substantial drop in performance occurs when coefficient optimization is removed. Based on this, it would have been useful to see whether simply equipping the non-LLM baselines with a similar coefficient optimizer would close the gaps that we currently see in Table 1. I am currently not convinced that the LLM is needed to facilitate the evolutionary search of equation templates.  \n\nMinor:\n1. The choice of describing the contents in Section 2.4 as \"Experience Management\" is mildly concerning. Given a large body of prior work in evolutionary computation, it would be appropriate to avoid fragmenting the literature for the benefit of differentiation and instead re-use existing terminology. If there is indeed an aspect that warrants a new term, please do correct me.  \n2. It was mildly odd to read that mathematical equations have been \"unreasonably\" effective in describing complex phenomena. Unclear why we think it is unreasonable.  \n\n[1] Meyerson et al., 2023. Language model crossover: Variation through few-shot prompting.  \n[2] Romera-Paredes et al., 2023. Mathematical discoveries from program search with large language models.  \n[3] Yang et al., 2024. Large Language Models as Optimizers.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "A7BHTPlNGP", "reviewer": "Reviewer_AFGC", "summary": "The paper presents LLM-SR (Large Language Model Symbolic Regression), a framework that employs LLMs to discover scientific equations by treating equations as programs. LLM-SR combines LLMs with evolutionary search techniques to generate, evaluate, and refine equations that describe complex scientific phenomena. The process includes hypothesis generation by the LLM, parameter optimization of proposed equations, and a dynamic experience management system to refine hypotheses iteratively. The authors demonstrate LLM-SR’s superior ability to discover accurate and generalizable equations compared to traditional symbolic regression methods, especially in out-of-domain settings.", "questions": "1. How would LLM-SR perform with even larger LLMs or specialized scientific models?\n2. How does LLM-SR perform with complex systems involving interactions across multiple equations?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The approach of using LLMs to generate symbolic equation structures combines programming and scientific prior knowledge, setting it apart from traditional symbolic regression methods.\n2. LLM-SR shows strong performance not only in-domain but also in out-of-domain scenarios, which is often challenging for symbolic regression techniques.\n3. By leveraging scientific prior knowledge and a structured experience buffer, LLM-SR reduces the search space, achieving good results with fewer iterations compared to traditional methods.", "weaknesses": "1. The framework relies heavily on the embedded scientific knowledge of LLMs, which may be biased, incomplete, or inaccurate for certain scientific domains, potentially limiting its robustness.\n2. Although the paper addresses different domains, it would benefit from more diverse and complex real-world datasets to further validate the method’s general applicability.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)