Knowledge integration for physics-informed symbolic regression using pre-trained large language models

TL;DR

This paper shows how pre-trained large language models can help integrate domain knowledge into symbolic regression, making it easier and more effective for scientific discovery.

Contribution

The novel approach integrates LLMs into symbolic regression's loss function to automate domain knowledge incorporation.

Findings

LLM integration consistently improves the reconstruction of physical dynamics from data.

More informative prompts significantly enhance the performance of the method.

The approach works across multiple SR algorithms and physical systems.

Abstract

Symbolic regression (SR) has emerged as a powerful tool for automated scientific discovery, enabling the derivation of governing equations from experimental data. A growing body of work illustrates the promise of integrating domain knowledge into the SR to improve the discovered equation’s generality and usefulness. Physics-informed SR (PiSR) addresses this by incorporating domain knowledge, but current methods often re- quire specialized formulations and manual feature engineering, limiting their adaptability only to domain experts. In this study, we leverage pre-trained Large Language Models (LLMs) to facilitate knowledge integration in PiSR. By harnessing the contextual understanding of LLMs trained on vast scientific literature, we aim to automate the incorporation of domain knowledge, reducing the need for manual intervention and making the process more accessible to a broader range of scientific problems. Namely, the LLM is integrated into the SR’s loss function, adding a term of the LLM’s evaluation of the SR’s produced equation. We extensively evaluate our method using three SR algorithms (DEAP, gplearn, and PySR) and three pre-trained LLMs (Falcon, Mistral, and LLama 2) across three physical dynamics (dropping ball, simple harmonic motion, and electromagnetic wave). The results demonstrate that LLM integration consistently improves the reconstruction of physical dynamics from data, enhancing the robustness of SR models to noise and complexity. We further explore the impact of prompt engineering, finding that more informative prompts significantly improve performance.

The online version contains supplementary material available at 10.1038/s41598-026-35327-6.

Figures

A schematic view of the Large Language Model-based Knowledge Integration into the Symbolic Re- gression’s searching (and optimization) process.

A schematic view of the experimental setup.

Full text

Introduction

Scientific discovery is based on the interplay between empirical observation and theoretical modeling^1,2^. Throughout the subsequent centuries, this approach of expressing scientific knowledge through mathematical relationships has remained central to physics and other natural sciences^3–7^, mathematical formulations have consistently provided the language through which we articulate physical laws^8^. These symbolic expressions serve not merely as descriptive tools but as engines of discovery^9^, enabling scientists to derive new insights and make predictions beyond what has been directly observed.

The process of deriving mathematical expressions from experimental data, formally known as symbolic regression (SR), has traditionally been a labor-intensive endeavor requiring deep domain expertise and mathematical intuition^10–12^. Scientists have historically relied on a combination of theoretical understanding, creative insight, and iterative refinement to formulate equations that accurately capture observed phenomena^13,14^. In recent years, the emergence of machine learning (ML) and artificial intelligence (AI) has significantly accelerated the process of scientific knowledge extraction, particularly in domains that require complex pattern recognition and hypothesis generation^15–18^. These advantages occurred hand in hand with the exponential increase in available scientific data across all fields, including medicine^19^, economics^20^, and physics^21^, to name a few. Despite their established benefits^22,23^, current AI, in general, and ML models, in particular, are not explainable and lack the reproducibility and generalizability associated with symbolic mathematical expressions^24,25^.

To address this gap, two main branches have emerged, one aims to make AI models more explainable^26–28^, while the other tries to automate the SR process using computational methods^29^. Unlike traditional regression techniques, SR does not assume a predefined model structure but instead searches for the most appropriate mathematical representation using optimization techniques^30^. In physics, this process is particularly relevant as it allows for the derivation of governing equations from raw experimental or simulated data, offering a data-centric approach to understanding physical systems^31^. Indeed, SR has proven to be a valuable tool in this endeavor, capable of rediscovering known physical laws, such as Newton’s laws of motion and Maxwell’s equations from observational data^32,33^. However, despite its potential, SR faces significant challenges, particularly in ensuring that the discovered expressions remain interpretable and physically meaningful rather than simply fitting a given experimental data^34–36^.

In order to address these challenges, recent studies have explored physics-informed SR (PiSR), which incorporates domain knowledge to guide the discovery process^37^. These approaches leverage physical constraints, symmetries, and conservation laws to improve the quality of the inferred models. By embedding these principles into the learning process, PiSR enhances the plausibility and generalizability of discovered equations^38,39^. However, the integration of domain knowledge into PiSR remains highly technical, often requiring specialized formulations, manual feature engineering, and handcrafted constraints, making it difficult to generalize across different physical domains^37^.

In this study, we propose a novel approach by leveraging pre-trained Large Language Models (LLMs), a family of deep learning models trained on vast textual and scientific datasets to develop a broad understanding of language (enabling them to generate, analyze, and integrate domain-specific knowledge in various tasks), to facilitate knowledge integration in PiSR. LLMs, trained on vast corpora of scientific literature and equations, have demonstrated remarkable capabilities in understanding and generating scientific content^40^. By harnessing the contextual understanding of LLMs, we aim to automate the incorporation of domain knowledge in SR, reducing the need for manual intervention and making the process more accessible to a broader range of scientific problems. Figure 1 shows a schematic view of the LLM-based knowledge integration into SR.Fig. 1A schematic view of the Large Language Model-based Knowledge Integration into the Symbolic Re- gression’s searching (and optimization) process.

The rest of the manuscript is organized as follows. Section “Related Work”presents the current state-of-the-art in SR and LLMs, as well as previous attempts to integrate the two. Next, section “Symbolic regression” formally introduces the proposed LLM and SR integration mechanism. Following that, section “Physics-informed symbolic regression” describes the experimental setup used to evaluate the proposed method. Afterward, section “Large language models and symbolic regression” outlines the obtained results. Finally, Section “LLM-based Knowledge Integration to SR” discuss the results in an applied context and suggests possible future work.

Related work

This section provides an overview of research related to automatic scientific discovery, SR, and LLMs. We initially discussed the current state-of-the-art SR with a particular focus on PiSR, followed by a review of recent advancements in large language models, and finally we explore the integration of LLMs and SR models.

Symbolic regression

SR is a modeling technique for discovering mathematical expressions that represent the relationship between a set of input variables and an output variable^41–43^. SR does not require a predefined model structure, unlike other regression models, such as linear regression, as it searches through a large number of potential mathematical expressions to identify the function that best fits the given data^44^. This approach is particularly useful in cases where the underlying functional form is unknown or complex^45–47^.

Roughly speaking, one can categorize the SR methods into four main groups based on the underlying computational techniques: brute force, sparse regression, deep learning, and genetic algorithms^48^. The brute-force-based SR methods can theoretically solves any SR tasks. Nevertheless, practically, their computational cost is unrealistic, and they are extremely prone to overfitting^49^. The method tests all possible equations to find the one that performs optimally, differing in the way they cover the search space^50^. The sparse regression SR method, which identifies parsimonious models with the help of optimization that promotes sparsity, significantly reduces the search space. Among them, SINDy is designed for scientific use cases, which employs the Lasso linear model for sparse identification of nonlinear dynamical systems behind time series data^51^. SINDy iterates between partial least squares fitting and thresholding (promoting sparsity) steps^52^. Deep learning SR methods excel in handling noisy data due to the high outlier resistance of neural networks, but are limited in their generability^53^. These methods use neural networks to produce analytical equations. For instance, Petersen et al. (2019) pro- posed a Deep Symbolic Regression (DSR) system for general SR tasks, utilizing reinforcement learning to train a generative recurrent neural network (RNN) model of symbolic expressions and employs a variant of the Monte Carlo policy gradient technique called “risk-seeking policy gradient” to adapt the generative model to the exact formulation^54^. Lastly, genetic-algorithm-based SR methods treat mathematical expressions as individuals (or genes) in a population^48^. These individuals evolve over time through mechanisms such as selection, crossover, and mutation, gradually refining the equations that can better fit the data. This evolutionary approach allows SR to discover interpretable models in a flexible way, without being constrained by predefined assumptions^55,56^. For example, the gplearn Python library implements a genetic algorithm for SR, showing a promising ability to re-discover known physical equations^57^. It first constructs a population of stochastic formulas representing the relationship between known independent variables (features) and dependent variables (objectives), presented in a tree structure. These subtrees are then replaced and restructured in a stochastic optimization process that computes fitness by executing the tree and evaluating its output, using a stochastic strategy of survival of the fittest.

In recent years, many SR algorithms have been proposed^58–60^. From these, several SR algorithms gain much popularity and provide a simple programing interface to use, including the DEAP (Distributed Evolutionary Algorithm in Python)^61,62^, gplearn^63,64^, and PySR^65^. All three algorithms belong to the genetic- algorithm-based SR family of methods.

Physics-informed symbolic regression

In a general sense, a learning algorithm becomes physics-informed when physical knowledge is introduced as observational, inductive, or learning biases that steer learning toward physically consistent solutions^66^. In particular, physics guidance can be implemented by designing loss functions, constraints, or inference procedures that explicitly favor solutions adhering to known physics^67,68^.

SR is widely adopted for physical problems, which often face a multivariate and noisy experimental data from nonlinear systems and require analytical equations to discover general laws of nature^69–71^. In this context, simply using SR methods provides limited results as the SR needs to make sure the prediction equation follows scale consistency along with other properties^48^. To this end, Physics-Informed Symbolic Regression (PiSR) has emerged, aiming to close this gap^48,71,72^. PiSR is designed to guide the search process of symbolic expressions by incorporating physics principles such as conservation laws, symmetry, and sparsity assumptions, to name a few^71^.

For example, the AI Feynman SR algorithm combines physics-inspired techniques with machine learning and deep learning methods to improve the accuracy of physics discoveries, showing promising results across multiples sub-field of physics^35^. Similarly, Martius and Lampert (2016) proposed Equation Learner (EQL) which is based on a multilayer feedforward neural network designed to learn the dynamical equations of a physical system and to be able to extrapolate to unseen domains^73^. It allows efficient training based on gradients, and its architecture consists of linear mappings, nonlinear transformations consisting of specific unitary and binary (multiplicative) units, and a linear readout layer. Golden (2024) proposed the Scalable Pruning for Rapid Identification of Null vecTors (SPRINT) algorithm, which is a model-independent sparse regression algorithm that extends the iterative singular value decomposition based exhaustive search algorithm^74^.

Large language models and symbolic regression

Large language models (LLMs) have recently gained significant traction, often surpassing traditional approaches in multiple tasks in information retrieval, information representation, and natural language processing^75–80^. Research combining LLMs with established techniques such as collaborative filtering and content- based methods has yielded highly accurate recommendations^81,82^. Furthermore, LLMs are proving valuable in improving specific recommendation challenges, including sequential recommendations^83,84^, long-tail recommendations^85^, and cold-start scenarios^86^.

Building on their success in several domains^87–90^, LLMs are also being explored for their potential to enhance knowledge integration in various downstream tasks^91^. Their ability to understand and generate human-like text, coupled with their vast pre-trained knowledge, makes them promising candidates for incorporating external knowledge into models for tasks like physical symbolic regression^92^. This involves leveraging LLMs to not only understand the underlying physical principles and relationships expressed in text but also to guide the search for symbolic equations that accurately describe observed data^93^. By integrating knowledge from scientific literature, textbooks, or even expert explanations, LLMs could potentially help overcome limitations of purely data-driven symbolic regression methods, leading to more interpretable and physically plausible models.

Sharlin and Josephson (2025) show that while existing SR methods mostly rely on techniques such as genetic algorithms and Markov chain Monte Carlo sampling, they do not formally utilize available domain knowledge^94^. Even recent Transformer-based SR methods are limited in the context they are used for since they learn pattern matching between the dataset and the mathematical expression they have been provided with, limiting their view of similar problems that can be used. In order to tackle this challenge, the authors used the GPT-4 LLM (and GPT-4o) model to propose expressions based on the data, which were then optimized and evaluated using external Python tools, and then the results were fed back to the LLM to propose better expressions. Specifically, they added two hints to the LLM: One is the initial prompt, which is to input data for the model to suggest expressions; the second hint is the iterative hint, which is to combine the data and feedback for the model to optimize the expression. Merler et al. (2024) investigated the integration of LLMs into the SR process, utilizing contextual learning (ICL) capabilities^92^. The authors proposed the ICSR method, which is inspired by the Optimization by Prompting (OPRO) framework. That is, by providing LLMs with previously tested equations and fitness scores, they are made to generate better candidate equations and iterate until convergence or reach the computational budget. Using LLaMA3 8B as the underlying LLM, the ICSR method performed comparably well in all current popular SR benchmarks and generates expressions with low average complexity. Li et al. (2024) focused on how to simplify the complexity of finding expressions reflecting the relationships between variables from observed data^95^. The authors automatically generate expressions by describing requirements in natural language to an LLM. In particular, the authors proposed the MLLM-SR method, which uses LLMs to generate, run SR methods, and repeat until automatically-generated conditions are met.

In all these studies, the LLM is combined in a pipeline before or after the SR’s equation prediction. In this study, we aim to integrate LLM into the SR’s search process itself.

LLM-based knowledge integration to SR

In this study, we propose the integration of SR with LLMs to address to provide a physically-informed context and straightforward knowledge integration of the expert user to the SR’s search process. We hypothesis that thanks to the information LLMs contain from a wide range of scientific texts they trained on, they can evaluate the equations produced by SR in terms of dimensional consistency, scientific validity, and domain knowledge. Using this capability, LLM can provide a guiding term to the SR’s loss function. Notably, we do not claim that LLMs inherently understands physics nor that the score is physics-informed but rather that the LLMs can operate as simple-to-use mechanism by users to introduce domain knowledge to SR’s optimization process.

To this end, let us assume an optimization system with two components, an SR and an LLM, denoted S and M, respectively. Given a dataset D ⊂ R^n×m^ with n samples and m features, we wish to find an analytical equation that fits the data while also aligning with pre-defined list of physical constraints, represented as a metric functions accepting an equation as an expression tree and aware of the LLM’s implicit knowledge. To this end, the SR component receives the dataset and a loss function L that expects a dataset, an analytical function (represented as an expression tree;^96,97^) and returns a value in R^+^ which indicates how well the analytical function fits and explains the provided dataset. Formally, we propose the following loss function:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L = w_{1} e + w_{2} s + w_{3} c$$\end{document}

where w1, w2, w3 ∈ [0*,* 1] ∧ w1 + w2 + w3 = 1, e is the mean square error^97^ between the prediction of the provided analytical function and the provided dataset, s is the size of the analytical equation in terms of nodes in the expression tree, and c is the LLM-based score of how well the analytical equation fulfills a pre-defined list of physical constraints.

Importantly, in order to ensure the reliability of the LLM-based score, which is not necessarily a floating point value but rather any string, we carefully design the LLM’s prompt. To be exact, we used the following prompt:

To ensure that the prompt would work consistently across different models, we adopted a few key design principles. First, we set the temperature to zero to guarantee deterministic outputs—ensuring that models like Mistral, Llama, and Falcon would always produce the same response for the same input. We introduced a “### ROLE” section to immediately frame the assistant’s focus on scientific reasoning, helping the model interpret the task correctly from the start. In the “### METRICS” section, we clearly defined three core evaluation criteria—dimensional consistency, simplicity, and physical realism—using well-spaced numerical scales so that all models could interpret them in the same way. We included a few carefully chosen examples (a classic kinematics equation, a dimensional mismatch, and an unnecessarily nested trigonometric form) to anchor the scoring logic. Finally, we enforced a strict Python-style list format for the outputs and emphasized the “ONLY” instruction to minimize the chance of extra text generation. With these measures in place, the [dim corr, simp, sim, feedback] lists returned by Mistral, Llama, and Falcon were not only structurally consistent but also reliably comparable. Importantly, the prompt included three metrics: dim corr, simp, and sim, which stand for the dimensional correction of the equation, its simplicity, and its physical realism, respectively.

To be exact, upon receiving the prompt, each LLM returns a concise, three-element list, including a score about physical dimensionality (c1) consistency, syntactic and structural simplicity (c2), and physical realism (c3). The LLM’s overall score is computed as follows: c : = 1 − (c1 + c2 + c3)*/*3.

Notably, the physics prior enters the optimization through the auxiliary term c, which evaluates whether a candidate symbolic expression respects pre-defined physical criteria. This is directly analogous to physics- informed objectives widely used across scientific machine learning, where physical knowledge is encoded as additional penalties/constraints to bias learning toward physically consistent solutions. Specifically, the LLM is prompted to return (i) dimensional consistency (dim_c_orr)and(ii)physicalrealism/plausibility(sim)*, *whichoperationa.

Experiments

In order to evaluate the proposed method, we designed an experiment across three in silico physical scenarios (free fall, simple harmonic motion, damped waves), using three symbolic regression implementations (PySR, DEAP-GP, gplearn) and three language models (Mistral 7B, Llama 2 7B, Falcon 7B), evaluating the robustness of the proposed method over all three axes. For each combination of these components, we conducted three experiments—benchmarking the performance, comparison between different knowledge integration provided to the LLM’s prompt, and noise robustness. For all of these experiments we recorded four main metrics: mean absolute error (MAE), mean square error (MSE), coefficient of determination (R^2^), and the expression tree distance. Figure 2 presents a schematic view of the experimental setup. Due to the stochastic nature of the SR algorithm, the results are shown as the mean value of n = 25 repetitions, balancing between computational time and statistical robustness.

Fig. 2A schematic view of the experimental setup.

Experiment components

In Silico physical experiments

As a representative set of physical problems, we adopted classical Newtonian movement cases with an in- creasing level of reconstruction complexity, as reflected by the number of nodes in their expression trees. All equations show non-linear dynamics, ensuring a real-world level of complexity. Namely, the three physical cases are: free fall of a ball, simple harmonic motion, and damped wave—were. First, consider a freely falling object near Earth’s surface under constant gravity. When dropped from rest at height h [m], its speed follows

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v = \sqrt {2gh} ,$$\end{document}

where g = 9*.81 m/*s^2^. Next, the simple harmonic motion takes the following form:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x(t) = A\cos \frac{k}{m} \, t + \varphi ,$$\end{document}

with A as amplitude, k the spring constant, m the mass, and φ the phase offset^98^. Finally, a damped plane-wave propagation at fixed position x. The field obeys the following equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E(t) \, = \, E_{0} \, e^{ - \alpha t/2} \cdot \cos \,kx - wt,$$\end{document}

where α is the damping coefficient^99^. From now on, we would refer to these equations as the “ground truth” (GT) equations.

In order to generate synthetic data for each of the in silico physical scenarios, we implemented the following procedure. Initially, we conducted N = 500 experiments, uniformly sampling five independent parameters for their initial condition, including: (1) mass, randomly selected from 0*.1 to 10 kg; (2) a characteristic length (e.g. drop height, spring amplitude, or radius) between 0.01 and 0.5 m; (3) initial height/displacement, varying from 1 to 100 m; (4) a linear drag/damping coefficient from 0 to 1 kg s^−1^; and (5) time, sampled randomly between zero and the physically realistic maximum of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt {{{2h} \mathord{\left/ {\vphantom {{2h} g}} \right. \kern-0pt} g}}$$\end{document} , where g = 9.81 m s^−*2^. Each experiment is run until reaching equilibrium or 1000 steps in time. This approach provides a broad and realistic representation of conditions^100,101^.

With parameters determined, we proceeded to generate ground truth data using analytical equations corresponding to each physical scenario. These calculations provided precise, idealized responses against which the model’s performance could be evaluated. For example, free fall velocity was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v = \sqrt {2gh}$$\end{document} to validate the model’s handling of square-root dependencies. Harmonic oscillator displacement was computed using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x\, = \,A{\mathrm{cos}}(\sqrt {k/mt\, + \,\phi } )$$\end{document} to assess the model’s capture of trigonometric functions. Damped wave magnitude was calculated according to E = E0e^*−αt/*2^ cos(kx − ωt) to test the model’s ability to model exponential decay and oscillation. The results of these calculations served as our ground truth^32^. Importantly, these calculations occur until an equilibrium is achieved or the first time a state is observed twice (for period dynamics). On average, this process took around 100 steps in time, resulting in around 50,000 samples for each case. Recognizing the inevitability of measurement imperfections in practical settings, we added zero-mean Gaussian noise to our ideal signals. The noise level was set to 1% of each signal’s variability, if not stated otherwise, resulting in an average signal-to-noise ratio (SNR) of approximately 40 dB, a value representative of typical laboratory instrumentation. This process yielded noisy observations, denoted by y, which incorporate a realistic level of measurement uncertainty^102,103^. With the five independent parameters defined, we constructed a predictor matrix X. The corresponding noisy signals were assigned as the target vector y. The data were retained in their original SI units, relying on the regression models to manage any required scaling, normalization, or transformations internally^100^.

Large language models

We used three open-source and locally-depolyable LLMs—Mistral 7B, Llama 2 7B, and Falcon 7B. We picked these models due to their relative small size, which allows to run them locally on standard computer system on the one hand and their reasoning and robust training process, on the other hand. While other LLMs such as Gemini 2 and GPT-o1 that requires API (Application Programming Interface) calls may produce better results, we decided to use local LLMs as they both cheaper and ensures sensitive data is not shared with 3rd party in practical usage^104–106^. Importantly, we set the sampling temperature to zero to guarantee deterministic, repeatable outputs and to ensure that differences in plausibility scores arise solely from model knowledge rather than random sampling noise^107–109^.

Symbolic regression models

We used three popular SR models: DEAP, PySR, and gplearn. A uniform population size of 100 and 50 generations (iterations) are used for all three SR models to balance between a robust search process and computational resources. Specifically, for PySR, we used the Huber loss function and enabled a complexity penalty that softly discourages extra operators, yielding more concise formulas. For DEAP, we controlled genetic variation with a crossover probability of 0.05, a mutation probability of 0.01, a tournament selection size of 3, and a maximum tree depth of 8 to prevent bloat. Finally, for gplearn, we applied a substantial complexity penalty on operator counts so that only genuinely informative structures were preserved. Across all three SR models, we utilized the following operator set: + , − , × , ÷ , exp, log, sin, and cos. Table 1 summarizes the used hyperparameters and their values for these models.

Table 1. Hyperparameter settings for symbolic regression tools.CategoryHyperparameterDEAPPySRgplearnGeneralPopulation size100100100Iterations /Generations505050Genetic settingsCrossoverprobability0.05——Mutationprobability0.01——Tournament size3——Tree SettingsTree–height limit8——Initial tree depth1——OperatorsBinary operators + , –, *, /,protected div + , –, *, /,∧add, sub,mul, divUnary operatorsneg, log,sin, cosneg, exp, log,sin, cosneg, exp, log,sin, cos

Notably, to prevent overfitting and improve computational efficiency, we used an early-stopping rule. Optimization stopped if the best composite loss improved by less than 0.1% over three generations^110^.

Experiment types

Benchmarking

For the benchmarking of the proposed LLM-integrated SR framework, we computed all 27 (3 LLMs, 3 SR, and 3 physical scenarios) configuration and compared them to the baseline configuration of each SR model in which no LLM model is used.

LLM prompt knowledge integration experiment

Due to the fact that the exact prompt used to query the LLM response can produce highly different re- sponses, and therefore numerical values for the SR optimization process, we explored the influence of different types of knowledge in prompts on the LLM’s influence on the SR’s produced equation. To this end, we consider eight configurations. Formally, prompt A provided no additional information, acting as a baseline. Prompt B included data-column descriptions only, testing if this reminder of physical quantities improved dimensional reasoning. Prompt C added a concise, natural-language “Experiment description”, exploring whether a broader physical context improved plausibility judgments^111^. Prompt D appended the GT equation, investigating whether this reiteration of the ground truth results in the SR optimization process to converge to it^112^. In addition to the single-cue prompts, we explored combinations to investigate synergistic or competing effects. Prompt E merged B and C to assess whether pairing raw data details with contextual framing yields improved performance. Prompt F merged B and D to test whether coupling data reminders with the true formula enhances dimensional correction, potentially at the risk of over-guiding the LLM^113^. Prompt G merged C and D to evaluate if contextual description combined with ground-truth disclosure facilitates reconciling intuition with explicit answer keys^114^. Finally, Prompt H fused B, C, and D, creating a “kitchen-sink” template to determine whether information overload hinders performance or if redundancy reinforces consistent scoring^115^. The prompt variant (A–H) used for the three physical scenarios is provided as supplementary material.

Noise robustness experiment

To assess how our LLM-integrated SR model holds up when measurements get noisy, we followed the experiment design proposed by^37^. Namely, we first followed the regular data generation procedure (see section “In silico physical experiments”) but without adding noise. We then added zero-mean Gaussian noise at five levels— 1%, 2%, 3%, 4%, and 5% of each signal’s standard deviation—applying it separately to the input features and to the target variable.

Expression tree distance metric

Evaluating equations only with numerical error values (e.g., mean square error or coefficient of determination) answers merely the question of “how well they fit the data”. However, our goal is to determine how close the discovered formulas are to the data numerically and how structurally similar they are to the real underlying formula. For this purpose, we use the expression tree distance method to compare each equation we obtain with the GT equations^116,117^.

We compute the expression tree distance by first converting equations to expression trees^53^, where operations are arranged hierarchically (e.g., ” + ”, ”-”, ””). Internal nodes represent operators with left-/right subtrees, while leaf nodes contain numbers or variables. The distance is measured by simultaneously comparing the trees from the root. Identical leaf node contents yield a distance of zero; differing variables yield one. For differing numbers, the distance is α|v1 − v2| for α ∈ [0,* 1], capped at one. A leaf/internal node mismatch yields a distance of one^118^. Formally, the expression tree distance is computed as follows: mismatched operators at matched internal nodes contribute a distance of one. For identical, commutative operators (e.g., ” + ” or ””), the subtree matches with the minimum total distance (direct or cross) is selected. For non-commutative operators (e.g. ”-”, ”/”), subtree matches are strictly left-to-left and right-to- right, and subdistances are summed. These local distances are summed and normalized by the number of node comparisons, yielding a distance in [0,* 1], and aggregated from root to leaves. A zero distance signifies structural equivalence (accounting for commutativity); higher values indicate structural deviations.

Results

Table 5 in supplementary material presents an evaluation of the LLM-integrated SR framework’s performance, in terms of MAE, MSE, 1 − R^2^, and expression tree distance, across the three distinct physical experiments. Across all experiments and SR algorithms, Mistral generally achieves the best results followed by LLaMA, and then Falcon. Moreover, PySR consistently outperforms the other two SR models, often achieving the best performance. DEAP also tends to perform respectably, while gplearn often exhibits slightly weaker performance.

Table 2 presents the average performance of each group of LLM integrated SR models under eight different prompt variants in three sets of physical experiments in terms of MAE, MSE, 1 − R^2^, and expression tree distance. Notably, the richer the information provided to the LLM, the better the performance of the SR model. Of all the LMM models, Mistral achieves the best results, followed by LLaMA, then Falcon. In addition, PySR outperforms the other two SR models. Importantly, for the cases where the GT equation is provided to the LLM (Cases D, F, G, and H), an expression tree distance of 1.00 is obtained, meaning a perfect reconstruction of the GT equation. Similarly, the same results are obtained when the variable descriptions, together with the experiment description, are provided. While the expression tree distance is 1.00, the other fitness metrics are not identical (despite the fact that the predicted equation is identical), as each experiment had slightly different samples, which cause small fluctuations.

Table 2. Prompt-sensitivity: MAE, MSE, 1−R^2^ and expression-tree score for each LLM–SR pair under eight prompt variants.IdxPromptLLMSRMAEMSE1 − R^2^Expression tree distanceANo contextLLaMADEAP0.1200.0300.080.05PySR0.1000.0200.060.02gplearn0.1500.0400.100.12FalconDEAP0.1400.0400.100.10PySR0.1200.0300.080.07gplearn0.1700.0500.120.15MistralDEAP0.1100.0300.070.03PySR0.0900.0200.050.01gplearn0.1400.0400.090.11BVariabledescriptionsLLaMADEAP0.1100.0280.070.04PySR0.0900.0180.050.01gplearn0.1400.0380.090.11FalconDEAP0.1300.0380.090.08PySR0.1100.0280.070.05gplearn0.1600.0480.110.14MistralDEAP0.1000.0280.060.02PySR0.0800.0180.040.01gplearn0.1300.0380.080.10CExperimentdescriptionLLaMADEAP0.1000.0250.060.03PySR0.0800.0160.040.00gplearn0.1300.0350.080.09FalconDEAP0.1200.0350.080.06PySR0.1000.0250.060.03gplearn0.1500.0450.100.13MistralDEAP0.0900.0250.050.01PySR0.0700.0160.030.00gplearn0.1200.0350.070.08DFormula at endLLaMADEAP0.0600.0080.030.01PySR0.0400.0050.020.00gplearn0.0900.0120.050.06FalconDEAP0.0800.0100.040.02PySR0.0600.0070.030.01gplearn0.1100.0150.060.10MistralDEAP0.0500.0060.020.00PySR0.0300.0040.010.00gplearn0.0800.0100.040.04EB + CLLaMADEAP0.0400.0040.010.00PySR0.0200.0020.010.00gplearn0.0700.0070.030.00FalconDEAP0.0600.0060.020.00PySR0.0400.0040.010.00gplearn0.0900.0100.040.00MistralDEAP0.0300.0030.010.00PySR0.0100.0010.000.00gplearn0.0600.0060.020.00FB + DLLaMADEAP0.0400.0040.010.00PySR0.0200.0020.010.00gplearn0.0700.0070.030.00FalconDEAP0.0600.0060.020.00PySR0.0400.0040.010.00gplearn0.0900.0100.040.00MistralDEAP0.0300.0030.010.00PySR0.0100.0010.000.00gplearn0.0600.0060.020.00GC + DLLaMADEAP0.0400.0040.010.00PySR0.0200.0020.010.00gplearn0.0700.0070.030.00FalconDEAP0.0600.0060.020.00PySR0.0400.0040.010.00gplearn0.0900.0100.040.00MistralDEAP0.0300.0030.010.00PySR0.0100.0010.000.00gplearn0.0600.0060.020.00HB + C + DLLaMADEAP0.0400.0040.010.00PySR0.0200.0020.010.00gplearn0.0700.0070.030.00FalconDEAP0.0600.0060.020.00PySR0.0400.0040.010.00gplearn0.0900.0100.040.00MistralDEAP0.0300.0030.010.00PySR0.0100.0010.000.00gplearn0.0600.0060.020.00

Table 3 summarizes the results of the noise experiment, dividing between noise addition to the source features, the target feature, for the cases of 1% to 5% standard deviation noise addition. On average, for the DEAP SR algorithm, there is an increase of the expression tree distance from 0.07 at 1% noise to 0.18 at 5% noise. In contrast, DEAP and gplearn start at lower baselines − 0.80 and 0.90 at 1% noise level and present an even steeper increase, especially under combined noise, where their expression tree distances are as small as 0.60 to 0.40 at the 5% noise level.Table 3. Robustness evaluation of SR models under different noise types and levels (1%–5%) for the expression- tree score.ExperimentSR ModelFeature noiseTarget noiseBoth noise1%2%3%4%5%1%2%3%4%5%1%2%3%4%5%DroppingballDEAP0.170.220.260.290.300.290.370.400.420.440.460.520.550.560.57PySR0.110.150.200.220.250.150.200.250.270.300.250.300.350.370.40gplearn0.200.250.300.320.350.310.400.420.430.450.480.540.570.580.60Simple harmonicmotionDEAP0.120.170.220.240.250.220.270.300.330.350.290.340.370.390.40PySR0.070.100.120.140.150.100.140.170.200.220.200.240.270.280.30gplearn0.100.150.170.190.200.200.240.270.300.330.280.320.350.370.40ElectromagneticwaveDEAP0.100.150.190.200.220.250.300.320.340.350.350.400.420.440.45PySR0.050.080.100.110.120.120.170.200.220.250.220.270.290.300.30gplearn0.080.120.150.170.180.220.270.290.300.320.320.370.390.400.40Note: Values are expression-tree scores.

Table 4 summarizes the results of the biased (non-zero mean) sensitivity analysis, where we perturb the data using ε ∼ N (µ, σ^2^) with µ = 1% of the mean value of each feature of the data), and report performance for 1% to 5% noise levels, divided between noise added to the source features, the target variable, and both. As expected, biased noise induces a stronger degradation than the mean-zero setting: across all three experiments, the average expression-tree distance for DEAP increases from 0.16 at 1% feature noise to 0.32 at 5% feature noise, while under combined noise it increases from 0.44 (1%) to 0.60 (5%). A similar trend is observed for gplearn, which shows a steep increase under combined noise, reaching distances around 0.53–0.73 at the 5% level depending on the experiment. In contrast, PySR remains comparatively more robust under biased noise, with lower distances across all noise types and levels (e.g., 0.10 to 0.23 on average for feature noise from 1 to 5%), although it also degrades noticeably under combined noise at higher noise levels.

Discussion and conclusion

In this study, we proposed a method of integrating LLM with SR for physical cases. The method introduces an LLM-based score for the SR loss function, allowing a simple method of introducing physical knowledge to the SR’s search process. To this end, we conducted a robust evaluation of the method using three SR algorithms (DEAP, gplearn, and PySR) together with three pre-trained LLM models (DEAP, Falcon, and Mistral) on three physical dynamics (dropping ball, simple harmonic motion, and electromagnetic wave). For each combination of these three, we computed the fitting metrics in terms of MAE, MSE, and 1 − R^2^ as well as the expression tree distance between the GT and the obtained equations.

Our results suggest that augmenting SR with an LLM-based evaluation term can improve reconstruction quality in the three physical scenarios and experimental settings studied, as indicated in Table 5 in supplementary material. The improvements are not uniform: they depend on the SR engine, the chosen LLM, and the prompt variant. We therefore interpret the LLM term as a configurable prior/regularizer that can help guide search, rather than as a guarantee of consistent gains or inherently more stable optimization across SR methods. This out- come aligns with previous studies that LLM tutors improve other data-driven model’s performance^119,120^. Specifically, we find that the improvement in the SR models’ performance is consistent, indicating that a reasonable usage of the LLM in the loss function, alongside the data-fitting terms results in a relatively stable optimization task^121^. Moreover, Mistral outperforms the other LLMs, which can be attributed to its comparably large size.

Moreover, exploring how different information about the experiment influences the performance of the LLM-integrated SR method, we reproduced the expected result that more information, on average, improves the method’s performance, as indicated by Table 3. This outcome is in line with recent studies showing that prompt engineering can function as a powerful inductive bias for downstream models^122,123^. Specifically, our prompt-sensitivity study indicates that providing clearer variable definitions and experiment context often correlates with improved performance of the LLM-integrated SR variants. This pattern is consistent with the view that prompts can act as an inductive bias for the LLM scoring term^124^. However, we do not directly measure optimization stability, so we avoid mechanistic claims about “stabilization” and instead describe the effect as a bias toward equations the LLM rates as more plausible under the provided context. Furthermore, when evaluated on our synthetic benchmarks with known ground truth, improvements are sometimes reflected in greater structural similarity to the GT equation; we treat this as evidence of improved recovery in these controlled settings, not as a general claim of “accurate discovery” of physical dynamics in real-world data.

Furthermore, we obtained the results of three physical experiments using data in which noise was de- liberately introduced to the target variable, following the approach of previous studies^35^. In real-world scenarios, noise can affect the target variable, the input variables, or both. While prior work has mostly focused on target noise, the latter two cases, though more prevalent in practice, are more challenging for SR, as noise in the input variables can compound the effects of target noise. Our noise analysis revealed two key findings, as revealed by Tables 3 and 4. First, in the synthetic noise experiments (Gaussian noise at the tested levels), LLM-integrated variants often degrade less than their corresponding non-LLM baselines, suggesting a potential robustness benefit under the specific noise model and conditions studied. Second, the method’s sensitivity to noise increases as the complexity of the underlying equation grows for all SR models (in terms of the size of the expression tree of the GT equation). Importantly, while LLM increases the computation time, its relative cost remains moderate and manageable. Specifically, the LLM term in- creases computation time, and the relative overhead is configuration-dependent: it is more modest when SR iterations are computationally heavy, but can be substantial when the SR loop is fast. In our case, the LLM-to-SR runtime ratio typically remains below 0.4, indicating that the LLM does not dominate the total runtime. Moreover, all experiments were executed on a Google Colab Premium environment, suggesting that these runtimes are achievable with accessible computing resources.Table 4. Robustness evaluation of SR models under different biased noise types and levels (1%–5%) for the expression-tree score.ExperimentSR ModelFeature noiseTarget noiseBoth noise1%2%3%4%5%1%2%3%4%5%1%2%3%4%5%Dropping ballDEAP0.200.260.310.350.370.360.450.490.520.540.550.620.660.680.70PySR0.140.190.240.270.310.200.270.320.350.390.320.380.440.470.51gplearn0.240.300.350.380.420.390.480.520.540.570.570.640.680.700.73Simple harmonicmotionDEAP0.150.210.270.300.320.280.340.380.410.440.360.420.460.490.52PySR0.090.130.160.190.210.140.190.230.270.300.250.300.340.370.40gplearn0.130.180.220.240.270.250.300.350.390.430.340.400.450.490.53ElectromagneticwaveDEAP0.130.190.240.260.280.310.370.400.430.460.420.480.520.550.58PySR0.070.110.140.160.180.160.220.260.290.320.280.340.380.400.43gplearn0.110.160.200.220.250.270.330.370.400.440.380.440.480.510.54Note: Values are expression-tree distances.

Based on these results, we recommend the following practical considerations for integrating LLM into SR search of physical dynamics. First, adopt a fully informative prompt that includes clear variable definitions with their physical dimensions and a detailed experiment description. This configuration—corresponding to Prompt E in Table 2, which show promising results, as it both fitted the data well and find the GT equation exactly for all LLM and SR combinations. Second, consider the predictions of different SR as their computation is relatively cheap, and they seem to provide different results. Observing the common parts of the obtained equations from these can be useful in partial usage.

This study has several limitations that should be acknowledged. First, the quality and consistency of LLM-generated scores remain a challenge. The models were used in a zero-shot setting without fine-tuning, which likely limited their effectiveness. Although structured prompting was applied to encourage more reliable outputs, the non-deterministic nature of LLMs led to occasional invalid or syntactically incorrect suggestions, requiring prompt revisions and increasing training time. Second, the quality of LLM score was not manually or systematically evaluated by human experts. Without human verification, it is challenging to gauge whether performance gains were due to genuinely helpful scores in terms of physical realism or the score captures more technical properties, such as the equation’s length. Third, our noise sensitivity analysis assumes a Gaussian noise with mean value equal to zero, which does not capture the full complexity of possible noise profiles, providing only an approximation to the noise robustness profile of the proposed method. Finally, the experiments were conducted on synthetically generated data with known ground truths. While this allows for controlled evaluation, it does not fully reflect the complexity and noise of real-world systems. For future work, we propose extending this investigation using more powerful, cloud-hosted LLMs, which may generate higher-quality and context-aware advice and thus better demonstrate the full potential of LLM performance, especially in more complex or higher-dimensional environments^104^, and incorporating human-in-the-loop systems for expert validation. Moreover, as the LLM component alters the optimization landscape, future works can investigate and develop possible SR models that benefits from the LLM’s feedback more comparable to current SR models.

Taken jointly, this study demonstrates that LLMs can be a straightforward way to introduce textual, physics-related context into SR via an auxiliary scoring term, potentially moving from fitting data toward generating more physically plausible hypotheses. When viewed in the larger context of automatic scientific discovery, where SR and LLMs have already taken a central role, this approach further deepens their connection to one another and the task. We emphasize that the LLM score is heuristic and does not guarantee physical validity; external validation remains necessary, especially for real-world datasets. Thus, more re- fined exploration of this relationship, taking into account the SR and LLM properties as well as the nature of the physical dynamics, should further reveal the optimal design of this method.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1