LLM-SR: Scientific Equation Discovery via Programming with Large Language Models

## Abstract

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 high-dimensional combinatorial and nonlinear hypothesis spaces. Traditional methods of equation discovery, commonly known as symbolic regression, largely focus on extracting equations from data alone, often neglecting the rich domain-specific prior knowledge that scientists typically depend on. 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 in an efficient manner. 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 physical understanding, which are then optimized against data to estimate skeleton parameters.
We demonstrate LLM-SR ’s effectiveness across three diverse scientific domains, where it
discovers physically accurate equations that provide significantly better fits to in-domain and out-of-domain data compared to the well-established symbolic regression baselines. Incorporating scientific prior knowledge also enables LLM-SR to search the equation space more efficiently than baselines 1 1 1 Code is available at: https://github.com/deep-symbolic-mathematics/LLM-SR .

## 1 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., ). As LLMs continue to evolve, researchers are exploring innovative ways to harness their potential for solving complex problems such as scientific discovery (Wang et al., ; AI4Science & Quantum, ).
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., ).
Additionally, by leveraging their ability to understand and reason with the help of programming and execution, LLMs have shown the potential to enhance automatic reasoning and problem-solving capabilities (Romera-Paredes et al., ; Madaan et al., ).
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 observations such that these equations are predictive of the underlying data, are interpretable, and generalize to unseen data.
Finding such equations offers several advantages over simply estimating a predictive model, as the resulting mathematical functions provide insights into the underlying physical processes, enable extrapolation beyond the observed data, and facilitate knowledge transfer across related problems (Schmidt & Lipson, ).
However, while evaluating the fit of a proposed equation is relatively straightforward, the inverse process of obtaining these mathematical equations from data is a challenging problem, known to be NP-hard (Virgolin & Pissis, ).
Current equation discovery methods encompass a wide variety of approaches from evolutionary search algorithms (Cranmer, ; Mundhenk et al., ; La Cava et al., ) to advanced deep learning methods using Transformers
(Biggio et al., ; Kamienny et al., ).
Most of the traditional symbolic regression techniques are built on top of Genetic Programming (GP) (Koza, ) evolutionary methods, representing mathematical equations as expression trees and searching the combinatorial space of possible equations through iterative mutation and recombination.
However, current 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. Thus, there is a need for equation discovery methods that effectively integrate prior scientific knowledge into the navigation of equation search space, a strategy akin to a scientist’s reliance on foundational knowledge when formulating hypotheses for scientific discovery.

Figure: Figure 1: The LLM-SR framework, consisting of three main steps: (a) Hypothesis Generation, where LLM generates equation program skeletons based on a structured prompt; (b) Data-driven Evaluation, which optimizes the parameters of each equation skeleton hypothesis and assesses its fit to the data; and (c) Experience Management, which maintains a diverse buffer of high-scoring hypotheses to provide informative in-context examples into LLM’s prompt for iterative refinement.
Refer to caption: extracted/5638465/fig_table_file/LLMSR-v9.jpg

To address these limitations, we introduce LLM-SR (Fig. 1), a novel framework that combines the strengths of LLMs, reliable optimizers, and evolutionary search for data-driven equation discovery.
At its core, LLM-SR is an iterative hypotheses refinement method that iteratively generates hypotheses, evaluates these hypotheses, and uses the evaluation signal to refine and search for better hypotheses. Specifically, given an initial equation hypothesis for the underlying data, LLM-SR utilizes LLMs to propose new equation hypotheses ( 1(a)), evaluates their fit on the data using off-the-shelf optimizers ( 1(b)), and uses this data-driven feedback and a carefully maintained dynamic memory of previous equations ( 1(c)) to 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 a given system to generate educated hypotheses for equation skeletons (i.e., equation forms with placeholder parameters for numeric coefficients and constants). The LLM’s adaptive and informed mutation and crossover capabilities (Meyerson et al., )
are then employed to refine the suggested equations in an iterative process. LLM-SR incorporates data-driven feedback into the search process by evaluating the fitness of the generated equations against the observed data, guiding the search towards more accurate equation skeletons.
By representing equations as programs, we take advantage of LLM’s ability to generate structured and executable code (Li et al., ; Shojaee et al., ) while providing a flexible and effective way to represent general mathematical relations. This representation also facilitates direct and differentiable parameter optimization to better optimize the coefficients or constants in the generated equations.
We evaluated LLM-SR using GPT-3.5-turbo and Mixtral-8x7B-Instruct (Jiang et al., ) as backbone LLMs
across four equation discovery problems spanning three scientific domains.
LLM-SR consistently outperforms state-of-the-art symbolic regression methods, discovering physically accurate equations with better fit and generalization in both in-domain (ID) and out-of-domain (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 analysis also highlights the crucial role of iterative refinement and evolutionary search in LLM-SR’s superior performance.
The major contributions of this work can be summarized as follows:

- •
We introduce LLM-SR, a novel method that combines the strengths of LLMs, off-the-shelf optimizers, and evolutionary search for data-driven equation discovery.
- •
We show that LLM-SR outperforms state-of-the-art SR methods by navigating the equation search space more efficiently and discovering more accurate equations with better out-of-domain generalization across 4 problems from 3 scientific domains.
- •
We demonstrate through a comprehensive ablation study that providing natural language descriptions of the problem and variables, data-driven feedback, and skeleton parameter optimization signals are vital 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 a concise and interpretable symbolic mathematical expression $\tilde{f}$ that closely approximates an unknown underlying equation $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$, which maps a $d$-dimensional input vector $\mathbf{x}$ to output $y$. Given a dataset $\mathcal{D}=\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{n}$ consisting of input-output pairs, SR methods seek to uncover the hidden mathematical relationship, such that $\tilde{f}(\mathbf{x}_{i})\approx y_{i}$ for all $i$. 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, ), prefix sequences (Petersen et al., ; Biggio et al., ), or context-free grammars (Brence et al., ), constructing a search space $\mathcal{F}$. These representations provide limited and structured search spaces, enabling evolutionary search approaches like genetic programming (GP) to explore and find candidate expressions. In contrast, our work employs program functions to directly map inputs $\mathbf{x}$ to output targets $y$:
def f(x): ... return y.
This offers a flexible and expressive way to represent equations, allowing for a diverse set of mathematical operations and step-by-step instructions to capture complex mathematical relations. However, the space of possible programs can be vast, with extremely sparse solutions, which can hinder the effectiveness of traditional search methods like GP in finding the desired equation programs (Devlin et al., ; Parisotto et al., ).
To tackle this challenge, we leverage the scientific knowledge and coding capabilities of LLMs to effectively navigate the program space for equation discovery.
LLMs have the ability to generate syntactically correct and knowledge-guided code snippets, guiding the efficient exploration of the large equation program optimization landscape.

Let $\pi_{\theta}$ denote a pre-trained LLM with parameters $\theta$. In this task, we iteratively sample a set of equation program skeletons from LLM, $\mathcal{F}=\{f:f\sim\pi_{\theta}\}$. The goal is to find the skeleton that maximizes the reward (evaluation score) $\text{Score}_{\mathcal{T}}(f,\mathcal{D})$ for a given scientific problem $\mathcal{T}$ and data observations $\mathcal{D}$. The optimization problem can be expressed as: $f^{*}=\arg\max_{f}\mathbb{E}_{f\in\mathcal{F}}\left[\text{Score}_{\mathcal{T}}
(f,\mathcal{D})\right]$,
where $f^{*}$ denotes the optimal equation program discovered at the end of the LLM-guided search.
LLM-SR prompts LLM to propose hypotheses given the problem specification,
and experience demonstrations containing previously discovered promising equations. In LLM-SR, LLM only generates hypotheses as equation program skeletons where the coefficients or constants of each equation are considered as a placeholder parameter vector that can be optimized. The parameters are optimized using robust off-the-shelf optimizers and evaluated for their fit. Promising hypotheses are then added to a dynamic experience buffer, which guides the equation refinement process, while non-promising hypotheses are discarded. Below we explain the key components of this framework, shown in Fig. 1, in more detail.

### 2.2 Hypothesis Generation

As shown in Fig. 1 (a), the hypothesis generation step relies on a pre-trained LLM to propose diverse and promising equation program skeletons. By leveraging the LLM’s extensive knowledge and generative capabilities, we aim to efficiently explore the vast space of potential equations while maintaining syntactic correctness and scientific plausibility.

Figure: Figure 2: Example of an input prompt for E. Coli bacterial growth with problem specifications, evaluation and optimization function, and the initial input equation example.
Refer to caption: extracted/5638465/fig_table_file/spec_bact7.jpg

#### 2.2.1 Prompt

To guide the LLM in generating relevant and well-structured equation program skeletons, we design prompts that provide essential context and specify the desired output format. Our prompt structure, shown in Fig. 2, consist of the following components:

Instruction.
We begin the prompt with a clear instruction to the LLM, requesting the completion of the last function’s body in the provided code. The instruction emphasizes the importance of considering the physical meaning and relationships of the input variables while completing the function. This guidance helps steer the LLM’s output towards scientifically meaningful and coherent suggestions.

Problem Specification.
Following the instruction, we include a concise description of the scientific problem at hand. This specification outlines the key variables, constraints, and objectives of the problem, providing the necessary context for the LLM to generate relevant equation program skeletons. By clearly defining the problem space, we enable the LLM to focus its efforts on the most pertinent aspects of the task.

Evaluation Function.
To ensure that the generated equation program skeletons can be effectively assessed, we include the evaluation function in the prompt. This function includes the criteria used to measure the quality and fitness of the proposed equations. By providing the evaluation function upfront, we enable the LLM to generate skeletons that are aligned with the desired performance objectives.

Experience Demonstration.
To further guide the LLM and provide concrete examples of equation program skeletons as well as their improvement trajectory, we incorporate experience demonstrations in the form of in-context examples (Fig. 2 only shows the initial example). These examples serve as a reference for the LLM to build upon and refine in its generated outputs.

#### 2.2.2 Hypothesis Sampling

At each iteration $t$, we sample a batch of $b$ equation skeletons $\mathcal{F}_{t}=\{f_{i}\}_{i=1}^{b}$ from the LLM $\pi_{\theta}$, where each equation skeleton $f_{i}$ is generated based on the constructed prompt $\mathbf{p}_{t}$: $f_{i}\sim\pi_{\theta}(\cdot|\mathbf{p}_{t})$. To encourage diversity in the generated equation program skeletons, we utilize stochastic temperature-based sampling, which allows for a balance between exploration (creativity) and exploitation (prior knowledge) in the hypothesis space.

After obtaining the sampled equation programs and before processing to the data-driven evaluation step, we execute each sampled equation program $f_{i}$ and discard those that fail to execute successfully.
By discarding such equation programs, we ensure that only valid and executable hypotheses are considered for further evaluation. Also, to maintain computational efficiency and prevent excessively time-consuming evaluations, we impose a maximum execution time threshold for each equation program. If the execution time of an equation program exceeds this threshold, we discard the corresponding hypothesis.

### 2.3 Data-driven Evaluation

Following the hypothesis generation step, we evaluate and score the generated equation skeleton hypotheses using observed data. In this step, demonstrated in Fig. 1 (b), we first optimize the parameters of each hypothesis with regard to the dataset, and then use the optimized parameters to evaluate and score each hypothesis.

#### 2.3.1 Hypothesis Optimization

Inspired by (Biggio et al., ), in LLM-SR, we decouple the equation discovery process into two steps: ($i$) discovering the equation program structures (skeletons) using the LLM, and ($ii$) optimizing the skeleton parameters/coefficients based on data. The LLM is responsible for generating equation skeletons and the core logic of the program, while the numeric values of the parameters are represented as placeholders in the form of a parameter vector params (as shown in Fig. 2). These placeholders are subsequently optimized to fit the observed data.
Let $\mathcal{F}_{t}$ denote the set of equation program skeletons generated by the LLM at iteration $t$.
Each equation program skeleton $f\in\mathcal{F}_{t}$ is a function that takes both input features $\mathbf{x}$ and a vector of learnable parameters params as input arguments and returns the predicted target values:
“def f(x, params): ... return y”.

To optimize the parameters for each skeleton hypothesis, we use two approaches: (1) numpy+BFGS and (2) pytorch+Adam. The BFGS algorithm (Fletcher, ), available in the scipy library (as represented in Fig. 2), is a gradient-free nonlinear optimization method. In contrast, the Adam optimizer (Kingma & Ba, ), available in pytorch, is a direct gradient-based optimization algorithm.
The program notion of equation skeletons employed in LLM-SR framework enables the use of differentiable programming techniques, which is not possible with other representations common in equation discovery such as expression trees or prefix notations.
The choice of the optimization approach also depends on the specific characteristics of the problem and the equation skeleton.
The numpy+BFGS method is often preferred when the number of parameters is relatively small, while pytorch+Adam is more suitable for larger-scale problems and complex equation structures that benefit from efficient gradient computation through differentiable programming.

#### 2.3.2 Fitness Assessment

After finding the optimized values of the skeleton parameters, we assess the fitness of each equation program hypothesis by measuring its ability to capture the underlying patterns in the data. To obtain the predicted target values $\hat{\mathbf{y}}$, we replace the placeholder parameter vector params in the equation program skeleton $f$ with the optimized parameter values $\texttt{params}^{*}$ and evaluate the equation skeleton function over the input data: $\hat{\mathbf{y}}=f(\mathbf{x},\texttt{params}^{*})$.
We then compute the negative Mean Squared Error (MSE) between the predicted and true target values as the fitness evaluation score: $s=\text{Score}_{\mathcal{T}}(f,\mathcal{D})=-\text{MSE}\left(\hat{\mathbf{y}},
\mathbf{y}\right)$, where $s$ represents the evaluation score or reward corresponding to the equation skeleton $f$.

### 2.4 Experience Management

To better navigate the search landscape and avoid local minima, LLM-SR contains an experience management step, shown in Fig. 1 (c), which maintains a diverse population of high-quality equation programs in the experience buffer and effectively samples from this population to construct informative prompts for the LLM in next iteration.
This step consists of two key components: a) experience manager for maintaining the experience buffer, and b) sampling from this experience buffer.

#### 2.4.1 Buffer Maintenance

Following the data-driven evaluation, the pairs of equation skeleton hypotheses and their corresponding scores, denoted as $(f,s)$, are stored in an experience buffer $\mathcal{P}_{t}$ for further refinement in subsequent iterations of the search process.
To avoid local optima, we follow (Cranmer, ; Romera-Paredes et al., ) and adopt an islands model, also known as a multi-population model for managing the experience buffer.
We initialize $m$ islands with a copy of the equation program example provided in the initial prompt (equation$_$v0 in Fig. 2).
These islands evolve independently, allowing for parallel exploration of different regions in the equation program space.
At each iteration $t$, the new equation program hypotheses $\mathcal{F}_{t}$, along with their scores, $\{(f,s):f\in\mathcal{F}_{t},s=-\text{Score}_{\mathcal{T}}(f,\mathcal{D})\}$, are added to the same island from which the in-context examples of prompts were sampled, if the new score $s$ is better than the current best.
Within each island, we further cluster the equation programs based on their signature, which is defined as their score. Equation programs with identical signatures are grouped together.
This clustering approach helps preserve diversity by ensuring that equation programs with different performance characteristics are maintained in the population.

#### 2.4.2 Experience Sampling

To construct informative prompts for the LLM, we sample equation programs from the experience buffer, updating the prompt with new experience demonstration in-context examples.
As per (Romera-Paredes et al., ), this process follows a two-stage sampling method: ($i$) selecting a random island from the $m$ available ones, and ($ii$) sampling $k$ equation programs from the selected island for inclusion as in-context examples in the prompt.
Cluster selection follows a Boltzmann sampling (Maza & Tidor, ), with a score-based probability of choosing cluster $i$:
$P_{i}=\frac{\exp\left(s_{i}/\tau_{c}\right)}{\sum_{i^{\prime}}\exp\left(s_{i}^
{\prime}/\tau_{c}\right)}$, where $s_{i}$ denote mean score of the $i$-th cluster, and $\tau_{c}$ is a temperature parameter.
After cluster selection, individual equation programs are sampled, favoring shorter programs, using the probability proportional to $\exp\left(-\tilde{l}_{i}/\tau_{p}\right)$, where $\tilde{l}_{i}$ is the normalized program length, and $\tau_{p}$ is a temperature parameter.
More details of these sampling steps are provided in App. B.
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 program hypotheses.

Figure: Algorithm 1 LLM-SR

Algorithm 1 presents a summarized pseudo-code of the LLM-SR framework. We first initialize the experience buffer population $\mathcal{P}_{0}$ using InitPop with simple linear initial equation skeleton example (as shown in Fig. 2 for E. coli growth problem). This equation skeleton serves as a starting point template for generating candidate equations, avoiding the need for manual prior knowledge incorporation. At each iteration $t$, SampleLLM (detailed in Algorithm 2), which represents Hypothesis Generation, samples $b$ equation program skeletons from the LLM using an updated prompt with $k$ in-context examples from $\mathcal{P}_{t-1}$. The sampled equations are evaluated using $\text{Score}_{\mathcal{T}}(f,\mathcal{D})$, and $\mathcal{P}_{t}$ is updated with the equation-score pairs $(f,s)$. In Algorithm 2, SampleDB function which represents Experience Sampling is first called to sample in-context examples from the experience buffer and update prompt accordingly. Then, LLM samples are collected with the updated prompt. After $T$ iteration, the best-scoring program $f^{*}$ from $\mathcal{P}_{t}$ and its corresponding score $s^{*}$ are returned as the optimal solution found for the problem.
LLM-SR iteratively refines equation skeletons using the LLM’s generative capabilities, guiding the search towards promising structures based on the evolving experience buffer population.

## 3 Experimental Setup

### 3.1 Benchmarks and Datasets

Current symbolic regression benchmarks (La Cava et al., ; Strogatz, ; Udrescu & Tegmark, ), mostly based on Feynman physics equations (Udrescu et al., ), may not effectively capture the true discovery process and reasoning capabilities of LLMs. Our findings show that LLM-SR rapidly achieves low Normalized MSE scores within very few iterations on Feynman problems, suggesting that LLMs have likely memorized these commonly known equations due to their prevalence in training data (check App. C for full results). Feynman equations are often commonly known, and easily memorized by LLMs, making them unsuitable evaluation benchmarks for this study.
To overcome this limitation, we introduce novel benchmark problems across three diverse scientific domains that are designed to challenge the model’s ability to uncover complex mathematical relations while leveraging its scientific prior knowledge.
These new benchmark problems are designed to simulate the conditions for scientific discovery,
ensuring that the equations cannot be trivially obtained through memorization of the language model. By focusing on diverse scientific domains and introducing custom modifications to known equations, we encourage LLM-SR to creatively explore the equation search space and identify mathematical relations. We next discuss these new benchmark problems.

Nonlinear Oscillators.
Nonlinear damped oscillators are ubiquitous in physics and engineering. The dynamics of oscillators are governed by differential equations that describe the relationship between the oscillator’s position, velocity, and the forces acting upon it.
Specifically, the oscillator’s motion is given by a second-order differential equation: $\ddot{x}+f(t,x,\dot{x})=0$, involving time ($t$), position ($x$), and nonlinear function $f(t,x,\dot{x})$ accounts for forces.
In this study, we explore two custom nonlinear forms: Oscillation 1: $\dot{v}=F\sin(\omega x)-\alpha v^{3}-\beta x^{3}-\gamma x\cdot v-x\cdot\cos(x)$; and Oscillation 2: $\dot{v}=F\sin(\omega t)-\alpha v^{3}-\beta x\cdot v-\delta x\cdot\exp(\gamma x)$, with $v=\dot{x}$ as velocity and ${\omega,\alpha,\beta,\delta,\gamma}$ as constants. These forms are designed to challenge the model’s ability to uncover complex equations, avoiding simple memorization of commonly known oscillators.
More details on the design and data generation of these oscillators are provided in App. D.1.

Bacterial Growth.
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 (Monod, ; Rosso et al., ). Discovering equations governing E. coli growth rate under different conditions is crucial for predicting and optimizing bacterial growth (Tuttle et al., ).
The bacterial population growth rate can be modeled using a differential equation with the effects of population density ($B$), substrate concentration ($S$), temperature ($T$), and $pH$ level, which is commonly formulated as $\frac{dB}{dt}=f(B,S,T,pH)=f_{B}(B)\cdot f_{S}(S)\cdot f_{T}(T)\cdot f_{pH}(pH)$.
To invoke the LLM’s prior knowledge about this problem, yet avoid memorization, we consider nonlinear arbitrary designs for $f_{T}(T)$ and $f_{pH}(pH)$ that share similar characteristics with the common models but are not trivial to recover from memory. More details on the specific differential equation employed and data for this problem are provided in App. D.2.

Material Stress Behavior.
Understanding the stress-strain behavior of materials under various conditions, particularly the relationship between stress, strain, and temperature, is crucial for designing and analyzing structures in various engineering fields. In this problem, we use a real-world dataset from (Aakash et al., ), containing tensile tests on aluminum
at six temperatures ranging from 20°C to 300°C. The data covers the elastic, plastic, and failure regions of the stress-strain curves, providing valuable insights into the material’s behavior. More details on the data and visualization of these stress-strain curves can be found in App. D.3.

### 3.2 Baselines

We compare LLM-SR against several state-of-the-art symbolic regression (SR) baselines, including GPlearn(^1^11[https://gplearn.readthedocs.io/en/stable/](https://gplearn.readthedocs.io/en/stable/)), Deep Symbolic Regression (DSR) (Petersen et al., ), Unified DSR (uDSR) (Landajuela et al., ), and PySR (Cranmer, ). GPlearn is a pioneering Genetic Programming (GP) based SR approach with the open-source gplearn package. DSR is a deep learning-based SR approach that employs reinforcement learning over symbolic expression sequence generations. uDSR extends DSR by incorporating large-scale pretraining with Transformers and neural-guided GP search at the decoding stage.
PySR(^2^22[https://github.com/MilesCranmer/PySR](https://github.com/MilesCranmer/PySR)) is a recent SR method that uses multi-island asynchronous GP-based evolution with the aim of scientific equation discovery.
For a fair comparison with LLM-SR, we allow all baselines to run for over $2$M iterations until convergence to their best performance.
More details on the implementation and parameter settings of each baseline are provided in App. A.

### 3.3 LLM-SR Configuration

We experiment with GPT-3.5-turbo and Mixtral-8x7B-Instruct as 2 different LLM backbones for LLM-SR. At each iteration, we sample $b=4$ equation skeletons per prompt from the LLM using a temperature $\tau=0.8$, optimize equation skeleton parameters via BFGS ($30$s timeout per program), and sample $k=2$ in-context examples from the experience buffer for the refinement.
We run LLM-SR variants for $2.5$K iterations in all experiments.
For more implementation details (incl. the experience buffer structure, prompt refinement strategy, and parallel evaluation), please refer to App. A.

## 4 Findings

### 4.1 LLM-SR Discovers more Accurate Equations

Table 1 shows the performance comparison
of LLM-SR (using GPT-3.5 and Mixtral backbones) against state-of-the-art symbolic regression methods on various scientific benchmark problems.
In these experiments, symbolic regression methods were allowed to run for considerably longer iterations (over $2$M) while LLM-SR variants ran for around $2$K iterations.
The table reports the Normalized Mean Squared Error (NMSE), where lower values indicate better performance.
Under the in-domain (ID) test setting, where the test set is sampled from the same domain as the data used for equation discovery, LLM-SR variants significantly outperform leading symbolic regression baselines across all benchmarks, despite running for far fewer iterations. Among the symbolic regression baselines, PySR and uDSR exhibit superior ID performance compared to DSR and GPlearn.

**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 Mean Squared Error. ’All’ column indicates average scores over all datasets. SR baselines ran for over $2$M iterations, while LLM-SR variants ran for around $2$K iterations.**
| Model | Oscillation 1 | Oscillation 2 | E. coli growth | Stress-Strain | All |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| ID$\downarrow$ | OOD$\downarrow$ | ID$\downarrow$ | OOD$\downarrow$ | ID$\downarrow$ | OOD$\downarrow$ | ID$\downarrow$ | OOD$\downarrow$ | ID$\downarrow$ | OOD$\downarrow$ |  |
| GPlearn | 0.0155 | 0.5567 | 0.7551 | 3.188 | 1.081 | 1.039 | 0.1063 | 0.4091 | 0.4894 | 1.298 |
| DSR (Petersen et al., ) | 0.0087 | 0.2454 | 0.0580 | 0.1945 | 0.9451 | 2.4291 | 0.3326 | 1.108 | 0.3361 | 0.9942 |
| uDSR (Landajuela et al., ) | 0.0003 | 0.0007 | 0.0032 | 0.0015 | 0.3322 | 5.4584 | 0.0502 | 0.1761 | 0.0964 | 1.409 |
| PySR (Cranmer, ) | 0.0009 | 0.3106 | 0.0002 | 0.0098 | 0.0376 | 1.0141 | 0.0331 | 0.1304 | 0.0179 | 0.3662 |
| LLM-SR (Mixtral) | 7.89e-8 | 0.0002 | 0.0030 | 0.0291 | 0.0026 | 0.0037 | 0.0162 | 0.0946 | 0.0054 | 0.0319 |
| LLM-SR (GPT-3.5) | 4.65e-7 | 0.0005 | 2.12e-7 | 3.81e-5 | 0.0214 | 0.0264 | 0.0210 | 0.0516 | 0.0106 | 0.0196 |

### 4.2 LLM-SR has better OOD generalization

Figure: Figure 3: Comparison of E. coli growth rate distributions from LLM-SR, PySR, and uDSR. LLM-SR aligns well with the ground truth, even for OOD data (unshaded regions), demonstrating better generalization than PySR and uDSR, which overfit the ID data (shaded regions and black points).
Refer to caption: extracted/5638465/fig_table_file/qual_bacterial2.jpg

The ability to generate equations that generalize well to domains beyond the training data, i.e., out-of-domain (OOD) data, is a crucial characteristic of effective equation discovery methods. To assess the generalization capability of the predicted equations, we evaluate LLM-SR and the symbolic regression baselines in the OOD test setting, focusing on their ability to predict systems outside the region encountered during training. Table 1 reveals that the performance gap between LLM-SR and the baselines is more pronounced in the OOD test setting. This observation suggests that equations discovered by LLM-SR exhibit superior generalization capabilities, likely attributable to the scientific prior knowledge embedded by LLMs in the equation discovery process.
For instance, on the E. coli growth problem, LLM-SR achieves an OOD NMSE of around 0.0037, considerably better than symbolic regression methods which all obtain OOD NMSE $>1$. Fig. 3 also compares the ground truth distribution of E. coli growth problem with the predicted distributions obtained from LLM-SR, PySR, and uDSR.
The shaded region and black points indicate in-domain (ID) data, while the rest represent out-of-domain (OOD).
Results show that the 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, leading symbolic regression methods like PySR and uDSR tend to overfit the observed data, with their predicted distributions deviating significantly from the true distributions in OOD regions.
This overfitting behavior highlights their limited ability to generalize beyond the training data and capture the true physical underlying patterns.
For more detailed results and analyses for other benchmark problems, check App. E.

### 4.3 LLM-SR Discovers Equations More Efficiently

Figure: Figure 4: Best score trajectories of LLM-SR with GPT-3.5 and Mixtral against SR baselines across different benchmark problems. LLM-SR discovers accurate equations more efficiently, requiring fewer iterations to outperform baselines, which fail to match LLM-SR even after $2$M iterations.
Refer to caption: extracted/5638465/fig_table_file/compare_curves-v7.png

Fig. 4 shows the performance trajectories of LLM-SR variants and symbolic regression baselines across different scientific benchmark problems, depicting the maximum fitting scores achieved over search iterations.
By leveraging scientific prior knowledge, LLM-SR needs to explore a considerably lower number of equation candidates in the vast optimization space compared to symbolic regression baselines that lack this knowledge. This is evident from the sharp drops in the error curves for LLM-SR variants, indicating they efficiently navigate the search space by exploiting 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 performance even after running for over $2$M iterations, as shown by their final results in Fig. 4.
The performance gap between LLM-SR and symbolic regression baselines also widens as iterations increase, particularly in the Oscillation 2 and E. coli growth problems, highlighting the effectiveness of LLM-SR’s iterative refinement process.

## 5 Analysis

### 5.1 Ablation Study

Fig. 5 presents an ablation study on the Oscillation 2  problem with GPT-3.5 backbone to investigate the impact of different components in the LLM-SR performance.

Problem Specification.
We create a “w/o Problem Specification” variant by removing the natural language description of the problem and variables of the system from the LLM’s input prompt. In this setup, the model operates solely as an adaptive evolutionary search agent without domain-specific information about the variables and the problem being studied. We find that this variant performs worse, with increases in the Normalized MSE from 2.12e-7 to 4.65e-5 for in-domain data and from 3.81e-5 to 7.10e-3 for OOD data. This highlights the significance of incorporating prior domain knowledge into the equation discovery process.

Figure: Figure 5: Ablation results on the Oscillation 2 problem, showing the impact of problem specification, iterative refinement, parameter optimization, and optimization methods on LLM-SR’s performance.
Refer to caption: extracted/5638465/fig_table_file/ablation7.png

Iterative Refinement.
The “w/o Iterative Refinement” variant is equivalent to the LLM sampling baseline, which generates more than $2$K samples for the initial LLM-SR’s prompt without computing feedback and updating in-context examples over the sampling span. The absence of iterative refinement leads to a substantial performance drop, with the Normalized MSE increasing to 1.01e-1 for in-domain data and 1.81e-1 for OOD data. This emphasizes the importance of the evolutionary search and iterative refinement process in LLM-SR’s success.

Coefficient Optimization.
We also create a “w/o Coeff Optimization” variant that requires the LLM to generate complete equations, including numerical coefficients and constants, in a single end-to-end step. This helps us study the role of parameter optimization using external optimizers in the equation discovery process. We find that this variant shows significantly worse results, with the Normalized MSE increasing to 3.75e-1 for in-domain data and 3.78e-1 for OOD data. These results indicate that the two-stage approach of generating equation skeletons followed by data-driven skeleton parameter optimization is essential for effective performance, as it helps navigate the intricate combinatorial optimization landscape involving both continuous parameter values and discrete equation structures.

Optimization method.
LLM-SR relies on direct and differentiable parameter optimization, a capability not present in current symbolic regression methods. We compare two different optimization frameworks: numpy+BFGS and torch+Adam. The numpy+BFGS variant yielded better-performing equations on both ID (2.12e-7) and OOD (3.81e-5) evaluations compared to torch+Adam (ID: 1.60e-6, OOD: 2.4e-4). However, our initial analysis indicated that this discrepancy could be primarily due to the LLM’s higher proficiency in generating numpy code compared to pytorch, rather than the superiority of one optimization method over the other.
Combining LLM-SR with LLM backbones that are better in generating pytorch code could potentially leverage differentiable parameter optimization to achieve better performance in future works.

### 5.2 Qualitative Analysis

Figure: Figure 6: 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 DSR, uDSR and PySR baselines. Shaded green terms denote recovered symbolic terms from true equations.
Refer to caption: extracted/5638465/fig_table_file/osc_equations-v13.jpg

Fig. 6 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 baseline models.
Moreover, LLM-SR provides explanations and reasoning steps based on the scientific knowledge about the problem, leading to more interpretable terms combined as the final equation.
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 physical characteristics.
In contrast, baselines (DSR, uDSR, PySR) generate equations lacking interpretability and understanding of the physical meanings or relations behind the terms. The equations appear as combination of mathematical operations and variables without clear connection to the problem’s underlying principles.
More detailed results and qualitative analysis on the discovered equations in other benchmark problems and across performance progression curves can also be found in App. E.

## 6 Related Work

LLMs and Optimization.
LLMs have demonstrated remarkable capabilities in various domains, but they are prone to generating incorrect or inconsistent outputs due to either having flawed implicit knowledge or lacking access to external facts (Madaan et al., ; Zhu et al., ). This limitation, often referred to as hallucination, hinders the direct application of LLMs in tasks that require high accuracy.
To mitigate these issues, researchers have explored combining LLMs with feedback or verification mechanisms, such as using the same model for feedback and self-refinement or employing external evaluators (Madaan et al., ; Yang et al., ; Haluptzok et al., ). Recent works have also coupled LLMs with evaluators in an iterative optimization loop, where LLMs act as evolutionary search agents (Lehman et al., ; Liu et al., ; Wu et al., ; Lange et al., ). LLMs leverage their prior knowledge, reasoning, and generation capabilities to perform adaptive mutation and crossover operations (Meyerson et al., ). This integration has shown promising results in applications like prompt optimization (Yang et al., ; Guo et al., ), neural architecture search (Chen et al., ; Zheng et al., ), and heuristic discovery (Romera-Paredes et al., ).
Most related to our work is FunSearch (Romera-Paredes et al., ) that combines LLMs with systematic evaluators to search for programs that push the boundaries in solving some established open mathematical problems.
Building upon these ideas, our LLM-SR framework employs LLM as an informed optimizer, leveraging its scientific prior knowledge and data-driven evaluators to discover mathematical equations underlying scientific observations.

LLMs for Scientific Discovery.
The use of LLMs in scientific discovery has gained significant attention in recent years (Wang et al., ). Researchers have explored the potential of LLMs across various scientific domains, such as drug discovery, biology, and materials design (AI4Science & Quantum, ). Recently, the ability of LLMs to propose plausible scientific hypotheses by leveraging their existing knowledge and their potential for data-driven discovery has been a topic of growing interest (Majumder et al., ; Zheng et al., ; Qi et al., ).
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 addresses this gap by introducing a novel approach that harnesses the scientific prior knowledge within LLMs and integrates it with data-driven evaluation. This approach showcases the potential of LLMs in scientific equation discovery, contributing to the rapidly growing body of research on the use of LLMs for scientific advancement.

Symbolic Regression.
Symbolic regression (SR) methods can be broadly categorized into search-based approaches, learning-based models, and hybrid methods. Search-based approaches mainly explore the space of mathematical expressions using evolutionary algorithms or reinforcement learning (Schmidt & Lipson, ; Cranmer, ; Petersen et al., ; Sun et al., ). Learning-based models, on the other hand, leverage large-scale synthetic data and Transformers to learn the mapping between numeric input observations and output mathematical expressions (Biggio et al., ; Kamienny et al., ). 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 (Shojaee et al., ; Udrescu & Tegmark, ; Mundhenk et al., ; Meidani et al., ). 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. Some works have attempted to address these issues by incorporating prior physical knowledge such as dimension analysis (Udrescu & Tegmark, ; Tenachi et al., ; Meidani & Barati Farimani, ) or using declarative bias and structures with pre-defined grammars
(Todorovski & Dzeroski, ; Todorovski & Džeroski, ). However, 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.

## 7 Discussion and Conclusion

In this work, we introduced LLM-SR, a novel framework that leverages the power of Large Language Models (LLMs) for scientific equation discovery.
We demonstrated the effectiveness of LLM-SR by evaluating it on diverse set of scientific problems, showcasing its ability to uncover physically meaningful and accurate equations.
Despite the promising results, LLM-SR has limitations that should be acknowledged. The computational resources associated with loading LLMs and generating a large number of equation programs can be substantial, especially when dealing with complex problems or large-scale datasets. Additionally, LLM-SR’s reliance on the scientific prior knowledge embedded within the LLMs may also be limited, biased, or incorrect in certain domains, potentially affecting the quality of the discovered equations.
However, we believe that the potential of using LLM-SR framework can go beyond what has been studied in this work, and that LLM-enabled scientific equation discovery can facilitate this task in various fields of science and engineering.
Moreover, with the rapid pace of LLM improvements, the computational cost is expected to decrease, making it more accessible and efficient.

Future research efforts can focus on several key areas to further enhance the capabilities of LLM-SR. Exploring the integration of more powerful or domain-specific language models can potentially improve the quality and relevance of the generated equations. Incorporating retrieval-augmented learning techniques, where the LLMs are augmented with literature-driven scientific knowledge, can also provide additional context and guidance during the equation discovery process. Developing methods to verify and guide the model towards generating equations that are consistent with established scientific principles can help ensure the physical validity of the discovered equations. Another potential future work is the development of more comprehensive benchmarks for evaluating LLM-based equation discovery methods. These benchmarks should either include a potentially empirical equation discovery process or be carefully designed to simulate such a process. This approach ensures that the equations generated by the models are not simply memorized but are derived through a genuine discovery process.
As LLMs continue to evolve and become more powerful, LLM-SR has the potential to become an indispensable tool for researchers and scientists, accelerating scientific discovery and innovation across various domains.
With the development of better benchmarks and evaluation protocols, we can unlock the full potential of LLM-based equation discovery and pave the way for accelerating scientific advancements.

## Appendix

## Appendix A Implementation Details

### A.1 SR Baselines

We compare LLM-SR against several state-of-the-art Symbolic Regression (SR) baselines, including GPlearn, Deep Symbolic Regression (DSR) (Petersen et al., ), Unified Deep Symbolic Regression (uDSR) (Landajuela et al., ), and PySR (Cranmer, ). GPlearn is a pioneering and standard GP-based SR approach that uses the open-source gplearn(^3^33[https://gplearn.readthedocs.io/en/stable/](https://gplearn.readthedocs.io/en/stable/)) package with parameters population size: $1000$, tournament size: $20$, and maximum generations: $2$M. DSR employs a RNN-based reinforcement learning search over symbolic expressions. uDSR extends DSR by incorporating large-scale pretraining and GP search at the decoding stage. Both DSR and uDSR approaches are implemented using deep-symbolic-optimization (DSO) (^4^44[https://github.com/dso-org/deep-symbolic-optimization](https://github.com/dso-org/deep-symbolic-optimization)) package with the standard default parameters: learning rate of $0.0005$, batch size of $500$, and a maximum of $2$M iterations.
PySR is also used with the state-of-the-art open-source pysr(^5^55[https://github.com/MilesCranmer/PySR](https://github.com/MilesCranmer/PySR)) package for SR that uses asynchronous GP-based evolution, implemented with a population size of $4,000$, up to $100$ sub-populations, and a maximum of $2$M iterations.
For a fair comparison, we allowed the SR baselines to run for over $2$M iterations until convergence to their best performance, while LLM-SR variants ran for only around $2$K iterations.
Each of SR basline models was subjected to a long execution time $5$ times, each time running for over $2$ million iterations. This extensive number of evaluations and search iterations ensures robust evaluation of the models’ capability to converge towards optimal solutions and effectively explore the vast equation space of each problem.

### A.2 LLM-SR

The LLM-SR framework was evaluated using two distinct backbone LLM models: GPT-3.5-turbo and Mixtral-8x7B-Instruct. For both variants, the framework underwent around $2.5$K iterations, a comparatively lower number of search iterations compared to symbolic regression baselines (around $2$M) due to the enhanced efficiency of this approach because of embedded scientific prior knowledge. The best results obtained from these iterations were then documented and reported.
Key parameters of the LLM-SR framework include:
Sampling Per Prompt: LLM-SR generates $b=4$ samples with temperature $\tau=0.8$ per prompt; Program Database Structure: The framework incorporates a program database consisting of $m=10$ islands. This database structure facilitates diverse program generation and efficient sampling.
Prompt Augmentation: To refine the prompt’s effectiveness, it includes $k=2$ best-shot in-context examples as experience demonstration sampled from the experience buffer. This strategy leverages the buffer’s accumulated knowledge to enhance prompt relevance and model performance over evolutionary mutations and the search process.
Execution Timeout: A critical operational parameter in our setup was the execution timeout set at $30$ seconds. If the execution of a program generated by the LLM exceeded this limit, the process was halted, and the evaluation score of None was returned. This constraint ensured timely progress and resource efficiency in the evaluation process.
Parallel Evaluation: In this framework, we deployed $e=4$ evaluators to operate concurrently. This parallelization allowed for rapid and efficient assessment of the generated programs per prompt.
As pointed out, for the LLM backbone models, a sampling temperature of $\tau=0.8$ was utilized.
This temperature setting was chosen to balance creativity (exploration) and adherence to the problem constraints and reliance on the prior knowledge (exploitation), based on preliminary experiments.

## Appendix B Additional Details of LLM-SR method

##### Hypothesis Generation and Data-driven Evaluation

Fig. 2 provided an example of specification for Bacterial growth problem. Here, Fig. 7 showcases illustrative examples of prompts and specifications tailored for the Oscillation 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 knowledge about the physical meaning and relations of variables to generate scientifically plausible hypotheses for new equation programs.

Figure: Figure 7: Example of input prompts program body for (a) Oscillation and (b) Stress-Strain problems, with problem specification, and the initial input equation program example (set as simple linear equation skeleton). For better readability, the evaluation function details are not included in this figure.
Refer to caption: extracted/5638465/fig_table_file/spec_prompts3.png

Fig. 8 also shows an example of prompt and specification for LLM-SR that prompts the model to 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.

During each prompting step, the language model generates $b=4$ distinct equation program skeletons using a generation temperature of $\tau=0.8$. These equation skeleton programs are concurrently evaluated to gather feedback. Evaluation is constrained by time and memory limits set at $T=30$ seconds and $M=2$GB , respectively. Equation programs that exceed these limits are disqualified and considered as discarded hypotheses by returning None scores.

Figure: Figure 8: An example of prompt structure, containing problem specification, evaluation and optimization function, and equation program with pytorch and Adam Optimizer
Refer to caption: extracted/5638465/fig_table_file/spec_torch.jpg

Experience Buffer Management.
The pairs of equation skeleton hypotheses and their corresponding scores, denoted as $(f,s)$, are then stored in an experience buffer $\mathcal{P}_{t}$ for further refinement in subsequent iterations of the search process.
To encourage diversity and avoid getting stuck in local optima, we follow (Cranmer, ) and adopt an islands model, also known as a multiple population or multiple-deme model for managing the experience buffer (Cranmer, ).
We initialize $m$ separate islands, each containing a copy of the equation program example provided in the initial prompt (equation$_$v0 in Fig. 2). These islands evolve independently, allowing for parallel exploration of different regions in the equation program space. At each iteration $t$, the newly generated equation program hypotheses $\mathcal{F}_{t}$ and their corresponding fitness scores $\{(f,s):f\in\mathcal{F}_{t},s=\text{MSE}(f,\mathcal{D})\}$ are added to the same island from which the prompts were sampled, if the score of new hypothesis $s$ is better than the current best-score in the island.
To maintain the quality and diversity of the experience buffer, we periodically reset the worst-performing islands. Every $T_{reset}$ iterations (e.g., every $4$ hrs), we identify the $m/2$ 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-performing equation program, obtained by randomly selecting one of the surviving $m/2$ islands and copying its highest-scoring equation program (favoring older programs in case of ties). This reset mechanism allows the framework to discard stagnant or unproductive regions of the equation program space and focus on more promising areas.
Within each island, we further cluster the equation programs based on their signature, which is defined as the equation program score. Equation programs with identical signatures are 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 the 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. Here, we use a two-stage sampling process. First, we randomly select an island from the $m$ available islands. Then, within the selected island, we sample $k$ equation programs (typically, $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 a cluster based on its evaluation score, favoring clusters with higher scores (i.e., higher-quality equation programs). Let $s_{i}$ 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_{i}$ of choosing cluster $i$ is given by:

$$ $P_{i}=\frac{\exp\left(\frac{s_{i}}{\tau_{c}}\right)}{\sum_{i^{\prime}}\exp \left(\frac{s_{i^{\prime}}}{\tau_{c}}\right)},\quad\tau_{c}=T_{0}\left(1-\frac {u\bmod N}{N}\right)$ $$

,
where $\tau_{c}$ is the temperature parameter, $u$ is the current number of equation programs in the island, and $T_{0}$ and $N$ are hyperparameters. This selection approach is known as the Boltzmann selection procedure (Maza & Tidor, ). Once a cluster is selected, we sample an equation program within that cluster, favoring shorter programs. Let $l_{i}$ denote the negative length of the $i$-th program within the chosen cluster (measured as the number of characters), and let $\tilde{l}_{i}=\frac{l_{i}-\min_{i^{\prime}}{l_{i^{\prime}}}}{\max_{i^{\prime}}
{l_{i^{\prime}}}+10^{-6}}$. We set the probability of selecting each equation program proportional to $\exp\left(\tilde{l}_{i}/\tau_{p}\right)$, where $\tau_{p}$ is a 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 a 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.

## Appendix C Limitation of Feynman Benchmark Problems

Figure: Figure 9: Trajectory of Normalized MSE score over iterations for LLM-SR (GPT-3.5) on Feynman benchmark equations versus new benchmark problems
Refer to caption: extracted/5638465/fig_table_file/feynman_curve4.png

The Feynman benchmark dataset (Udrescu & Tegmark, ), while commonly used for evaluating symbolic regression and equation discovery methods, may not effectively capture the true discovery process and reasoning capabilities of LLMs.
Fig.9 compares the performance of LLM-SR with GPT-3.5 backbone across different benchmark problems, illustrating the best score trajectory of Normalized Mean Squared Error against the number of iterations for each problem. For the Feynman benchmark equations, LLM-SR rapidly achieves low MSE scores within very few iterations, suggesting that LLMs may have memorized or become overly familiar with these equations due to their prevalence in training data. Qualitative examples in Fig. 10, 11, and 12 further support this notion, as the LLM’s one-pass responses to some Feynman problems not only exhibit functional accuracy but also mirror the exact form of the corresponding physics expressions.

Figure: Figure 10: An example of LLM response to Feynman I.37.4 problem, demonstrating model memorization without iterative search.
Refer to caption: extracted/5638465/fig_table_file/qual_feyn_I37_4_v3.jpg

Figure: Figure 11: An example of LLM response to Feynman I.26.2 problem, demonstrating model memorization without iterative search.
Refer to caption: extracted/5638465/fig_table_file/qual_feyn_I26_2_v3.jpg

Figure: Figure 12: An example of LLM response to Feynman I.12.11 problem, demonstrating model memorization without iterative search.
Refer to caption: extracted/5638465/fig_table_file/qual_feyn_I12_11_v3.jpg

In contrast, the newly designed benchmark problems, which challenge the model with less familiar components, require more reasoning and refinement steps for LLM-SR to discover well-fitting equations. The increased number of iterations in the refinement process demonstrates that these new problems effectively simulate the discovery process rather than simply relying on retrieval from memory. This observation highlights the necessity for introducing novel benchmark problems that can better assess the true reasoning capabilities of LLMs in equation discovery while leveraging its scientific prior knowledge, as the commonly used Feynman equations may not provide an accurate evaluation due to the models’ over-exposure and memorization of these exact physic equations.

## Appendix D Additional Details on Benchmark Problems and Datasets

### D.1 Nonlinear Oscillator Equations

**Table 2: Parameter values for Oscillator datasets.**
| Dataset | Time range | initial values | F | $\alpha$ | $\beta$ | $\delta$ | $\gamma$ | $\omega$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Oscillator 1 | (0, 50) | {x=0.5, v=0.5} | 0.8 | 0.5 | 0.2 | _ | 0.5 | 1.0 |
| Oscillator 2 | (0, 50) | {x=0.5, v=0.5} | 0.3 | 0.5 | 1.0 | 5.0 | 0.5 | _ |

Figure: Figure 13: Phase diagrams of trajectories corresponding to custom oscillators: (a) Oscillator 1 and (b) Oscillator 2
Refer to caption: extracted/5638465/fig_table_file/oscillation_phase_v2.jpg

In this work, we simulate two nonlinear oscillators using ’solve_ivp’(^6^66[https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)) in scipy library to generate data. The parameters and initial values of these simulations are provided in Table 2, and the equations used to generate data for this task are repeated here for ease of reference:

$$ $\dot{v}=F\sin(\omega x)-\alpha v^{3}-\beta x^{3}-\gamma x\cdot v-x\cos(x)$ $$

and

$$ $\dot{v}=F\sin(\omega t)-\alpha v^{3}-\beta x\cdot v-\delta x\cdot\exp(\gamma x).$ $$

Fig. 13 illustrates the phase plane diagrams of these nonlinear damped oscillators. These oscillators have nonlinear terms for driving force, restoring force, and damping force, making them challenging problems for equation identification.
To assess the generalization capability of the predicted equations, we partition the data into training, in-domain validation, and out-of-domain validation sets based on the trajectory time in simulations. Specifically, we utilize the time interval $T=[0,20)$ to evaluate the out-of-domain generalization of the equations.

### D.2 E. coli Growth Rate Equations

The custom equation used to generate data for this task is:

$$ $\frac{dB}{dt}=\mu_{max}B\left(\frac{S}{K_{s}+S}\right)\frac{\tanh\left(k(T-x_{ 0})\right)}{1+c(T-x_{decay})^{4}}\exp\left(-\left|\text{pH}-\text{pH}_{opt} \right|\right)\sin\left(\frac{(\text{pH}-\text{pH}_{min})\pi}{\text{pH}_{max}- \text{pH}_{min}}\right)^{2},$ $$

where temperature and pH dependency are approximated by complex nonlinear functions with new unseen operators, the effect of population density is assumed to be linear (single bacterial population), and $\left(\frac{S}{K_{s}+S}\right)$ denotes Monod equation (Monod, ) for substrate concentration model. Fig. 14 shows the numeric behavior of custom designed temperature and pH functions, $f_{T}(T)$ and $f_{pH}(pH)$, along with some previous forms studied in literature (Rosso et al., ).

Figure: Figure 14: Scheme of some common temperature and pH models in bacterial growth, as well as the custom model behavior.
Refer to caption: extracted/5638465/fig_table_file/pH_dependency_v2.jpg

### D.3 Material Stress Behavior Analysis

Figure: Figure 15: Stress-strain curves of Aluminium 6061-T651 under various temperatures (data from (Aakash et al., ))
Refer to caption: extracted/5638465/fig_table_file/stress_gt_v2.png

Fig. 15 depicts the stress-strain curves of Aluminium 6061-T651 at different temperatures, showing the elastic, plastic, and failure regions of material and highlighting the significant influence of temperature on the material’s response.
The stress-strain relationship for aluminum alloys can be described using various constitutive models (Samantaray et al., ), such as the Johnson-Cook model (Johnson & Cook, ) which incorporates the effects of strain, strain rate, and temperature on the stress. A simplified version of J-C equation (by neglecting the effect of strain rate) is:

$$ $\sigma=\left(A+B\varepsilon^{n}\right)\left(1-\left(\frac{T-T_{r}}{T_{m}-T_{r} }\right)^{m}\right),$ $$

where $\sigma$ is the stress, $\varepsilon$ is strain, $T$ is the current temperature, $T_{r}$ is the reference temperature, $T_{m}$ is the melting temperature, and $A$, $B$, $C$, $n$, $m$ are material constants.

For this problem, the data represents the tensile behavior of material under uniaxial tension for 6 different temperatures. We allocate the data corresponding to $T=200^{\circ}$C for use as the validation set.

## Appendix E Further Results and Visualizations

##### Performance 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 given the best-shot sampled in-context equation examples provided as experience demonstration sampled from the experience buffer.

Figure: Figure 16: Performance Trajectory of LLM-SR along with the best-scoring found equation programs over iterations on the Oscillation 2 problem.
Refer to caption: extracted/5638465/fig_table_file/osc2_annotated_curve_v2.jpg

Fig. 16 shows the NMSE values for the Oscillation 2 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 while ground truth equation for this dataset is
$\dot{v}=0.3\sin(t)-0.5v^{3}-x\cdot v-5.0x\cdot\exp(0.5x),$
LLM-SR has discovered the equation $\dot{v}=0.3\sin(t)-0.5v^{3}-x\cdot v+5.0(1-\exp(x)).$
By evaluating the different terms in these two forms, we observe that in fact $5.0(1-\exp(x))\approx-5.0x\cdot\exp(0.5x)$ for $x\in(-2,2)$, which is the approximate range of displacement in this dataset.

Figure: Figure 17: Performance Trajectory of LLM-SR along with the best-scoring found equation programs over iterations on the Oscillation 1 problem.
Refer to caption: extracted/5638465/fig_table_file/osc1_annotated_curve_v2.jpg

Fig. 17 illustrates the evolution of equations for the Oscillation 1 problem. It is evident that the model attempts to incorporate various nonlinear terms corresponding to driving, restoring, and damping forces, as evidenced by comments or variable names within the code, aiming to enhance accuracy. An intriguing observation emerges wherein the model identifies a trivial solution for the nonlinear oscillation problem, exploiting the flexibility in its code representations. As depicted in Fig.  17, the last step discovered equation yielding very low error turns out to be a trivial solution based on the physical interpretation of variables in this context. This solution was discovered through the utilization of the np.gradient function in the numpy library, employing numerical differentiation of $0.5v^{2}$ with respect to displacement ($x$), leading to an approximation of acceleration: $-6.67\times\frac{d}{dx}(-0.076v^{2})=0.5\frac{d}{dx}(v^{2})=\frac{dv}{dt}$. Discovery of this edge case by at last iterations is a very interesting observation

Similarly, Fig.  18 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 ($pH_opt$ and $T_opt$) from the early iterations which comes from model prior knowledge about the bell-shaped effect of these variables on the growth rate. To enhance accuracy, the model necessitates exploration and incorporation of 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 $\frac{S}{K+S}$. Additionally, the figure demonstrates LLM-SR’s comprehension that different components of the function should be multiplied together in the final step, underscoring how prior knowledge of the problem structure can guide LLM-SR’s evolutionary steps.

Figure: Figure 18: Performance Trajectory of LLM-SR along with the best-scoring found equation programs over iterations on the E. coli growth problem.
Refer to caption: extracted/5638465/fig_table_file/bg22_annotated_curve_v2.jpg

Fig. 19 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. An additional significant observation for this problem is that stress-strain relationships often exhibit piece-wise behavior (as it can also be observed in Fig. 15 representing data obesrvations), which traditional symbolic regression models struggle to identify. 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 sigmoid function to model piece-wise relations. This differentiability and smooth approximation of if-else conditions are particularly helpful for the parameter optimization step, providing smooth functions for the optimizer to navigate.

Figure: Figure 19: Performance Trajectory of LLM-SR along with the best-scoring found equation programs over iterations on the Stress-Strain problem.
Refer to caption: extracted/5638465/fig_table_file/stress_annotated_curve_v2.jpg

##### Qualitative Analysis

Fig.  20 and Fig. 21 depict the equation programs identified by LLM-SR and other Symbolic Regression (SR) methods for the E. coli growth and the stress-strain problems, respectively. The diverse range of equation forms identified by different SR baselines reflects the challenges posed by the datasets. Notably, in both datasets, the SR methods yield either lengthy or highly nonlinear equations that deviate from the prior knowledge of the systems, as evidenced by the poor out-of-domain (OOD) performance scores in Table 1. In contrast, LLM-SR finds flexible equation programs that are more interpretable and aligned with the domain-specific scientific priors of the systems.

Figure: Figure 20: Discovered equations from LLM-SR and other SR baseline methods for E. coli bacterial growth rate problem.
Refer to caption: extracted/5638465/fig_table_file/bacterial_equations-v5.jpg

Figure: Figure 21: Discovered equations from LLM-SR and other SR baseline methods for Stress-Strain problem.
Refer to caption: extracted/5638465/fig_table_file/ss_equations-v4.jpg

Fig. 22 and Fig. 23 offer a 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 exhibit a notable degree of alignment with both the in-domain and out-of-domain regions of these systems. This alignment suggests that LLM-SR effectively captures the underlying patterns and dynamics of the data, enabling it to better generalize to unseen data points beyond the training domain.

Figure: Figure 22: Qualitative evaluation of the performance of LLM-SR on Oscillation 2 problem compared to uDSR and PySR baselines. Plots show the target acceleration $\ddot{x}$ with respect to time ($t$), displacement ($x$), and velocity ($v$).
Refer to caption: extracted/5638465/fig_table_file/qual_osc2_v2.jpg

Figure: Figure 23: Qualitative evaluation of LLM-SR performance for Stress-Strain problem compared to uDSR and PySR baselines.
Refer to caption: extracted/5638465/fig_table_file/qual_stress_v2.jpg

## Peer Reviews

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
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.

2. Providing a real-world example of how LLM-SR assists in scientific problem solving could make the results more compelling.
soundness: 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.

2. Robust evaluation: The authors have carefully designed benchmark problems that align well with the real-world applications.

3. 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.

4. Well written: The paper is well written with clear data presentation.
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.
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.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. **Clarification on the Necessity of Complex Mechanisms**:
   - *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?
   - *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.

2. **Verification of LLM Prior Knowledge Utilization**:
   - *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?
   - *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.

3. **Evaluation on Real-World Datasets with Noise and Outliers**:
   - *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?
   - *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.

4. **Impact of LLM Programming Performance on Results**:
   - *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?
   - *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.

5. **Interpretability of Results**:
   - *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?
   - *Suggestion*: Presenting case studies or examples where the equations discovered by LLM-SR provide clearer or more scientifically insightful interpretations would reinforce this claim.

6. **Broader Comparison with Hybrid Methods**:
   - *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?
   - *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.
soundness: 4
strengths: 1. **Innovative Application of LLMs for Scientific Equation Discovery**:
   - 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.

2. **Experience Buffer and Iterative Optimization**:
   - 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.

3. **Comprehensive Multi-Domain Benchmarking**:
   - 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.

4. **Promising Experimental Results**:
   - 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.
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.

**Main Contributions:**
1. **Novel Framework**: LLM-SR combines LLMs' strengths with evolutionary algorithms to navigate complex equation discovery.
2. **Integration of Scientific Priors**: The method leverages LLMs' embedded scientific knowledge for more efficient hypothesis generation.
3. **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.
4. **Performance Advantage**: LLM-SR reportedly outperforms traditional symbolic regression methods, showing better accuracy and generalization, especially in out-of-domain tests.
5. **Ablation Study**: Demonstrates the importance of various components, such as problem specification and iterative refinement, in the overall performance of LLM-SR.
weaknesses: **Complexity and Necessity of the Method**:
The 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.

**Utilization and Verification of Prior Knowledge**:
Although 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.

**Breadth and Realism of Experimental Design**:
The 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.

**Impact of Programming Capabilities on Experimental Results**:
The 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.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. L243 briefly described how population initialization is performed. Is the linear equation skeleton manually written? How many skeletons are provided at initialization?  
2. In L149, should the expectation be over the data points in $\mathcal{D}$?  
3. 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).
soundness: 3
strengths: 1. To the best of my knowledge, framing SR as equation _template_ discovery seems novel and a useful contribution.
2. The analyses presented to demonstrate LLM memorization of the Feynman dataset are quite convincing.
3. The new benchmark problems are well-motivated and look like sound contributions.
4. The paper is, overall, clearly written and easy to follow, and visualizations are used well in various places.
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.
weaknesses: My main concerns with the work are twofold:
1. 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.  
2. 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.  
  (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.  
  (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.  

Minor:
1. 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.  
2. It was mildly odd to read that mathematical equations have been "unreasonably" effective in describing complex phenomena. Unclear why we think it is unreasonable.  

[1] Meyerson et al., 2023. Language model crossover: Variation through few-shot prompting.  
[2] Romera-Paredes et al., 2023. Mathematical discoveries from program search with large language models.  
[3] Yang et al., 2024. Large Language Models as Optimizers.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. How would LLM-SR perform with even larger LLMs or specialized scientific models?
2. How does LLM-SR perform with complex systems involving interactions across multiple equations?
soundness: 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.
2. LLM-SR shows strong performance not only in-domain but also in out-of-domain scenarios, which is often challenging for symbolic regression techniques.
3. 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.
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.
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.
2. 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.

Meta-review: This paper proposes LLM-SR, a framework that leverages large language models for symbolic regression by integrating program synthesis, numerical optimization, and evolutionary search to discover accurate and generalizable scientific equations. Most of the reviewers agreed that the proposed approach is novel and technically solid. Moreover, this paper is written with clear data presentation and effective use of visualizations. 

During the review process, Reviewer 89Tg raised ethical concerns regarding dataset usage. The authors addressed this in their rebuttal by explaining the licensing and availability of the datasets, with detailed information provided in the paper's appendix. The reviewer confirmed that their concerns had been addressed and suggested the paper for the best paper award. As a result, I believe further ethical review is unnecessary. Nevertheless, I strongly recommend that the authors further discuss dataset usage in the final version to enhance clarity and transparency.

All the reviews tend to accept the paper during the reviewer-author discussion phase. Therefore, I recommend this paper to the ICLR 2025 conference.

Decision: Accept (Oral)