Paper deep dive
MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models
Boxiao Wang, Runxiang Wang, Kai Li, Chongming Li, Zhiwei Chen, Yifan Zhang, Jian Cheng
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/3/2026, 3:09:58 AM
Summary
The paper introduces MOT-SR, a Multi-Objective Tool-augmented Symbolic Regression framework that uses Large Language Models to discover analytical equations from data. It addresses limitations of existing LLM-based SR methods by integrating external analytical tools for variable dependency analysis and employing a multi-objective evaluation module that optimizes for accuracy, complexity, and generalization via a dynamic Pareto front. The system utilizes two collaborative LLM modules: a Meta Strategy Generator for selecting tools and synthesizing strategies, and an Equation Generator for producing candidate equations. MOT-SR outperforms existing methods on 40 standard tasks and demonstrates high reliability in modeling extreme mass-ratio inspiral (EMRI) orbital dynamics.
Entities (10)
Relation Signals (7)
MOT-SR → optimizes → accuracy, complexity, and generalization
confidence 95% · ...jointly optimizing for accuracy, complexity, and generalization via a multi-objective evaluation module...
Meta Strategy Generator → partof → MOT-SR
confidence 95% · MOT-SR employs two collaborative LLM modules: a Meta Strategy Generator...
Equation Generator → partof → MOT-SR
confidence 95% · ...and an Equation Generator, which produces new candidate equations accordingly.
MOT-SR → uses → Large Language Models
confidence 95% · MOT-SR employs two collaborative LLM modules: a Meta Strategy Generator... and an Equation Generator...
MOT-SR → appliedto → EMRI orbital modeling
confidence 90% · We further validate MOT-SR on extreme mass-ratio inspiral (EMRI) orbital modeling...
MOT-SR → outperforms → existing SR methods
confidence 90% · Across 40 standard tasks, MOT-SR outperforms existing SR methods in accuracy, generalization, and efficiency.
LLM-SR → comparedwith → MOT-SR
confidence 85% · MOT-SR outperforms existing SR methods... including LLM-SR.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Symbolic Regression (SR) aims to discover analytical equations from observational data and plays a central role in scientific modeling. While recent Large Language Model (LLM) based approaches show promise, they face two limitations. First, they lack data analysis mechanisms for uncovering variable dependencies, which reduces the efficiency of equation discovery. Second, most methods rely on single-objective evaluation focused solely on fitting error. This neglect of structural complexity and generalization often causes models to converge prematurely to local optima, limiting their ability to explore the broader equation space. We propose Multi-Objective Tool-augmented Symbolic Regression (MOT-SR), a unified framework that integrates external analytical tools to extract structural priors and guide equation generation, while jointly optimizing for accuracy, complexity, and generalization via a multi-objective evaluation module that maintains a dynamic Pareto front. MOT-SR employs two collaborative LLM modules: a Meta Strategy Generator, which selects tools and synthesizes structural optimization strategies based on Pareto-optimal equations, and an Equation Generator, which produces new candidate equations accordingly. The system operates in a closed-loop manner, continuously refining both strategies and equation structures. Across 40 standard tasks, MOT-SR outperforms existing SR methods in accuracy, generalization, and efficiency. We further validate MOT-SR on extreme mass-ratio inspiral (EMRI) orbital modeling, an important problem in space-based gravitational-wave astronomy where small local errors can accumulate substantially over long-term evolution. The discovered interpretable correction achieves the lowest trajectory-level integration error on held-out configurations. These results demonstrate the potential of MOT-SR to enable reliable modeling of long-horizon scientific dynamics.
Tags
Links
- Source: https://arxiv.org/abs/2607.29561v1
- Canonical: https://arxiv.org/abs/2607.29561v1
Trouble viewing inline? Open PDF directly →
Full Text
132,170 characters extracted from source content.
Expand or collapse full text
MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models Boxiao Wang National Key Laboratory of Cognition and Decision Intelligence for Complex Systems Institute of Automation, Chinese Academy of Sciences Beijing, China wangboxiao2026@ia.ac.cn Runxiang Wang University of the Chinese Academy of Sciences School of Advanced Interdisciplinary Sciences Beijing, China wangrunxiang2023@ia.ac.cn Kai Li National Key Laboratory of Cognition and Decision Intelligence for Complex Systems Institute of Automation, Chinese Academy of Sciences Beijing, China kai.li@ia.ac.cn Chongming Li School of Astronomy and Space Science, University of the Chinese Academy of Sciences Beijing, China lichongming24@mails.ucas.ac.cn Zhiwei Chen Relativistic Astrophysics, Institute for Theoretical Physics Goethe University Frankfurt Frankfurt am Main, Germany zchen@itp.uni- frankfurt.de Yifan Zhang National Key Laboratory of Cognition and Decision Intelligence for Complex Systems Institute of Automation, Chinese Academy of Sciences Beijing, China yfzhang@nlpr.ia.ac.cn Jian Cheng National Lab of Pattern Recognition Institute of Automation, Chinese Academy of Sciences Beijing, China jcheng@nlpr.ia.ac.cn Abstract Discovering compact and interpretable equations from observa- tional data is essential for understanding complex scientific sys- tems. Symbolic regression (SR) provides a general computational framework for this task. While recent Large Language Model (LLM) based SR approaches show promise, they face two key limitations. First, they lack dedicated data analysis mechanisms for uncovering variable dependencies, which reduces the efficiency of equation discovery. Second, most methods rely on single-objective evalua- tion focused solely on fitting error. This neglect of structural com- plexity and generalization often causes models to converge prema- turely to local optima, limiting their ability to explore the broader equation space. We propose Multi-Objective Tool-augmented Sym- bolic Regression (MOT-SR 1 ), a unified framework that integrates external analytical tools to extract structural priors and guide equa- tion generation, while jointly optimizing for accuracy, complexity, and generalization via a multi-objective evaluation module that maintains a dynamic Pareto front. MOT-SR employs two collab- orative LLM modules: a Meta Strategy Generator, which selects tools and synthesizes structural optimization strategies based on Pareto-optimal equations, and an Equation Generator, which pro- duces new candidate equations accordingly. The system operates in a closed-loop manner, continuously refining both strategies and equation structures. Across 40 standard tasks, MOT-SR outperforms existing SR methods in accuracy, generalization, and efficiency. We further validate MOT-SR on extreme mass-ratio inspiral (EMRI) 1 Code is available athttps://github.com/wswbx/MOT-SR. orbital modeling, an important problem in space-based gravitational- wave astronomy where small local errors can accumulate substan- tially over long-term evolution. The discovered interpretable cor- rection achieves the lowest trajectory-level integration error on held-out configurations. These results demonstrate the potential of MOT-SR to enable reliable modeling of long-horizon scientific dynamics. Keywords Symbolic Regression, Extreme Mass-Ratio Inspirals, Gravitational Wave Modeling, Large Language Models 1 Introduction Symbolic Regression (SR) [34] aims to discover underlying math- ematical equations from data and has long been recognized as a key methodology in scientific discovery. It has been widely applied across disciplines, from identifying physical laws [ 33,45] and mod- eling chemical systems [7,12], to analyzing dynamic processes in biological or economic systems [ 51,60]. By generating compact and interpretable equations, SR enables both accurate prediction and deep insight into system behavior. SR has long been recognized as an NP-hard problem [ 59], moti- vating diverse algorithmic developments. Early approaches based on genetic programming [ 9,49] evolve equations via mutation and crossover. Reinforcement learning [44] models SR as a sequential decision-making process. Recently, Transformer-based models have enabled end-to-end learning from data to equations [ 4,20,63]. With arXiv:2607.29561v1 [cs.LG] 31 Jul 2026 Wang et al. Diversity-Aware Parent Selection LLM Generation NMSE NMSE Meta Strategy Generator Equation Generator & Evaluation Non-Dominated Sorting & Pruning Observed Data 푋,푦 (train-ID) (train-OOD) Multi-Objective Evaluation AST Length (Complexity) Search Strategy Update Frontier ( ) 풫 퓉+1 Iterative Loop (Generation 푡 → 푡 + 1 ) Data Analysis Prompting Correlation FFT/WaveletCausal Scientific Toolbox Residual Calculation LLM Reasoning Pareto Frontier Buffer ( ) 풫 퓉 푒 1 ,푒 2 ,...,푒 푁 Population of Equations Figure 1: Overview of MOT-SR. Given observed data(푋, 푦), MOT-SR maintains a Pareto-frontier buffer풫 푡 of candidate equations and multi-objective scores. At iteration푡, a Meta Strategy Generator analyzes residuals of Pareto-optimal candidates, applies scientific tools, and outputs search strat- egy. Following this search strategy, an Equation Generator & Evaluation module performs diversity-aware parent selec- tion, generates new equations, and evaluates them byNMSE on퐷 ID train and퐷 OOD train (training-derived in-domain (ID)/out-of- domain (OOD) split) and by AST length for complexity. Non- dominated sorting and pruning update the frontier to풫 푡+1 , forming a closed-loop scientific equation discovery process. the rise of large language models (LLMs), methods such as LLM- SR [52] and LaSR [18] leverage LLM’s in-context learning capabil- ities and scientific priors to perform symbolic reasoning and equa- tion generation. Despite encouraging progress, existing LLM-based SR methods face two key limitations. First, they typically lack systematic analy- sis of variable dependencies and data distributions, relying instead on problem descriptions as context. This often results in poorly constrained search spaces, which limits both the efficiency and directionality of equation exploration. Second, most approaches adopt a single-objective evaluation, typically minimizing fitting er- ror, while overlooking other critical factors such as equation com- plexity, generalization, and diversity. This can lead to overfitting and premature convergence to locally optimal solutions. To address these challenges, we propose Multi-Objective Tool- augmented Symbolic Regression (MOT-SR), a unified framework (Figure 1) inspired by the human scientific modeling process. Scien- tists typically begin by analyzing data using a variety of tools, eval- uate candidate hypotheses from multiple perspectives, and itera- tively refine their modeling direction accordingly. MOT-SR oper- ationalizes this human-like modeling paradigm through two core mechanisms and a cooperative, LLM-driven evolution system. At its core, MOT-SR first employs a tool-augmented analysis mechanism that invokes a suite of analytical tools to extract vari- able relationships from multiple complementary dimensions, such as linearity, periodicity, and causality. These insights are converted into interpretable priors and structural constraints that guide equa- tion generation. In parallel, MOT-SR adopts a multi-objective eval- uation mechanism that jointly assesses candidate equations across three dimensions: accuracy on the in-domain (ID) training subset 퐷 ID train , generalization on the out-of-domain (OOD) training subset 퐷 OOD train , and structural complexity measured via abstract syntax tree (AST) length. A Pareto front is maintained using non-dominated sorting to preserve high-quality equations that represent optimal trade-offs. To realize these mechanisms in an adaptive and iterative man- ner, MOT-SR incorporates two cooperating LLMs. Specifically, the Meta Strategy Generator analyzes residual patterns of the current Pareto-optimal candidates and leverages external scientific tools (e.g., correlation, FFT/wavelet, and causal discovery) to produces a data-driven search strategy. Guided by this strategy, the Equa- tion Generator & Evaluation module performs diversity-aware par- ent selection and generates new candidate equations, which are evaluated and integrated into the Pareto-frontier buffer via non- dominated sorting and pruning. We evaluate MOT-SR on five benchmarks spanning physics, chem- istry, biology, and materials science, using both the open-source LLaMA-3.1 [21] and the commercial GPT-4o mini [41]. Across all tasks, MOT-SR consistently outperforms traditional SR methods and recent LLM-based baselines in terms of accuracy, generaliza- tion, and equation compactness. Beyond benchmark evaluation, we apply MOT-SR to extreme mass-ratio inspiral (EMRI) orbital evolution, a challenging prob- lem in space-based gravitational-wave astronomy where small lo- cal errors accumulate over many orbital cycles. Using 58 configu- rations for equation discovery and 30 held-out configurations for a posteriori evaluation, MOT-SR identifies a compact symbolic cor- rection that achieves the lowest trajectory-level integration error— approximately three orders of magnitude lower than a neural resid- ual baseline and 26.8 times lower than LLM-SR. These results demon- strate that MOT-SR can recover transferable correction structures that remain reliable across unseen physical configurations and long horizons, highlighting the potential of interpretable equation dis- covery for demanding scientific dynamical systems. 2 Preliminaries In SR, the learning task typically starts with a dataset consisting of input-output pairs: 퐷 = (x 푖 , 푦 푖 ) 푛 푖=1 ,x 푖 ∈ ℝ 푑 , 푦 푖 ∈ ℝ, wherex 푖 denotes a푑-dimensional input vector and푦 푖 is the cor- responding scalar output. The goal is to discover an explicit ana- lytical equation푓 (⋅)such that the predicted outputŝ푦 푖 = 푓 (x 푖 ) closely approximate the ground truth푦 푖 . To assess the quality of a candidate equation, the normalized mean squared error (NMSE) is defined as NMSE(푓 , 퐷) = 1 푛 푛 ∑ 푖=1 ( 푓 (x 푖 ) − 푦 푖 휎 푦 ) 2 .(1) where휎 푦 is the standard deviation of the target values across the dataset퐷. This metric reflects the equation’s predictive accuracy, normalized by the variance of the outputs. Beyond fitting accu- racy, SR also values simplicity and generalization, seeking equa- tions that are not only accurate but also compact and transferable to unseen domains. MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models Our work builds on LLM-SR, a framework that leverages LLMs to generate symbolic equations through iterative optimization. Its core pipeline includes: •Equation skeleton generation: Structured prompts contain task-specific information (e.g., variable names, optimiza- tion goals, example equations), guiding the LLM to pro- duce physically plausible equation skeletons. •Parameter optimization: The skeletons’ free parameters are optimized (e.g., via BFGS [14] ) and scored using NMSE. •Feedback: High-quality equations are retained and reused as in-context examples, enabling iterative refinement through feedback-driven generation. While promising, LLM-based SR methods face key limitations: they lack systematic modeling of variable dependencies, leading to structurally under-informed equations; and they rely on a single- objective evaluation focused solely on fitting error, neglecting fac- tors such as complexity and generalization. These shortcomings diminish their effectiveness in solving more challenging tasks. 3 Method To address the above limitations, we proposeMOT-SR(Multi-Objective Tool-augmented Symbolic Regression), a unified framework that integrates external data analysis tools, multi-objective evaluation, and cooperative LLMs to enhance equation quality and search ef- ficiency. MOT-SR first extracts structural priors by analyzing vari- able relationships with diverse analytical tools. It then introduces a multi-objective evaluation mechanism that jointly considers fitting error, equation complexity, and generalization, dynamically main- taining a Pareto front through non-dominated sorting. Finally, two complementary LLMs work in tandem: one generates structural refinement strategies, while the other synthesizes candidate equa- tions accordingly, forming a closed-loop process that continuously guides and improves equation discovery. 3.1 Tool-Augmented Variable Analysis This component constructs structural priors by quantifying diverse variable relationships using a suite of carefully designed analytical tools. Details are provided in AppendixI. Linear correlation toolsare essential for identifying domi- nant variables and constructing interpretable regression structures. MOT-SR integrates several complementary methods to assess lin- ear dependencies from multiple statistical angles. ThePearson cor- relation coefficient[42] quantifies pairwise linear associations, es- pecially effective for Gaussian-like data.Simple linear regression andresidual variance analysis[38] evaluate predictive capacity and error stability.PCA-based explained variance[ 19] identifies the key directions of structural variance. These tools jointly establish a solid basis for linear trend detection. Nonlinear dependency toolsare employed to capture com- plex interactions essential for modeling nonlinear systems. MOT- SR integrates three complementary methods: theSpearman rank correlation[ 54] measures monotonic associations based on rank, offering robustness to noise and non-Gaussian data;mutual infor- mation[ 50] measures the overall statistical dependency between variables without assuming any parametric form; andmutual in- formation regression[ 43] quantifies the marginal contribution of each variable conditioned on others. Together, these tools guide whether nonlinear or higher-order terms should be introduced. Time-frequency analysis toolshelp detect periodicity and transient dynamics, which frequently occur in oscillatory and mul- tiscale systems. MOT-SR employs two complementary methods: Fast Fourier Transform[8] identifies dominant global frequency components, while thewavelet transform energy spectrum[11] cap- tures localized, non-stationary fluctuations. These insights support the inclusion of periodic terms (e.g., sin) in candidate equations. Causal inference toolsare used to identify whether one vari- able may influence another in a predictive or explanatory sense. MOT-SR adoptsGranger causality[16] for detecting temporal causal influence in linear time-series, andConvergent Cross Mapping[55] for identifying latent causality in nonlinear systems with potential delays. These methods provide structural signals that enhance the interpretability and explanatory power of generated equations. Dynamic complexity toolshelp assess intrinsic system rich- ness and redundancy, guiding the pruning of over-specified com- ponents in equations. MOT-SR employs three methods: theLya- punov exponent[ 61], which measures sensitivity to initial condi- tions and indicates chaotic behavior; thecorrelation dimension[17], which estimates the system’s effective degrees of freedom; andDy- namic Time Warping (DTW)[48], which evaluates time-shifted sim- ilarity between variable trajectories. Together, they offer structural cues for constructing compact and robust equations. Distribution consistency toolsevaluate whether input vari- ables behave uniformly across different input regions. MOT-SR uses theKolmogorov–Smirnov (KS) test[ 35] to detect distributional shifts between subdomains, informing the use of region-dependent structures to reflect local variations in the data distribution. Rather than introducing new tools, the key innovation of MOT- SR lies in enabling LLMs to autonomously invoke and coordinate these tools to extract dependency patterns among variables. These insights are distilled into concise guidance that informs variable selection and function composition, thereby improving responsive- ness to data characteristics and enhancing the scientific plausibil- ity of generated equations. By adaptively combining outputs from heterogeneous analyses, for instance by linking correlation mea- sures with periodicity detection, MOT-SR supports more targeted equation discovery. The impact of this tool-augmented analysis on equation generation is substantiated by the case studies in Appen- dix J.3. Looking ahead, MOT-SR opens the possibility for LLMs to synthesize new tools, further expanding the scope of SR research. 3.2 Multi-Objective Evaluation To enhance search efficiency and model quality, MOT-SR adopts a multi-objective evaluation scheme that jointly considers predic- tive accuracy, generalization, and structural simplicity. A Pareto- based selection strategy maintains a diverse set of non-dominated candidate equations, improving robustness and exploration of the solution space. We define the ID/OOD regionswithinthe training set퐷 train as a spatial split (Figure 2). The detailed configuration of the ID-OOD partition for퐷 train is provided in AppendixC. Evaluation Metrics.Unlike prior LLM-based SR methods that opti- mize only fitting error, MOT-SR employs a three-fold metric: Wang et al. OOD (n=4800) ID (n=5200) OOD (n=1977) ID (n=2023) Figure 2: Illustration of the training-derived ID/OOD split used in MOT-SR. The shaded region denotes퐷 ID train defined by per-dimension percentile bounds, and the white region denotes퐷 OOD train = 퐷 train ∖ 퐷 ID train , shown on two representative tasks. (1) In-Domain (ID) Accuracy.Measures interpolation perfor- mance within the central region of the input space. The in-domain subset퐷 ID train is selected from the middle intervals along each input dimension (Figure2). Accuracy is quantified using NMSE on퐷 ID train , denoted as NMSE ID . (2) Out-of-Domain (OOD) Generalization.Measures extrap- olation performance in the peripheral regions퐷 OOD train = 퐷 train ∖ 퐷 ID train , which correspond to the white regions in Figure 2. NMSE on this subset is denoted as NMSE OOD and reflects robustness to dis- tributional shifts, a key requirement for scientific applications. It is worth noting that NMSE OOD is calculated on the training-derived region퐷 OOD train and does not involve any held-out test data. (3) Equation Complexity.To promote interpretability, equa- tion complexity is measured via the size of its abstract syntax tree (AST) [ 40]. In MOT-SR, each equation is represented as a Python function and parsed into an AST, where each node corresponds to a variable, operator, or function call. The total node count provides a low-overhead estimate of structural complexity, guiding the search toward compact equations. Although Python AST length may dis- proportionately penalize composite function calls, we adopt it as a computationally inexpensive proxy for expression-tree complex- ity; its empirical ranking consistency and associated trade-offs are discussed in Appendix L. Pareto-Guided Multi-Objective Optimization.Based on the three met- rics, MOT-SR formulates equation discovery as a multi-objective optimization problem over a candidate set풫= 푓 1 , 푓 2 , ... , 푓 푀 , aiming to balance fitting accuracy, extrapolation, and structural simplicity: 푓 ∗ =arg min 푓 ∈풫 ( NMSE(푓 , 퐷 ID train ),NMSE(푓 , 퐷 OOD train ),ASTLen(푓 ) ) A function푓 푎 dominates푓 푏 if it performs no worse across all objectives and better on at least one. The Pareto front풫 ∗ con- sists of all non-dominated candidates in풫. To improve the qual- ity of the front, MOT-SR filters out overly simple candidates with low equation complexity that may mislead the search by applying upper-bound thresholds on NMSE(푓 , 퐷 ID train )and NMSE(푓 , 퐷 OOD train ). Remaining equations are ranked via non-dominated sorting to con- struct the current Pareto-optimal set. Details of the dynamic thresh- old construction are provided in Appendix C.4. 3.3 Cooperative LLMs for Equation Evolution To integrate variable-aware analysis with multi-objective optimiza- tion, MOT-SR employs a cooperative framework consisting of two LLMs: theMeta Strategy Generator휋 stg and theEquation Generator 휋 eq . The former extracts structural search strategies, while the lat- ter generates candidate equations accordingly, jointly driving the population toward Pareto-optimality (Algorithm1). Meta Strategy Generator.휋 stg formulates equation search strate- gies by synthesizing variable-level insights and structural abstrac- tion. At each iteration, given the current Pareto front풫 ∗ 푡−1 ,휋 stg autonomously selects a subset of tools풜⊆Toolbox, applies them to the dataset퐷 train , and summarizes key relationships as natural language descriptionsℛ varrel to guide downstream generation. In parallel,휋 stg performs structural abstraction over풫 ∗ 푡−1 , pro- ducing complementary guidanceℛ struct through: 1)Commonality Extraction, which identifies frequent substructures that define pre- vailing symbolic motifs; 2)Disparity Analysis, which diagnoses re- gional residual patterns to uncover structural weaknesses; and 3) Blind Spot Discovery, which detects unexplored symbolic compo- nents to encourage diversity. The final strategy풮 푡 = (ℛ varrel ,ℛ struct ) is passed to equation generator휋 eq to steer the next round of gen- eration. Equation Generator.Upon receiving the strategy풮 푡 from휋 stg , 휋 eq is responsible for synthesizing a new batch of candidate equa- tions풫 푡 . To promote structural diversity and prevent premature convergence, MOT-SR incorporates astructure-diversity–guided par- ent sampling module, which explicitly prioritizes structurally dis- tinct candidates during equation generation. Specifically, each equation푓 ∈풫 ∗ 푡 on the Pareto front is parsed into an abstract syntax tree (AST), from which a set of symbolic subtrees푆(푓 )is extracted. The pairwise structural dissimilarity be- tween two equations is computed via a subtree-overlap metric de- fined as: SyntaxDiv(푓 푖 , 푓 푗 ) = − |푆(푓 푖 ) ∩ 푆(푓 푗 )| |푆(푓 푖 )| . This score reflects the relative uniqueness of푓 푖 ’s structure com- pared to푓 푗 . For each equation, we compute the average structural diversity score: Score div (푓 푖 ) = 1 푁−1 ∑ 푗≠푖 SyntaxDiv(푓 푖 , 푓 푗 ). A softmax sampling procedure is then applied to construct a par- ent set풫 parent , biased toward structurally diverse equations. These parents serve as in-context examples, which, together with the strategy풮 푡 , are fed into휋 eq to generate new equation candidates: 푓 ′ ∼ 휋 eq (풮 푡 ,풫 parent ). This procedure facilitates the generation of structurally diverse and generalizable expressions, driving the con- tinuous evolution and expansion of the equation population풫. 4 Experiments 4.1 Datasets To assess MOT-SR’s performance, we adopt two sets of challenging datasets. The first includes four standard benchmarks from LLM- SR, spanning nonlinear oscillatory systems (Oscillation 1 & 2), where Oscillation 1 focuses on periodic signal composition and Os- cillation 2 introduces cross-variable interactions and non-periodic MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models Model Oscillation 1Oscillation 2E. coli growthStress-Strain Acc avg-0.001 (%)↑ NMSE↓ Acc avg-0.001 (%)↑ NMSE↓ Acc avg-0.1 (%)↑ NMSE↓ Acc avg-0.1 (%)↑ NMSE↓ GPlern0.110.09720.050.20000.761.002328.430.3496 PySR3.800.00037.020.00022.800.406870.600.0347 RAG-SR39.471.49e-60.430.02822.040.275476.280.0282 uDSR1.780.00020.360.08561.120.505959.150.0639 LaSR (Llama-3.1)2.790.74851.090.03103.440.134971.840.0320 LLM-SR (Llama-3.1)12.672.55e-58.204.70e-51.360.581576.210.0333 LLM-SR (4o-mini)11.122.07e-58.664.51e-53.240.086371.280.0491 MOT-SR (Llama-3.1)100.001.27e-1599.451.70e-106.600.020885.020.0150 MOT-SR (4o-mini)99.991.42e-1399.574.25e-106.320.017886.330.0144 Table 1: Overall performance of MOT-SR and baseline methods on four benchmarks. disturbances to increase modeling difficulty, anE. coli Growth task modeling multivariate biological dynamics with nonlinear cou- plings [37,46], and aStress-Straintask from materials science featuring piecewise nonlinear deformation behavior [1]. The second isLSR-Synth–Chemistryfrom LLM-SRBENCH [53], which consists of 36 tasks derived from a chemical kinetics base equation with progressively increasing symbolic complexity. It is specifically designed to evaluate a model’s ability to generalize across nested, unseen, and semantically rich expressions. In this work, we focus on LSR-Synth–Chemistry rather than the complete LLM-SRBENCH suite because the other three datasets are constructed variants of the four benchmark problems already included in our evaluation. By concentrating on LSR-Synth–Chemistry, which is both novel and complementary, we ensure a comprehensive yet non-redundant assessment of model performance while also tak- ing into account the practical constraints of available hardware re- sources. Full dataset descriptions are provided in Appendix Ffor completeness. 4.2 Baselines We compare MOT-SR against a range of representative baselines from both classical and LLM-based SR methods. For the four stan- dard tasks in the LLM-SR benchmark, we include GPleaRn, a clas- sical genetic programming-based SR method; PySR [18], which com- bines evolutionary search with symbolic compression; uDSR [27], which replaces DSR’s RNN policy with a pretrained Transformer and neural-guided decoding; RAG-SR [ 63], which augments equa- tion generation with structure retrieval; and LLM-SR [52]. On the more challenging LSR-Synth–Chemistry, we compare MOT-SR with leading LLM-enhanced methods, including SGA [32], which com- bines LLM-based hypothesis generation with physics-informed pa- rameter optimization via bilevel search, and LaSR [18], which ex- tracts abstract symbolic concepts from prior equations to guide hybrid LLM-evolutionary equation generation. 4.3 Evaluation Metrics We evaluate different methods using three metrics: (1) Accuracy to tolerance휏, denoted as Acc all (휏)and Acc avg (휏), and (2) Normalized Mean Squared Error (NMSE). Acc all (휏)measures task-level correct- ness by requiring all test points to satisfy the relative error bound 휏, i.e., Acc all (휏) =1 ( max 1≤푖≤푁 test | ̂푦 푖 −푦 푖 푦 푖 | ≤ 휏 ) . Acc avg (휏)computes ModelAcc all-0.1 (%)↑NMSE↓ SGA8.330.0458 LaSR27.772.77e-04 LLM-SR66.668.01e-06 MOT-SR86.113.85e-07 Table 2: Comparison on LSR-Synth–Chemistry. the proportion of test points that meet the same criterion, defined as Acc avg (휏) = 1 푁 test ∑ 푁 test 푖=1 1 (| ̂푦 푖 −푦 푖 푦 푖 | ≤ 휏 ) . We further evaluate the quality of the final Pareto front using Hypervolume (HV) [65] and Inverted Generational Distance (IGD) [64], two standard indicators in multi-objective optimization (see Ap- pendixHfor details). HV measures the volume dominated by the obtained solution set with respect to a fixed reference point, captur- ing both convergence and diversity. It reflects the overall coverage of the objective space and favors solution sets that are both well- converged and diverse. IGD computes the average distance from the ground truth equations to its nearest counterpart in the gener- ated front, emphasizing approximation accuracy with respect to the true Pareto-optimal equations. For LLM-SR, we construct a pseudo Pareto front by extracting the nondominated set at each generation and aggregating them over 100 generations, yielding a hindsight-aggregated (best-of-run) front for a stringent compari- son. 4.4 MOT-SR Configuration For fair comparison, we adopt the same LLMs across all methods: LLaMA-3.1-8B and GPT-4o-mini, covering both lightweight and high-performance scenarios. In each iteration, MOT-SR generates four candidate equations for evaluation. The total number of iter- ations is set to 2000 for standard benchmarks and 1000 for LSR- Synth–Chemistry, following the LLM-SRBench protocol. Traditional baselines are allowed more iterations to ensure convergence. Ad- ditional details, including prompt design and sampling configura- tions, are provided in Appendix MandD. 5 Findings 5.1 MOT-SR Achieves the Best Overall Performance As shown in Table1, MOT-SR consistently outperforms both clas- sical and LLM-based SR baselines, achieving significantly lower NMSE and higher accuracy across all benchmarks. For example, on Wang et al. Oscillator1 10 14 10 12 10 10 10 8 10 6 10 4 10 2 NMSE (log scale) Oscillator2 Bactgrow 10 2 10 1 10 0 NMSE (log scale) Stress-Strain PySR RAG-SR MOT-SR LLM-SR ID OOD (a) Generalization across domains under ID and OOD. 10 14 10 12 10 10 10 8 10 6 10 4 10 2 10 0 10 2 NMSE (log scale) NMSE Comparison on LSR Synth Chemistry LLM-SR MOT-SR ID OOD (b) Boxplot of generalization on 36 chemistry tasks. Figure 3: Comparison of generalization performance. Oscillator 1 with LLaMA, MOT-SR reaches an NMSE of1.27×10 −15 , far surpassing LLM-SR’s2.55 × 10 −5 . With both LLaMA and GPT- 4o-mini, MOT-SR attains over 90% accuracy on several datasets. On the more challenging LSR-Synth–ChemistRy (Table2), MOT- SR achieves the best overall performance with an average NMSE of3.85 × 10 −7 and 86.1% accuracy, significantly outperforming other LLM-based models. These results validate the effectiveness of MOT-SR’s tool-augmented multi-objective framework in enabling the discovery of more accurate equations. 5.2 MOT-SR Shows Stronger Generalization We evaluate generalization on the held-out test set under in-domain and out-of-domain conditions. In this section (Figure3), “ID/OOD” denotes the test partitions퐷 ID test and퐷 OOD test , obtained by applying percentile-based bounds from퐷 train (Sec. 3.2) to the test inputs. These partitions are used only for final reporting and are not in- volved in parameter fitting or any discovery-loop objectives. As shown in Figure 3a, MOT-SR consistently achieves the lowest NMSE across all four standard tasks under both퐷 ID test and퐷 OOD test , significantly outperforming LLM-SR and other baselines. For in- stance, on the test OOD split of Oscillator 1 (퐷 OOD test ), MOT-SR reaches an NMSE of6.20 × 10 −14 , nearly eleven orders of magnitude lower than LLM-SR (1.4 × 10 −3 ). Figure 3bfurther shows its superior me- dian NMSEs across 36 chemistry tasks in LSR-Synth, under both test ID and test OOD conditions. MOT-SR’s generalization advantage stems from two key mech- anisms. First, our multi-objective optimization explicitly incorpo- rates performance on a training-derived OOD region, guiding the 0500100015002000 Iteration 10 13 10 10 10 7 10 4 10 1 NMSE (log scale) Oscillation 1 0500100015002000 Iteration 10 11 10 9 10 7 10 5 10 3 10 1 NMSE (log scale) Oscillation 2 0500100015002000 Iteration 10 2 10 1 10 0 NMSE (log scale) E. coli Growth 0500100015002000 Iteration 10 2 10 1 NMSE (log scale) Stress-Strain MOT-SR (Llama-3.1)MOT-SR (GPT-4o-mini)LLM-SR (Llama-3.1)LLM-SR (GPT-4o-mini) Figure 4: Training convergence comparison. search toward expressions that maintain low error under distribu- tional shifts. Second, our meta-strategy module integrates variable- level analysis and structural diagnostics: the former identifies sta- ble input dependencies via analytical tools, while the latter detects effective substructures and potential failure modes across the pop- ulation. Together, these components help MOT-SR uncover sym- bolic relationships that generalize beyond distribution-specific pat- terns. 5.3 MOT-SR Improves Discovery Efficiency We compare the convergence behavior of MOT-SR and LLM-SR across four benchmark tasks. As shown in Figure4, MOT-SR not only reduces error more rapidly but also converges to lower final NMSE values. In most cases, it outperforms LLM-SR’s best results (at 2000 iterations) within the first 1000 iterations. This efficiency gain stems from MOT-SR’s meta-strategy design. By extracting key variable dependencies through analytical tools and summariz- ing structural patterns from the Pareto front, MOT-SR narrows the search space and avoids redundant exploration, steering the model toward high-quality equations with fewer iterations. These results highlight its efficiency in navigating the symbolic search space while maintaining both accuracy and generalization. We further provide the computational cost analysis in Appendix J.2. 5.4 MOT-SR Exhibits Superior Multi-objective Optimization Quality We compare MOT-SR and LLM-SR on Oscillator 1 using HV and IGD. As shown in Figure5a, MOT-SR achieves faster growth in HV compared to LLM-SR, indicating earlier and more effective discovery of diverse, high-quality non-dominated solutions. Sim- ilarly, Figure 5bshows that MOT-SR consistently reduces IGD at a faster rate, with better convergence toward the true Pareto front. LLM-SR, in contrast, plateaus earlier with higher IGD, reflecting a tendency to remain in suboptimal regions. These results demon- strate that MOT-SR significantly outperforms baselines in terms of convergence speed, solution diversity, and optimization quality. 5.5 Ablation Study To assess the contributions of key components in MOT-SR, we conduct ablation experiments on the Oscillator 1 benchmark with MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models 0500100015002000 Iterations 0.0 0.5 1.0 HV MOT-SR LLMSR (a) HV 0500100015002000 Iterations 0.0 0.5 1.0 IGD MOT-SR LLMSR (b) IGD Figure 5: Comparison of HV and IGD of different models. Full MOT-SRw/o Dataw/o Structw/o Strategy 10 14 10 12 10 10 10 8 10 6 NMSE (log scale) 1.27e-15 6.20e-14 1.55e-08 3.99e-07 1.18e-10 1.70e-07 9.76e-07 1.53e-05 ID OOD ID OOD Figure 6: Ablation results on the Oscillation 1 problem. the LLaMA backbone. Our analysis focuses on two mechanisms: the multi-objective optimization mechanism and the meta strategy generator. 5.5.1 Impact of Multi-objective Optimization.We assess the im- pact of multi-objective optimization by replacing it with a single- objective variant (w/o MultiObj) that optimizes only fitting error, excluding complexity and generalization. In contrast,w/o Multi- Objrecovers only 1 of 9 symbolic terms, compared to 4 of 6 for MOT-SR: Ground Truth ̇푣 = 0.8sin(푥) − 0.5푥 ⋅ 푣 − 0.5푣 3 − 0.2푥 3 − 푥cos(푥). MOT-SR w/o MultiObj ̇푣 = 푎 0 tanh(푎 1 푥) − 푎 2 tanh(푎 3 푣) − 푎 4 푥 − 푎 5 푣 − 푎 6 푥푣− 푎 7 (푥 − 푣) − 푎 8 (푥 + 푣) − 푎 9 푣 − (1 − 푎 9 )푥. MOT-SR ̇푣 = −푎 0 푥 + 푎 1 푣 − 푎 2 푥 3 + 푎 3 푣 3 − 푎 4 푥푣+ 푎 5 sin(푎 6 푥). This indicates multi-objective optimization is essential for discov- ering compact and interpretable equations. In contrast, single-objective optimization often leads to redundancy and overfitting. 5.5.2 Effectiveness of the Meta Strategy Generator.To further in- vestigate the contribution of the meta strategy generator, we con- duct a series of ablation studies by progressively removing its core components: the data analysis submodule (w/o Data), the structure analysis submodule (w/o Struct), and the entire module (w/o Strat- egy). As shown in Figure 6, removing the data analysis module im- pairs strategy refinement based on variable relationships, leading to a notable performance drop. Excluding the structure module re- sults in a noticeable decline in generalization and accuracy, as the model can no longer extract structural patterns from prior equa- tions to guide generation. Removing the entire meta strategy gen- erator results in prompt-only generation without feedback, pro- ducing the lowest accuracy and stability among all variants. These findings demonstrate that the data-driven and structure- guided feedback mechanisms offer complementary strengths in guiding MOT-SR’s equation generation. Evaluations of the abla- tion settings on other benchmark datasets, are provided in Appen- dixG.1, and further ablation results on the role of different tools in the toolset are reported in AppendixG.2. 5.6 Case Study: EMRI Orbital Evolution Correction Problem and scientific objective.Space-based gravitational-wave detection is a major frontier in modern astronomy, opening the millihertz band to compact systems inaccessible to ground-based detectors [47]. Extreme mass-ratio inspirals (EMRIs), in which a compact object inspirals around a much more massive black hole, are key sources in this band and probe the spacetime and physi- cal processes near the central black hole [2]. Their many orbital cycles require accurate evolution models, since small local errors can cause substantial long-term trajectory and phase deviations. We present MOT-SR as the first application of SR to EMRI orbital- evolution modeling and use this task to examine whether symbolic equation discovery can address a practical modeling problem in gravitational-wave astrophysics. Pn5AAK is widely used for its computational efficiency but may lose accuracy over long evolu- tions, while self-force calculations are more accurate but costly and currently limited in scope. We study eccentric equatorial orbits in Schwarzschild spacetime usingFast Self-forced Inspirals (FSI) [ 57] as the reference andPn5AAK[22], denoted PN5 below, as the approximation. This comparison quantifies accumulated or- bital deviations and identifies where PN5 remains reliable. We focus on correcting the evolution of the dimensionless semi- latus rectum푝. Given the eccentricity푒and mass ratio휂, the learn- ing target is the discrepancy between the local evolution rates pre- dicted by FSI and PN5, defined asΔ ̇ 푝(푝, 푒, 휂) = ̇ 푝 FSI (푝, 푒, 휂)− ̇ 푝 PN5 (푝, 푒, 휂). The discovered residual is added to the PN5 evolution rate and integrated over time to obtain the corrected trajectory푝 corr (푡). This setting evaluates whether an equation fitted to evolution-rate data remains reliable when repeatedly applied within a dynamical solver. Full details of the residual construction, unit conversion, equation search, and numerical integration are provided in Appendix E. Dataset construction.We generate Schwarzschild eccentric-EMRI trajectories with central black hole masses between10 4 and10 7 푀 ⊙ , secondary masses between5and100푀 ⊙ , initial eccentricities0 ≤ 푒 0 ≤ 0.2, and initial semi-latus recta satisfying6 + 2푒 0 + 0.05 ≤ 푝 0 ≤ 12. Each case is simulated with both FSI and PN5, and the resulting evolution rates are aligned at the same orbital states be- fore the residual labels are computed. Among the 100 generated cases, Cases 89–100 contain only one or two valid states and are removed before evaluation. The remaining 88 cases are divided at the trajectory level: 58 complete trajectories are used for equation discovery and parameter fitting, and 30 complete trajectories are held out for testing. No trajectory contributes samples to both sub- sets. The equation structure, fitted constants, and candidate selec- tion are fixed before evaluation on the 30 held-out cases. Baselines and evaluation.We compare against N, an adapted neural residual baseline motivated by the neural UDE formulation Wang et al. Evaluation groupNNMOT-SR LLM-SR (a) Overall evaluation Discovery cases (58)1.623.48e-032.65e-01 Held-out cases (30)1.201.17e-033.13e-02 All valid cases (88)1.482.69e-031.86e-01 (b) Held-out cases by trajectory length Medium (500–2,000) 4.47e-011.32e-041.09e-02 Long (2,000–10,000)1.47 8.18e-05 4.19e-03 Extra-long (> 10,000)2.811.40e-038.56e-02 Table3:Meancase-levelNMSEoftheintegrated푝(푡)trajecto- ries. Panel (a) reports results on the discovery, held-out, and complete sets. Panel (b) groups the held-out cases by trajec- tory length. introduced by Keith et al. [23]. N uses a feed-forward ResidualMLP with hidden dimensions[128, 128, 64]and SiLU activations to pre- dictΔ ̇ 푝from(푝, 푒, 휂). The predicted residual is added to the PN5 evolution rate before integration. We also report LLM-SR as an al- ternative LLM-based symbolic correction. Following the standard a posteriori evaluation protocol for sci- entific equation discovery, we assess the complete integrated tra- jectory rather than point-wise residual fitting. For each case, we ap- ply the NMSE defined in Eq. (1) to the predicted trajectory푝 corr (푡) and the FSI reference푝 FSI (푡). Results.Table3reports the mean case-level NMSE. Lower val- ues indicate closer agreement with the FSI reference trajectory. On the 30 held-out trajectories, MOT-SR achieves a mean NMSE of1.17 × 10 −3 , approximately three orders of magnitude lower than N and26.8times lower than LLM-SR. MOT-SR also con- sistently achieves the lowest mean NMSE among the three meth- ods on the 58 discovery cases and across all 88 valid cases, demon- strating its accuracy under complete-trajectory integration. More importantly, the evaluation split is constructed at the case level: each held-out trajectory represents a distinct EMRI configuration that is entirely excluded from equation discovery. The shared sym- bolic correction discovered from 58 cases remains effective when integrated over 30 unseen configurations. This result provides ev- idence that MOT-SR generalizes across EMRI cases, instead of re- lying on samples drawn from trajectories encountered during dis- covery. To examine error accumulation over longer integration hori- zons, Table3groups the held-out trajectories with more than 500 valid integration points by trajectory length. The length-stratified results reveal clear differences in long-term integration behavior. The mean NMSE of N increases monotonically from4.47 × 10 −1 on medium trajectories to1.47on long trajectories and2.81on extra-long trajectories. This trend indicates that small errors in the neural residual accumulate when the fitted correction is repeatedly applied during integration. In comparison, MOT-SR maintains a mean NMSE between8.18 × 10 −5 and1.40 × 10 −3 and achieves the lowest error in every trajectory-length group. 7 8 p case42 FSI NIT MOT-SR LLM-SR N 0.0000.0010.0020.003 t (yr) 10 5 10 1 | p p FSI |/| p FSI | (a) Case 42 7 8 9 p case53 FSI NIT MOT-SR LLM-SR N 0.0000.0050.0100.0150.020 t (yr) 10 5 10 0 | p p FSI |/| p FSI | (b) Case 53 Figure7:Aposterioriintegrationfortwoheld-outEMRIcon- figurations, comparing FSI with N, LLM-SR, and MOT-SR. Quantitative results are based on all 30 held-out cases. This behavior is consistent with the training design of MOT-SR. MOT-SR divides풟 train into ID and OOD subsets and evaluates can- didate equations for both ID accuracy and OOD generalization dur- ing discovery. This internal generalization criterion discourages equations whose accuracy is restricted to a narrow subset of the discovery configurations, yielding a correction that remains reli- able under repeated long-term integration. Compared with LLM- SR, MOT-SR reduces the mean NMSE by factors of82.6,51.2, and 61.1on the medium, long, and extra-long trajectories, respectively. The consistent advantage over LLM-SR shows that producing a symbolic expression alone does not fully realize the potential of SR for this task; the ID/OOD-aware discovery process of MOT-SR fur- ther exploits the ability of symbolic equations to capture correction structures shared across EMRI configurations. Collectively, these results show that MOT-SR preserves low integration error across different trajectory lengths and unseen EMRI configurations. Figure 7provides two illustrative examples. In both cases, the trajectory produced by MOT-SR remains closer to the FSI refer- ence over the plotted integration horizon, while the neural resid- ual leads to a visibly larger deviation. Together with the aggre- gate results, these examples provide evidence of trajectory-level transfer across unseen EMRI configurations within the considered Schwarzschild equatorial eccentric regime. 6 Conclusion We introduced MOT-SR, a unified symbolic regression framework combining tool-augmented variable analysis, multi-objective eval- uation, and cooperative LLMs for scientific equation discovery. Across standard benchmarks, MOT-SR consistently improves predictive accuracy, generalization, equation compactness, and search effi- ciency over classical and LLM-based baselines. Beyond benchmark tasks, MOT-SR discovers a compact symbolic correction for EMRI orbital evolution that substantially reduces long-horizon trajectory- level integration error and remains effective across held-out config- urations. These results demonstrate the potential of tool-augmented, multi-objective equation discovery for building interpretable and transferable models of complex scientific dynamics. Limitations and Ethical Considerations MOT-SR currently relies on a fixed set of evaluation objectives and manually specified analytical tools. Future work will intro- duce domain-specific objectives, automate tool discovery and re- finement, and extend MOT-SR to high-dimensional scientific data MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models and more complex dynamical systems. This work involves no hu- man participants or personal data. The benchmark and astrophysi- cal data are used solely for scientific equation discovery and model evaluation. References [1]B. S. Aakash, JohnPatrick Connors, and Michael D Shields. 2019. Stress-strain data for aluminum 6061-T651 from 9 lots at 6 temperatures under uniaxial and plane strain tension.Data in Brief25 (Aug 2019), 104085.doi:10.1016/j.dib.2019. 104085 [2]Pau Amaro-Seoane, Jeff Andrews, Manuel Arca Sedda, Abbas Askar, Quentin Baghi, Razvan Balasov, Imre Bartos, Simone S Bavera, Jillian Bellovary, Christo- pher PL Berry, et al. 2023. Astrophysics with the laser interferometer space antenna.Living Reviews in Relativity26, 1 (2023), 2. [3]Johan Andersson. 2000. A survey of multiobjective optimization in engineering design.Department of Mechanical Engineering, Linktjping University. Sweden (2000), 38. [4]Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. 2021. Neural Symbolic Regression that Scales. arXiv:2106.06427[cs.LG]https://arxiv.org/abs/2106.06427 [5]L. Biggio*, T. Bendinelli*, A. Neitz, A. Lucchi, and G. Parascandolo. 2021. Neural Symbolic Regression that Scales. InProceedings of 38th International Conference on Machine Learning (ICML 2021) (Proceedings of Machine Learning Research, Vol. 139). PMLR, 936–945. *equal contribution.https://proceedings.mlr.press/ v139/biggio21a.html [6]Sandra C Cerda-Flores, Arturo A Rojas-Punzo, and Fabricio Nápoles-Rivera. 2022. Applications of multi-objective optimization to industrial processes: a literature review.Processes10, 1 (2022), 133. [7]Jindou Chen, Jidong Tian, Liang Wu, ChenXinWei, Xiaokang Yang, Yaohui Jin, and Yanyan Xu. 2025. KinFormer: Generalizable Dynamical Symbolic Re- gression for Catalytic Organic Reaction Kinetics. InInternational Conference on Representation Learning, Y. Yue, A. Garg, N. Peng, F. Sha, and R. Yu (Eds.), Vol. 2025. 67058–67080.https://proceedings.iclr.c/paper_files/paper/2025/file/ a76b693f36916a5ed84d6e5b39a0dc03-Paper-Conference.pdf [8]James W. Cooley and John W. Tukey. 1965. An Algorithm for the Machine Calculation of Complex Fourier Series.Math. Comput.19 (1965), 297–301.doi:10. 1090/S0025-5718-1965-0178586-1 [9]Miles Cranmer. 2023. Interpretable Machine Learning for Science with PySR and SymbolicRegression.jl. arXiv: 2305.01582[astro-ph.IM]https://arxiv.org/ abs/2305.01582 [10]Laure Crochepierre, Lydia Boudjeloud-Assala, and Vincent Barbesant. 2022. In- teractive Reinforcement Learning for Symbolic Regression from Multi-Format Human-Preference Feedbacks. InIJCAI 2022- 31st International Joint Conference on Artificial Intelligence. Vienne, Austria. https://hal.science/hal-03695471 [11]Ingrid Daubechies. 1992.Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, USA. [12]Song Deng, Junjie Wang, Li Tao, Su Zhang, and Hongwei Sun. 2023. EV charg- ing load forecasting model mining algorithm based on hybrid intelligence.Com- puters and Electrical Engineering112 (2023), 109010. doi:10.1016/j.compeleceng. 2023.109010 [13]Mengge Du, Yuntian Chen, and Dongxiao Zhang. 2023. DISCOVER: Deep iden- tification of symbolically concise open-form PDEs via enhanced reinforcement- learning. arXiv: 2210.02181[cs.LG]https://arxiv.org/abs/2210.02181 [14]Roger Fletcher. 1987.Practical Methods of Optimization(2nd ed.). John Wiley & Sons, Chichester, New York. [15]Jenna C Fromer and Connor W Coley. 2023. Computer-aided multi-objective optimization in small molecule discovery.Patterns4, 2 (2023). [16]C. W. J. Granger. 1969. Investigating Causal Relations by Econometric Models and Cross-spectral Methods.Econometrica37, 3 (1969), 424–438. http://w. jstor.org/stable/1912791 [17]Peter Grassberger and Itamar Procaccia. 1983. Measuring the strangeness of strange attractors.Physica D: Nonlinear Phenomena9, 1 (1983), 189–208. doi:10. 1016/0167-2789(83)90298-1 [18]Arya Grayeli, Atharva Sehgal, Omar Costilla-Reyes, Miles Cranmer, and Swarat Chaudhuri. 2024. Symbolic Regression with a Learned Concept Library. arXiv:2409.09359[cs.LG]https://arxiv.org/abs/2409.09359 [19]I. T. Jolliffe. 2002.Principal Component Analysis(2nd ed.). Springer‑Verlag, New York, NY, USA.doi:10.1007/b98835 [20]Pierre-Alexandre Kamienny, Stéphane d’Ascoli, Guillaume Lample, and François Charton. 2022. End-to-end symbolic regression with transformers. arXiv:2204.10532[cs.LG]https://arxiv.org/abs/2204.10532 [21]Paul Kassianik, Baturay Saglam, Alexander Chen, Blaine Nelson, Anu Vel- lore, Massimo Aufiero, Fraser Burch, Dhruv Kedia, Avi Zohary, Sajana Weerawardhena, Aman Priyanshu, Adam Swanda, Amy Chang, Hyrum Anderson, Kojin Oshiba, Omar Santos, Yaron Singer, and Amin Kar- basi. 2025. Llama-3.1-FoundationAI-SecurityLLM-Base-8B Technical Report. arXiv: 2504.21039[cs.CR]https://arxiv.org/abs/2504.21039 [22]Michael L Katz, Alvin JK Chua, Lorenzo Speri, Niels Warburton, and Scott A Hughes. 2021. Fast extreme-mass-ratio-inspiral waveforms: New tools for milli- hertz gravitational-wave data analysis.Physical Review D104, 6 (2021), 064047. [23]Brendan Keith, Akshay Khadse, and Scott E. Field. 2021. Learning orbital dy- namics of binary black hole systems from gravitational wave measurements. Wang et al. Phys. Rev. Research3 (Nov 2021), 043101. Issue 4.doi:10.1103/PhysRevResearch. 3.043101 [24]Michael Kommenda, Andreas Beham, Michael Affenzeller, and Gabriel Kron- berger. 2015. Complexity measures for multi-objective symbolic regression. In International Conference on Computer Aided Systems Theory. Springer, 409–416. [25]J.R. Koza. 1990. Genetically breeding populations of computer programs to solve problems in artificial intelligence. In[1990] Proceedings of the 2nd International IEEE Conference on Tools for Artificial Intelligence. 819–827.doi:10.1109/TAI.1990. 130444 [26]Jiří Kubalík, Erik Derner, and Robert Babuška. 2021. Multi-objective symbolic regression for physics-aware dynamic modeling.Expert Systems with Applica- tions182 (2021), 115210. [27]Mikel Landajuela, Chak Shing Lee, Jiachen Yang, Ruben Glatt, Claudio P San- tiago, Ignacio Aravena, Terrell Mundhenk, Garrett Mulcahy, and Brenden K Petersen. 2022. A Unified Framework for Deep Symbolic Regression. In Advances in Neural Information Processing Systems, S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Eds.), Vol. 35. Curran Associates, Inc., 33985–33998. https://proceedings.neurips.c/paper_files/paper/2022/file/ dbca58f35bddc6e4003b2d80e42f838-Paper-Conference.pdf [28]Mikel Landajuela, Brenden K. Petersen, Soo K. Kim, Claudio P. Santiago, Ruben Glatt, T. Nathan Mundhenk, Jacob F. Pettit, and Daniel M. Faissol. 2021. Improv- ing exploration in policy gradient search: Application to symbolic optimization. arXiv:2107.09158[cs.LG]https://arxiv.org/abs/2107.09158 [29]Wenqiang Li, Weijun Li, Linjun Sun, Min Wu, Lina Yu, Jingyi Liu, Yanjie Li, and Song Tian. 2023. Transformer-based model for symbolic regression via joint supervised learning. InInternational Conference on Learning Representa- tions.https://api.semanticscholar.org/CorpusID:259298765 [30]Xi Lin, Zhiyuan Yang, and Qingfu Zhang. 2022. Pareto Set Learning for Neural Multi-objective Combinatorial Optimization. arXiv:2203.15386[cs.LG]https: //arxiv.org/abs/2203.15386 [31]Thibaut Lust and Jacques Teghem. 2012.The multiobjective multi- dimensional knapsack problem: a survey and a new approach.In- ternational Transactions in Operational Research19, 4 (2012), 495–520. arXiv: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1475-3995.2011.00840.x doi:10.1111/j.1475-3995.2011.00840.x [32]Pingchuan Ma, Tsun-Hsuan Wang, Minghao Guo, Zhiqing Sun, Joshua B. Tenen- baum, Daniela Rus, Chuang Gan, and Wojciech Matusik. 2024. LLM and Simula- tion as Bilevel Optimizers: A New Paradigm to Advance Physical Scientific Dis- covery. InProceedings of the 41st International Conference on Machine Learning (Proceedings of Machine Learning Research, Vol. 235), Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (Eds.). PMLR, 33940–33962. https://proceedings.mlr.press/ v235/ma24m.html [33]Nour Makke and Sanjay Chawla. 2024. Data-driven discovery of Tsallis-like distribution using symbolic regression in high-energy physics.PNAS Nexus 3, 11 (10 2024), pgae467. arXiv: https://academic.oup.com/pnasnexus/article- pdf/3/11/pgae467/60816181/pgae467.pdfdoi:10.1093/pnasnexus/pgae467 [34]Nour Makke and Sanjay Chawla. 2024. Interpretable scientific discovery with symbolic regression: a review.Artificial Intelligence Review57 (01 2024). doi:10. 1007/s10462-023-10622-0 [35]Frank J. Massey. 1951. The Kolmogorov-Smirnov Test for Goodness of Fit.J. Amer. Statist. Assoc.46, 253 (1951), 68–78.http://w.jstor.org/stable/2280095 [36]Matteo Merler, Katsiaryna Haitsiukevich, Nicola Dainese, and Pekka Marttinen. 2024. In-Context Symbolic Regression: Leveraging Large Language Models for Function Discovery. InProceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 4: Student Research Workshop). Associa- tion for Computational Linguistics, 589–606. doi:10.18653/v1/2024.acl-srw.49 [37]Jacques Monod. 1949. THE GROWTH OF BACTERIAL CULTURES.Annual Review of Microbiology3, Volume 3, 1949 (1949), 371–394.doi:10.1146/annurev. mi.03.100149.002103 [38]Douglas C. Montgomery, Elizabeth A. Peck, and G. Geoffrey Vining. 2013.In- troduction to Linear Regression Analysis(5 ed.). Wiley, Hoboken, NJ. [39]T. Nathan Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio P. San- tiago, Daniel M. Faissol, and Brenden K. Petersen. 2021.Symbolic Regression via Neural-Guided Genetic Programming Population Seeding. arXiv: 2111.00053[cs.NE]https://arxiv.org/abs/2111.00053 [40]Iulian Neamtiu, Jeffrey S. Foster, and Michael Hicks. 2005. Understanding source code evolution using abstract syntax tree matching. InProceedings of the 2005 International Workshop on Mining Software Repositories(St. Louis, Mis- souri)(MSR ’05). Association for Computing Machinery, New York, NY, USA, 1–5. doi:10.1145/1083142.1083143 [41]OpenAI. 2024. GPT‑4o mini: advancing cost‑efficient intelligence.https: //openai.com/index/gpt-4o-mini-advancing-cost-efficient-intelligence/ . Ac- cessed: 2025‑07‑29. [42]Karl Pearson. 1895. Note on Regression and Inheritance in the Case of Two Parents.Proceedings of the Royal Society of London58 (1895), 240–242.http: //w.jstor.org/stable/115794 [43]Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay. 2011. Scikit‑learn: Machine Learning in Python.Journal of Machine Learning Research12 (2011), 2825–2830. [44]Brenden K. Petersen, Mikel Landajuela, T. Nathan Mundhenk, Claudio P. San- tiago, Soo K. Kim, and Joanne T. Kim. 2021. Deep symbolic regression: Re- covering mathematical expressions from data via risk-seeking policy gradients. arXiv:1912.04871[cs.LG]https://arxiv.org/abs/1912.04871 [45]Julia Reuter, Hani Elmestikawy, Fabien Evrard, Sanaz Mostaghim, and Berend van Wachem. 2023. Graph Networks as Inductive Bias for Genetic Program- ming: Symbolic Models for Particle-Laden Flows. InGenetic Programming, Gisele Pappa, Mario Giacobini, and Zdenek Vasicek (Eds.). Springer Nature Switzerland, Cham, 36–51. [46]L Rosso, J. R. Lobry, S Bajard, and J. P. Flandrois. 1995. Convenient Model To Describe the Combined Effects of Temperature and pH on Microbial Growth. Applied and Environmental Microbiology61, 2 (Feb 1995), 610–6.doi:10.1128/ aem.61.2.610-616.1995 [47]Wen-Hong Ruan, Zong-Kuan Guo, Rong-Gen Cai, and Yuan-Zhong Zhang. 2020. Taiji program: Gravitational-wave sources.International Journal of Modern Physics A35, 17 (2020), 2050075. [48]H. Sakoe and S. Chiba. 1978. Dynamic programming algorithm optimization for spoken word recognition.IEEE Transactions on Acoustics, Speech, and Signal Processing26, 1 (1978), 43–49.doi:10.1109/TASSP.1978.1163055 [49]Michael Schmidt and Hod Lipson. 2009.Distilling Free-Form Nat- ural Laws from Experimental Data.Science324, 5923 (2009), 81– 85. arXiv:https://w.science.org/doi/pdf/10.1126/science.1165893doi:10.1126/ science.1165893 [50]C. E. Shannon. 1948. A mathematical theory of communication.The Bell System Technical Journal27, 3 (1948), 379–423. doi:10.1002/j.1538-7305.1948.tb01338.x [51]Hao Shi, Weili Song, Xinting Zhang, Jiahe Shi, Cuicui Luo, Xiang Ao, Hamid Arian, and Luis Seco. 2024. AlphaForge: A Framework to Mine and Dynamically Combine Formulaic Alpha Factors. arXiv: 2406.18394[q-fin.CP]https://arxiv. org/abs/2406.18394 [52]Parshin Shojaee, Kazem Meidani, Shashank Gupta, Amir Barati Farimani, and Chandan K Reddy. 2025. LLM-SR: Scientific Equation Discovery via Program- ming with Large Language Models. arXiv: 2404.18400[cs.LG]https://arxiv.org/ abs/2404.18400 [53]Parshin Shojaee, Ngoc-Hieu Nguyen, Kazem Meidani, Amir Barati Fari- mani, Khoa D Doan, and Chandan K Reddy. 2025. LLM-SRBench: A New Benchmark for Scientific Equation Discovery with Large Language Models. arXiv: 2504.10415[cs.CL]https://arxiv.org/abs/2504.10415 [54]C. Spearman. 1904. The Proof and Measurement of Association between Two Things.The American Journal of Psychology15, 1 (1904), 72–101.http://w. jstor.org/stable/1412159 [55]George Sugihara, Robert May, Hao Ye, Chih hao Hsieh, Ethan Deyle, Michael Fogarty, and Stephan Munch. 2012.Detecting Causality in Complex Ecosystems.Science338, 6106 (2012), 496– 500.arXiv: https://w.science.org/doi/pdf/10.1126/science.1227079 doi:10.1126/science.1227079 [56]Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. 2021. SymbolicGPT: A Generative Transformer Model for Symbolic Regression. arXiv: 2106.14131[cs.LG]https://arxiv.org/abs/2106.14131 [57]Maarten van de Meent and Niels Warburton. 2018. Fast self-forced inspirals. Classical and Quantum Gravity35, 14, Article 144003 (July 2018), 144003 pages. arXiv: 1802.05281[gr-qc]doi:10.1088/1361-6382/aac8ce [58]Martin Vastl, Jonáš Kulhánek, Jiří Kubalík, Erik Derner, and Robert Babuška. 2024. SymFormer: End-to-End Symbolic Regression Using Transformer-Based Architecture.IEEE Access12 (2024), 37840–37849. doi:10.1109/ACCESS.2024. 3374649 [59]Marco Virgolin and Solon P. Pissis. 2022. Symbolic Regression is NP-hard. arXiv:2207.01018[cs.NE]https://arxiv.org/abs/2207.01018 [60]Ylva Wahlquist, Jesper Sundell, and Kristian Soltesz. 2024. Learning pharmaco- metric covariate model structures with symbolic regression networks.Journal of Pharmacokinetics and Pharmacodynamics51, 2 (2024), 155–167.doi:10.1007/ s10928-023-09887-3 [61]Alan Wolf, Jack B. Swift, Harry L. Swinney, and John A. Vastano. 1985. Determin- ing Lyapunov exponents from a time series.Physica D: Nonlinear Phenomena 16, 3 (1985), 285–317. doi:10.1016/0167-2789(85)90011-9 [62]Shunyu Yao, Fei Liu, Xi Lin, Zhichao Lu, Zhenkun Wang, and Qingfu Zhang. 2025. Multi-objective Evolution of Heuristic Using Large Language Model. arXiv: 2409.16867[cs.AI]https://arxiv.org/abs/2409.16867 [63]Hengzhe Zhang, Qi Chen, Bing XUE, Wolfgang Banzhaf, and Mengjie Zhang. 2025. RAG-SR: Retrieval-Augmented Generation for Neural Symbolic Regres- sion. InThe Thirteenth International Conference on Learning Representations. https://openreview.net/forum?id=NdHka08uWn MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models [64]Eckart Zitzler, Kalyanmoy Deb, and Lothar Thiele. 2000. Comparison of Multiob- jective Evolutionary Algorithms: Empirical Results.Evolutionary Computation 8, 2 (2000), 173–195.doi:10.1162/106365600568202 [65]E. Zitzler and L. Thiele. 1999. Multiobjective evolutionary algorithms: a com- parative case study and the strength Pareto approach.IEEE Transactions on Evolutionary Computation3, 4 (1999), 257–271.doi:10.1109/4235.797969 Wang et al. Appendix A Related Work A.1 Symbolic Regression Traditional SR methods mainly rely on evolutionary algorithms, reinforcement learning [44], and Transformers [5]. For instance, genetic programming [25] formulates equation discovery as an evolutionary search over tree-based representations, refining structures via mutation and crossover. Reinforcement learning-based symbolic regression was first introduced by Petersen et al. [44] and has since developed into various policy-optimization frameworks [10,13,28,39]. More recently, Transformer-based models [20,29,56,58,63] have been applied to SR, leveraging large-scale pretraining to improve equation generation. However, these models typically lack mechanisms to incorporate physical priors. With advances in natural language processing, LLM-based SR methods such as LLM-SR [52], LaSR [18], and ICSR [36] have emerged. LLM-SR utilizes scientific priors encoded in LLMs to generate plausible equation forms, followed by data-driven parameter fitting. LaSR introduces abstract concept generation to guide hypothesis construction, while ICSR formats training data as in-context prompts to induce function generation. Yet these approaches still rely heavily on pre-trained knowledge and lack explicit modeling of variable relationships or structured reasoning over data, which limits both their search efficiency and generalization ability. A.2 Multi-objective Optimization Multi-objective optimization is a fundamental research direction in the optimization community and has been widely applied in areas such as combinatorial optimization [30,31,62] and industrial design [3,6,15]. In the context of SR, Kommenda et al. systematically compared various complexity measures for multi-objective SR and proposed a new metric that preserves semantic information while improving search efficiency [24]. Kubalík et al. introduced a physics-aware multi-objective approach [26] that jointly optimizes model accuracy and physical consistency, enhancing interpretability and reliability. MOT-SR is the first framework to integrate LLMs into multi-objective SR, enabling data-driven equation discovery with improved accuracy, simplicity, and generalization. B Algorithmic Pseudocode for MOT-SR Algorithm 1MOT-SR 1:Input:Training set퐷; Maximum iterations푇;푛samples per iteration; Initial population풫 0 (optional); Meta Strategy Generator휋 stg ; Equation Generator휋 eq ;Toolbox: Data analysis toolset 2:Output:Approximate Pareto front풫 ∗ 3:Initialize풫 0 4:for푡 = 1, ... , 푇do 5:푅푒푠 ←EvaluateEquationResiduals(풫 ∗ 푡−1 , 퐷) 6:ℛ varrel ← 휋 stg .DataAnalysis(Toolbox, 푅푒푠,풫 ∗ 푡−1 ) 7:ℛ struct ← 휋 stg .GenerateStructuralPrompts(풫 ∗ 푡−1 ) 8:풮 푡 ← (ℛ varrel ,ℛ struct ) 9:풫 parent ←ParentSelection(풫 ∗ 푡−1 ) 10:풫 푡 ←풫 ∗ 푡−1 11:for푖 = 1, ... , 푛do 12:푓 ∼ 휋 eq (풮 푡 ,풫 parent ) 13:푠푐표푟푒 ←MultiObjectiveEvaluation(푓 , 퐷) 14:풫 푡 ←풫 푡 ∪ 푓 , 푠푐표푟푒 15:end for 16:풫 ∗ 푡 ← PopulationManagement ( 풫 푡 ) 17:end for 18:풫 ∗ ←풫 푇 C Data Partition and Thresholding Protocol for ID/OOD Evaluation C.1 Motivation and Protocol Definition In scientific discovery, a valid symbolic model must not only fit the observed data (interpolation) but also capture the underlying physical laws to predict system behaviors in unobserved regions (extrapolation). To rigorously evaluate this capability, we adopt apercentile-based masking strategyto partition the global data distribution풟 total into In-Domain (ID) and Out-of-Domain (OOD) subsets. Unlike random splitting, which results in independent and identically distributed (i.i.d.) subsets, our protocol enforces a strictly defined distributional shiftbetween training (ID) and testing (OOD), thereby assessing the model’s robustness against extrapolation risks. MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models C.2 Partition Methodology For each task, letx∈ ℝ 푑 denote the input variables and푦the target. The protocol is defined based onper-dimensionmarginal percentiles of the training inputs in퐷 train . Percentile candidates.We iterate over percentile values 푝 ∈ 10, 11, ... , 49.(2) Per-dimension central intervals.For each input dimension푗 ∈ 1, ... , 푑, we compute lower and upper bounds from the training inputs: ℓ 푗 (푝) =percentile(푋 푗 , 푝),푢 푗 (푝) =percentile(푋 푗 , 100 − 푝),(3) where푋 푗 denotes the set of training values of the푗-th coordinate. Training-derived ID region as a Cartesian product.We define the training-derived ID subset as the Cartesian product of these per- dimension central intervals: 퐷 ID train (푝) = (x, 푦) ∈ 퐷 train ∶ ℓ 푗 (푝) ≤ 푥 푗 ≤ 푢 푗 (푝), ∀푗 .(4) Geometrically,퐷 ID train (푝)forms a central hyper-rectangle in the input space (Figure 2). Training-derived OOD region as the complement within퐷 train .The training-derived OOD subset is defined as the spatial complement of the ID regionwithin the training set: 퐷 OOD train (푝) = 퐷 train ∖ 퐷 ID train (푝).(5) This subset corresponds to peripheral regions along at least one dimension, inducing an explicit distributional shift relative to퐷 ID train (푝). Selecting the split percentile.Among candidate percentiles, we select푝 ⋆ to balance the sample counts in the two regions: 푝 ⋆ =arg max 푝∈10,...,49 min ( | 퐷 ID train (푝) | , | 퐷 OOD train (푝) | ) ,(6) which encourages a non-trivial OOD subset while retaining sufficient ID samples for stable fitting. Finally, we set 퐷 ID train ∶= 퐷 ID train (푝 ⋆ ),퐷 OOD train ∶= 퐷 OOD train (푝 ⋆ ).(7) C.3 Usage in MOT-SR and No-Leakage Guarantee MOT-SR uses the above training-derived split as follows: •Parameter fitting (discovery loop).Equation constants are fittedonlyon퐷 ID train . •Generalization feedback (discovery loop).Extrapolation feedback is computed on퐷 OOD train to guide multi-objective selection, but this subset isneverused for parameter fitting. •Final evaluation (locked test).Benchmark test sets (e.g.,퐷 ID test and퐷 OOD test ) are held out and usedonlyfor final reporting. Therefore, the “OOD” objective used inside the discovery loop refers strictly to the training-derived region퐷 OOD train (Eq.5), and does not involve any held-out test data. C.4 Dynamic NMSE Thresholding for Pre-Filtering To avoid flooding the Pareto selection with clearly underfitting candidates, we apply a looseupper-boundNMSE threshold in each generation. For every candidate푖, we compute its worst-case training error across the two training-derived regions: 푒 푖 =max ( NMSE (푖) ID ,NMSE (푖) OOD ) ,(8) and define the current best worst-case error in the population as 푒 ⋆ =min 푖 푒 푖 .(9) We then construct a dynamic upper bound 푇 = 10 0.5log 10 (푒 ⋆ ) = √ 푒 ⋆ ,(10) so that푇is on the order of √ 푒 ⋆ in log-scale. Any candidate satisfying max ( NMSE (푖) ID ,NMSE (푖) OOD ) > 푇(11) is treated as clearly underfitting; in implementation we assign it a very large objective value so it cannot enter the Pareto front. This step does not penalize low-NMSE solutions, and only removes candidates whose errors are orders of magnitude worse than the current best. Wang et al. D MOT-SR Configuration and Language Model Details D.1 MOT-SR Configuration We implement MOT-SR with both open-source and commercial LLM backbones:LLaMA-3.1-8B-InstructandGPT-4o-mini. TheLLaMA- 3.1model is quantized and deployed locally on NVIDIA H100 80GB GPUs for efficient inference. AndGPT-4o-miniis accessed via the OpenAI API, providing high-quality reasoning without requiring local resources. Following the LLM-SRBench protocol, the maximum number of optimization iterations is set to 2000 for the four main benchmark datasets, and 1000 for the LSR-Synth–Chemistry dataset. In each iteration, the Meta Strategy Generator selects 3 tools from the tool set to analyze variable-level relationships and generates a natural language strategy prompt based on both data characteristics and symbolic structural patterns. Both theMeta Strategy Generatorand Equation Generatorin MOT-SR utilize LLM decoding with top-푘 = 30, top-푝 = 0.3, and temperature 0.6. In each iteration, the Equation Generator produces 4 candidate expressions. During Pareto frontier construction, we apply an NMSE-based filtering mechanism to eliminate trivial expressions: for each candidate, we compute the base-10 logarithm of its minimal NMSE (across ID and OOD), and discard it if it exceeds10 ⌈0.5⋅min log 10 NMSE⌉ , where the minimum is taken over the current candidate pool. D.2 Large Language Models MOT-SR employs two backbone LLMs:LLaMA-3.1-8B-InstructandGPT-4o-mini, covering both open-source and API-accessible com- mercial models. The LLaMA-3.1 model is locally quantized to 4-bit precision and deployed on NVIDIA H100 80GB GPUs, requiring ap- proximately 8GB of VRAM during inference, thus allowing execution on consumer-grade hardware with appropriate quantization. The GPT-4o-mini model is accessed via OpenAI’s API, offering high-quality reasoning with minimal deployment overhead. E Symbolic Equation Discovery for EMRI Orbital Evolution E.1 EMRI Orbital-Evolution Dataset Here, we provide the preprocessing details used to construct the symbolic regression dataset. Each case is simulated using both FSI and PN5. The resulting evolution rates are aligned at the same orbital states(푝, 푒, 휂)before the residual labels are computed. If the output points of the two models do not coincide, they are aligned through interpolation or recomputation. The original data retain time and evolution rates in seconds. To express the rates with respect to the dimensionless time푡/푀, we define 푀 sec = 퐺푀 푐 3 .(12) For푋 ∈ FSI,PN5, the dimensionless semi-latus-rectum evolution rate is ̇ 푝 푋 = 푀 sec ̇ 푝 sec 푋 ,(13) where ̇ 푝 sec 푋 = 푑푝/푑푡 sec . The symbolic regression target is then constructed as Δ ̇ 푝 = 푀 sec ( ̇ 푝 sec FSI − ̇ 푝 sec PN5 ) .(14) The resulting regression samples have the form (푝, 푒, 휂) ⟼ Δ ̇ 푝,휂 = 휇 푀 .(15) Among the 100 generated cases, Cases 89–100 contain only one or two valid states and are removed before evaluation. The remaining 88 cases are split at the trajectory level: 58 complete trajectories are used for equation discovery and parameter fitting, while 30 complete trajectories are held out for testing. No trajectory contributes samples to both subsets. The equation structure, fitted constants, and candidate selection are fixed before evaluation on the held-out cases. Only푝,푒, and휂are used as inputs to the symbolic regression model. The masses푀and휇determine the mass ratio휂, while푀additionally defines the time-conversion factor푀 sec . Extrinsic parameters, including redshift, sky location, source orientation, distance, and initial phases, are excluded because they do not determine the orbit-averaged evolution rate considered here. E.2 Formulating EMRI Orbital Evolution as Symbolic Regression Predicting a local correction instead of a trajectory.Directly regressing the full trajectory푝(푡)would entangle the local modeling error with numerical integration and the selected initial conditions. We instead search for a local and configuration-independent symbolic residual 푓 휃 ∶ (푝, 푒, 휂) ⟼ Δ ̇ 푝,(16) where푓 휃 is an explicit symbolic expression and휃denotes its fitted constants. The symbolic hypothesis is trained to approximate the rate discrepancy in Equation ( 14): 푓 휃 (푝, 푒, 휂) ≃ ̇ 푝 FSI (푝, 푒, 휂) − ̇ 푝 PN5 (푝, 푒, 휂).(17) This formulation asks the symbolic model to identify the systematic part of the approximation error shared across EMRI configurations, rather than memorizing a separate time-dependent correction for each trajectory. MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models Corrected orbital-evolution rate.Given a candidate symbolic expression, the corrected PN5 evolution rate is reconstructed as ̇ 푝 corr (푝, 푒, 휂) = ̇ 푝 PN5 (푝, 푒, 휂) + 푓 휃 (푝, 푒, 휂).(18) When the orbital evolution is integrated in physical time, the corresponding rate is ̇ 푝 sec corr = ̇ 푝 PN5 + 푓 휃 (푝, 푒, 휂) 푀 sec .(19) The symbolic correction is therefore repeatedly evaluated at the current orbital state during numerical integration, rather than being applied once as a post-processing adjustment. Symbolic fitting and multi-objective evaluation.For a set of residual samples풟, the point-wise prediction error of a candidate expression is measured by NMSE Δ ̇ 푝 ( 푓 휃 ,풟 ) = ∑ 푖 [ 푓 휃 (푝 푖 , 푒 푖 , 휂 푖 ) − Δ ̇ 푝 푖 ] 2 ∑ 푖 [ Δ ̇ 푝 푖 −Δ ̇ 푝 풟 ] 2 .(20) Within the 58-trajectory discovery set, MOT-SR applies the same training-derived ID/OOD partition used in the general framework. The free constants are optimized using the ID subset, while the OOD subset supplies extrapolation feedback. Candidate expressions are then jointly evaluated according to ( NMSE ID ,NMSE OOD ,ASTLen ) ,(21) so that the search favors expressions that balance local fitting accuracy, cross-region generalization, and structural compactness. A posteriori trajectory evaluation.A low point-wise residual error does not necessarily guarantee an accurate long-term orbit, because small systematic errors can accumulate when the learned expression is repeatedly evaluated inside the dynamical solver. We therefore perform ana posteriorievaluation by inserting Equation (18) into the orbital-evolution procedure and integrating from the initial conditions of each complete trajectory. This produces the corrected trajectory푝 corr (푡). For each EMRI case, trajectory-level accuracy is measured against the FSI reference using NMSE 푝(푡) = ∑ 푘 [ 푝 corr (푡 푘 ) − 푝 FSI (푡 푘 ) ] 2 ∑ 푘 [ 푝 FSI (푡 푘 ) − 푝 FSI ] 2 .(22) The equation structure, fitted constants, and final candidate selection are fixed before integration on the 30 held-out configurations. Hence, Equation (22) evaluates whether a single symbolic correction discovered from the 58 discovery cases remains reliable under repeated integration across unseen EMRI systems and different trajectory lengths. F Datasets F.1 Nonlinear Oscillator Oscillatory systems with nonlinear damping are foundational in physics and engineering for modeling the motion of objects subjected to restoring and dissipative forces. These systems are governed by second-order differential equations of the form̈푥 + 푓 (푡, 푥, ̇푥) = 0, where the nonlinear function푓encapsulates dynamic interactions among position, velocity, and possibly time. To evaluate a model’s ability to recover such complex dynamics, two synthetic oscillator tasks are used. The first system is defined by: ̇푣 = 퐹sin(휔푥) − 훼푣 3 − 훽푥 3 − 훾푥푣 − 푥cos(푥) (퐹 = 0.8, 훼 = 0.5, 훽 = 0.2, 훾 = 0.5, 휔 = 1.0) The second system follows: ̇푣 = 퐹sin(휔푡) − 훼푣 3 − 훽푥푣 − 훿푥exp(훾푥) (퐹 = 0.3, 훼 = 0.5, 훽 = 1.0, 훿 = 5.0, 훾 = 0.5, 휔 = 1.0) with initial conditions푥 = 0.5,푣 = ̇푥 = 0.5, and simulation time푡 ∈ [0, 50]. These equations exhibit rich nonlinear structures and variable couplings, making them ideal benchmarks for testing symbolic reasoning and generalization beyond simple oscillatory behavior. F.2 Bacterial Growth Accurately capturing the dynamics of E. coli proliferation under varying environmental conditions is of critical importance in areas such as biotechnology and microbiological risk assessment. This benchmark reflects realistic yet challenging biological modeling, where growth depends on multiple interacting factors. The dataset models population change via: 푑퐵 푑푡 = 푓 퐵 (퐵) ⋅ 푓 푆 (푆) ⋅ 푓 푇 (푇) ⋅ 푓 pH (pH) where퐵is bacterial density,푆is nutrient concentration,푇is temperature, and pH denotes acidity. Each term accounts for a distinct physi- ological influence on growth. Wang et al. The explicit expression used is: 푑퐵 푑푡 = 휇 max 퐵 ( 푆 퐾 푆 + 푆 ) ( tanh푘(푇 − 푥 0 ) 1 + 푐(푇 − 푥 decay ) 4 ) exp ( −|pH−pH opt | ) sin ( (pH−pH min )휋 pH max −pH min ) 2 This formulation introduces complex nonlinearities and multimodal interactions across environmental axes, challenging models to inte- grate biological structure while avoiding rote memorization. F.3 Material Stress Behavior To evaluate symbolic models under realistic experimental conditions, this benchmark focuses on the stress-strain response of Aluminum 6061-T651 under thermal influence. The dataset records tensile strength measurements at six different temperatures, ranging from room temperature to 300°C, simulating diverse material states. In contrast to synthetic equations, this task lacks an explicit ground-truth formula, requiring data-driven inference of hidden physical laws. It serves as a test of a model’s capacity to identify empirical regularities in noisy, high-variance regimes. A widely used approximation of this behavior is given by: 휎 = ( 퐴 + 퐵휀 푛 ) ( 1 − ( 푇 − 푇 푟 푇 푚 − 푇 푟 ) 푚 ) where휎denotes stress,휀is strain,푇is the temperature,푇 푟 is a reference point, and푇 푚 is the melting point. The coefficients퐴,퐵,푛, and푚 are empirically determined for the alloy. This benchmark emphasizes the need for robust symbolic reasoning in the absence of prior symbolic templates, bridging theory-driven modeling with experimental data interpretation. F.4 LSR-Synth–Chemistry TheLSR-Synth–Chemistrydataset, introduced as part of the LLM-SRBench benchmark by Shojaee et al., is designed to evaluate symbolic regression models on chemically motivated yet synthetically constructed reaction kinetics. It comprises 36 differential equation discovery tasks, each modeling the time evolution of a reactant concentration퐴(푡)using novel, data-driven expressions. Each equation features a distinct combination of symbolic components—ranging from classic kinetic motifs (e.g., first- and second-order decay terms like−푘퐴(푡),−푘퐴(푡) 2 ) to synthetic nonlinearities. These include exponential decays such as exp(−푘 푠 푡), logarithmic forms like log(퐴(푡) + 1), square roots, and oscillatory terms such as sin(휔퐴(푡))and cos(⋅). In addition, rational expressions like 퐴(푡) 2 1+훽퐴(푡) 4 introduce challenges in terms of singularity avoidance and numerical stability. The dataset emphasizes symbolic diversity and parametric variability, making it a rigorous testbed for evaluating model generaliza- tion, interpretability, and robustness across structurally distinct formulations. Each task was carefully validated for analytical solvability, numerical stability, and scientific plausibility via expert review. G Additional Ablation Studies G.1 Ablation Study on Remaining Benchmarks We further validated the key components via ablations on the remaining three benchmarks, with the results shown in Figure8to Figure10. The experimental results are highly consistent with the findings presented in the main text. First, for the single-objective optimization variant (w/o MultiObj), a significant degradation in out-of-domain (OOD) generalization ability was observed across all benchmarks when compared to the full model. Although the in-domain (ID) fitting error of this variant is comparable to that of the full MOT-SR on some tasks, its performance on OOD data systematically confirms the tendency of single-objective optimization to overfit the ID data. This tendency consequently hinders the discovery of universally applicable scientific laws. Second, the ablation experiments on the internal components of the Meta Strategy Generator show that removing either the data analysis module (w/o Data) or the structure analysis module (w/o Struct) leads to a distinct decline in model performance. More critically, the degree of performance degradation caused by removing both modules simultaneously (w/o Strategy) exceeds the impact of removing either single module. This phenomenon clearly reveals the functional complementarity and synergistic effect between the data-driven prior knowledge and the structure-driven feedback mechanism. Both are indispensable for achieving MOT-SR’s final performance. G.2 Ablation Study on Data Analysis Tools To evaluate the contributions of the data analysis tools in MOT-SR, we conducted a series of ablation experiments. Our framework integrates six categories of analysis tools, and we set up six control groups (w/o tool-1 to 6), each corresponding to the removal of one category of tools. As shown in the Figure 11, removing any single category of analysis tools in this test leads to a decline in the model’s final performance. It is noteworthy that although there is only one tool in the sixth category, its removal still resulted in a final error nearly an order of magnitude higher than that of the full MOT-SR. This result strongly demonstrates that external analysis tools can effectively enhance a MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models large language model’s insight into the underlying structure and variable relationships within data, enabling it to capture deeper trends behind the data points and thereby generate more accurate symbolic equations. Full MOT-SR w/o Multi-Objectivew/o Structw/o Dataw/o Strategy 10 10 10 9 10 8 10 7 10 6 10 5 10 4 10 3 NMSE (log scale) 1.70e-10 2.16e-07 8.10e-09 2.32e-06 1.63e-06 7.72e-05 8.26e-09 3.36e-07 1.85e-06 9.01e-05 Oscillation 2 IDOOD Figure 8: Ablation study on Oscillation 2. Full MOT-SR w/o Multi-Objectivew/o Structw/o Dataw/o Strategy 10 2 10 1 10 0 NMSE (log scale) 2.08e-02 1.45e-02 3.06e-02 1.04e+00 7.21e-02 2.21e-02 6.37e-01 8.48e-01 7.24e-01 9.23e-01 E. coli Growth IDOOD Figure 9: Ablation study on E. coli growth. Wang et al. Full MOT-SR w/o Multi-Objectivew/o Structw/o Dataw/o Strategy 10 3 10 2 10 1 NMSE (log scale) 1.50e-02 2.04e-03 1.54e-02 6.04e-02 1.98e-02 2.02e-02 1.93e-02 1.21e-02 2.01e-02 2.33e-02 Stress-Strain IDOOD Figure 10: Ablation study on Stress-strain. Full MOT-SRw/o tool-1w/o tool-2w/o tool-3w/o tool-4w/o tool-5w/o tool-6 10 15 10 14 10 13 10 12 10 11 10 10 10 9 NMSE (log scale) 1.27e-15 3.13e-09 5.97e-11 1.38e-12 1.48e-12 5.55e-11 7.41e-14 Figure 11: Ablation study of the Data Analysis Tools on Oscillator 1. H HV and IGD Definition F.1 HV.The hypervolume (HV) measures the volume in the objective space that is dominated by the approximate Pareto front and bounded by a reference point. It reflects the convergence and diversity of the solutions. The HV is formally defined as follows: HV(푃, 푟 ∗ ) =VOL ( ⋃ v∈푃 [푣 1 , 푟 ∗ 1 ] × [푣 2 , 푟 ∗ 2 ] × ⋯ × [푣 푚 , 푟 ∗ 푚 ] ) , where푃denotes the approximate Pareto front obtained by the symbolic regression algorithm,v= (푣 1 , ... , 푣 푚 ) ⊤ represents an objective vector, VOL(⋅)indicates the Lebesgue measure, and푟 ∗ = (푟 ∗ 1 , ... , 푟 ∗ 푚 ) ⊤ is the reference point. To eliminate the impact of varying scales and units among objectives in HV calculation, we normalize each objective to the interval[0, 1]. For an objective value푓 푖 (푥), the normalized value is: 푓 ′ 푖 (푥) = 푓 푖 (푥) − 푧 ideal 푖 푧 nadir 푖 − 푧 ideal 푖 , where 푧 ideal 푖 =min v ∈푃 푣 푖 , 푧 nadir 푖 =max v∈푃 푣 푖 MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models After normalization, each objective is scaled to the range[0, 1]. Based on this, the reference point is set as푟 ∗ = (1.05, ... , 1.05) ⊤ . F.2 IGD.The Inverted Generational Distance (IGD) measures both the convergence and diversity of the predicted front by computing the distance from points on the true Pareto front to their nearest counterparts in the predicted set. The formula is as follows: IGD(푃, 푃 ∗ ) = 1 |푃 ∗ | ∑ p ∗ ∈푃 ∗ min p∈푃 푑(p,p ∗ ), where푃is the predicted Pareto front,푃 ∗ is a set of reference points sampled from the true Pareto front, and푑(p,p ∗ )represents the Euclidean distance between vectorspandp ∗ in objective space. The computation of the Inverted Generational Distance (IGD) metric conventionally depends on the availability of a known true Pareto front, which is not always accessible in practical scenarios. In the context of symbolic regression evaluation, however, the IGD metric benefits from an inherent advantage, as the true Pareto front can be explicitly defined by the ground truth equation푓 true . This allows for an exact assessment of the quality of generated equations relative to the ground truth, without reliance on an approximate Pareto front. Accordingly, the IGD formulation in this study is simplified as follows: IGD(푃, 푓 true ) =min 푓 ∈푃 푑(푓 , 푓 true ) Afterwards, all objectives are normalized to the[0, 1]range to ensure consistency with the true objective values. I Data Analysis Tools in MOT-SR I.1 Linear Correlation Tools Linear correlation tools provide foundational insights into how variables are related through linear or approximately linear patterns. These tools are instrumental for identifying primary dependencies, detecting dominant axes of variation, and establishing structural priors for symbolic equation modeling. MOT-SR incorporates four linear-correlation-based modules: •Pearson Correlation Coefficient: It quantifies the degree of linear association between two variables푥and푦using the classical Pearson correlation coefficient푟. The input consists of two real-valued vectors, the output includes the correlation score and the associated푝-value indicating statistical significance. A high absolute value of푟(close to 1) suggests a strong linear trend, while values near 0 indicate little to no linear dependency. This measure is especially suitable for Gaussian-like distributions and helps identify candidate variables for symbolic terms with additive or multiplicative linear effects. •Simple Linear Regression: This module fits a univariate linear model푦 = 푎 + 푏푥using least squares estimation and returns the regression coefficients, the푅 2 value (explained variance), and associated statistics such as푝-values and standard errors. The inputs are one-dimensional arrays푥and푦, and the output reflects how well a straight line fits the data. The푅 2 value is particularly informative, indicating how much of the variance in 푦can be explained by푥. It provides a predictive perspective on the relationship beyond correlation, supporting model selection based on explanatory power. •Residual Variance: This tool computes the variance of residuals푦 − ̂푦after applying linear regression. The lower the residual variance, the better the linear model captures the relationship between푥and푦. It takes the same inputs as Simple Linear Regression and internally reuses the linear regression output. This metric emphasizes prediction error dispersion, offering a direct measure of the noise or unmodeled nonlinear structure in the data. It is particularly valuable for pruning noisy or unstable variables from candidate model terms. •PCA Explained Variance: Principal Component Analysis (PCA) is applied to the joint space of푥and푦, and the proportion of variance explained by the first principal component is reported. The input is a two-dimensional matrix formed by stacking푥and푦, the output is a scalar value indicating the strength of the shared linear structure. This approach is robust to scaling and rotation and can capture the dominant direction of variation when the relationship is more general than univariate regression. It aids in identifying latent linear couplings and guides the selection of structurally informative variables. I.2 Nonlinear Dependency Tools Nonlinear dependency tools are employed to capture complex, non-monotonic interactions between variables, which are essential for accurately modeling nonlinear systems. Unlike linear tools, these methods can reveal relationships that are not adequately described by simple linear transformations. MOT-SR incorporates three nonlinear-dependency-based modules: •Spearman Rank Correlation:This tool measures the strength and direction of monotonic relationships between two variables using their rank values. It computes Spearman’s휌, a non-parametric counterpart to Pearson’s푟, by evaluating how well the relationship between 푥and푦can be described by a monotonic function. The input consists of two numerical vectors, the output includes the correlation Wang et al. coefficient and its푝-value. This method is particularly effective when data exhibit nonlinear but monotonic trends, such as saturating growth or sigmoid-like behavior. •Mutual Information:This tool quantifies the total amount of information shared between two variables, regardless of the specific functional form of their relationship. It estimates mutual information by discretizing the input variables푥and푦into bins and calculating the joint entropy. The output is a scalar metric representing dependency strength. Unlike correlation coefficients, mutual information can capture both linear and highly nonlinear dependencies, making it suitable for detecting complex statistical associations in symbolic modeling. •Mutual Information Regression Score:This variant employs thescikit-learnimplementation of mutual information to assess the relevance of푥in predicting푦using a regression-based formulation. The method internally estimates how much knowing푥reduces uncertainty about푦. It accepts one-dimensional arrays as input and outputs a scalar score. This score is robust to arbitrary nonlinearities and discontinuities, making it valuable for feature selection in nonlinear symbolic regression tasks. I.3 Time-Frequency Analysis Tools Time-frequency analysis tools are essential for detecting periodicities, oscillations, and transient dynamics that are characteristic of nonlin- ear and multi-scale systems. These tools enable symbolic regression models to incorporate periodic or time-varying terms where appropriate. MOT-SR employs two such modules: •Fast Fourier Transform Frequency Difference:This module applies the discrete Fourier transform to the input signals푥and푦 to extract their respective power spectra. It identifies the dominant frequency component for each variable and returns the absolute difference between them. The input consists of two one-dimensional arrays representing time-series data, and the output includes the dominant frequencies of both variables as well as their difference. This tool captures global periodic patterns and is useful for detecting whether the signals share synchronized or harmonically related structures. A small frequency difference suggests potential functional alignment via sinusoidal or oscillatory terms. •Wavelet Energy Correlation:This module leverages discrete wavelet decomposition to capture localized energy features of푥and푦at multiple temporal scales. It decomposes both signals into several levels using a specified wavelet basis (e.g., Daubechies 4), computes the energy at each level, and measures the Pearson correlation between the resulting energy vectors. The input includes two one-dimensional arrays and optional parameters for wavelet type and decomposition level. The output is a correlation coefficient representing how similarly the energy of the two signals is distributed across scales. Unlike Fourier-based methods, this tool excels at capturing transient and non- stationary dependencies, offering robust structural priors for equations involving local periodicity or bursts. I.4 Causal Inference Tools Causal inference tools aim to uncover directional relationships between variables, particularly whether the historical behavior of one variable contributes to or influences another. These tools help MOT-SR to build models that not only fit the data well but also respect temporal or structural causality, thereby improving interpretability and generalization. Two modules are implemented: •Granger Causality:This method assesses whether past values of a variable푥improve the prediction of a variable푦in a multivariate time series setting. The inputs are two equal-length time series arrays, and the test evaluates multiple lagged regression models to determine if푥Granger-causes푦. The output includes푝-values for different lag orders, and the smallest푝-value is used as the primary metric. A statistically significant result implies that the past of푥contains information predictive of푦, suggesting a directional dependency. This method is especially valuable for systems where delayed effects are prominent, such as in control dynamics or feedback loops. •Convergent Cross Mapping:CCM is a nonlinear causal discovery technique grounded in dynamical systems theory. It tests whether the state of variable푥can be reconstructed from the historical trajectory of푦, indicating that푥leaves an imprint on푦. The input consists of two time series and parameters for embedding dimension and library sampling. CCM constructs a manifold from the delay embedding of푦, uses it to cross-predict푥, and evaluates the reconstruction accuracy (typically using a cross-map skill휌). A high휌suggests that푥 causally influences푦. Unlike Granger causality, CCM does not rely on linearity or temporal precedence, making it ideal for identifying nonlinear, feedback-driven causality in complex systems. I.5 Dynamic Complexity Tools Dynamic complexity tools capture nuanced structural and temporal irregularities in time-series data, particularly useful for distinguishing chaotic, nonlinear, or asynchronous behaviors. MOT-SR leverages three such tools: •Lyapunov Exponent Difference:This tool estimates the maximal Lyapunov exponents (MLEs) of two time-series variables푥and푦, then computes their absolute difference. Lyapunov exponents characterize how sensitive a system is to initial conditions. A high MLE suggests chaotic behavior, whereas near-zero or negative values indicate stable dynamics. The input consists of two real-valued time series. The tool automatically determines the embedding dimension and time lag for phase space reconstruction. The resulting metric reflects how similarly (or differently)푥and푦behave in terms of dynamical predictability. MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models •Correlation Dimension Difference:This tool quantifies and compares the fractal (correlation) dimensions of푥and푦. The correlation dimension serves as a complexity measure, indicating how densely a system’s trajectory fills its phase space. The input includes two real- valued sequences and a fixed embedding dimension (usually2). The metric is the absolute difference between their estimated correlation dimensions. This captures differences in dynamic complexity and is useful for identifying structural mismatches in multivariate nonlinear processes. •Dynamic Time Warping Distance:DTW measures the alignment cost between푥and푦by allowing local nonlinear stretching or compression in time. Unlike Euclidean distance, DTW is robust to phase shifts and unequal pacing in temporal evolution. The input is a pair of univariate time series, and the output is a scalar DTW distance. This tool is particularly effective in identifying temporal patterns that share similar shapes but occur at different rates or phases. I.6 Distribution Consistency Tools This category includes statistical tools that assess whether two variables share similar probability distributions. Such tools are especially useful when validating whether a generated or transformed signal preserves the underlying distributional structure of the original data. •Kolmogorov–Smirnov Test:The Kolmogorov–Smirnov (KS) test is a non-parametric method that quantifies the maximum distance between the empirical cumulative distribution functions (ECDFs) of two datasets푥and푦. The input to this tool is a pair of real-valued vectors, and it returns the KS statistic (a scalar metric) along with the corresponding푝-value indicating statistical significance. The KS statistic captures how different the two distributions are—larger values indicate more significant deviations. This test is sensitive to both location and shape differences in the distributions and is thus useful for evaluating whether symbolic transformations (e.g., derived equations) maintain statistical fidelity to the data source. J Results Presentation J.1 Equation Recovery and Pareto Front Evolution Fig12ato Fig12dare sampled from four test cases in which MOT-SR successfully identified the ground-truth equations on the LSR-Synth– Chemistry dataset. Each figure illustrates the evolution of NMSE of the Pareto front equation set across iterations. For visualization, we present the structural skeletons of the equations discovered by MOT-SR. Highlighted segments in red indicate terms that match exactly with the corresponding ground-truth equation. It is noteworthy that the size of the Pareto front typically expands rapidly during the early stages, forming a diverse equation population. As iterations proceed, this population undergoes continuous updates. Once sufficient structural experience has accumulated, the Pareto front enters a phase of rapid convergence and eventually stabilizes to the correct target equation. Below we present comparisons between MOT-SR and LLM-SR on the same tasks: CRK14: Ground truth: −푘퐴(푡) + 푘 푝 sin(휔퐴(푡)) MOT-SR: 휃 0 퐴+휃 1 sin(휃 2 퐴) LLM-SR: 휃 0 퐴 ( 1 + (휃 1 퐴) 2 휃 2 ) 1 + 1+ 퐴 2 휃 5 1+ (휃 1 퐴) 2 휃 2 + 휃 4 퐴 휃 5 + 퐴 +휃 3 퐴 CRK19: Ground truth: −푘퐴(푡) 2 + 푘 푝 sin(휔퐴(푡)) MOT-SR: 휃 0 sin(휃 1 퐴)+휃 2 퐴 2 LLM-SR: −휃 0 퐴 −휃 1 퐴 2 + 휃 2 푡 + 휃 3 푡 2 + 휃 4 퐴 3 CRK21: Ground truth: −푘퐴(푡)푒 −푘 푠 푡 + 푘 푝 sin(휔퐴(푡)) MOT-SR: 휃 0 퐴푒 −휃 1 푡 +휃 2 sin(휃 3 퐴) Wang et al. LLM-SR: 휃 0 푡 +휃 1 퐴푒 −휃 2 푡 + 휃 3 퐴 휃 4 푡 + 1 + 휃 5 퐴 휃 6 1 + 휃 7 퐴 휃 6 CRK36: Ground truth: −푘퐴(푡) + 푘 푞 퐴(푡)log(훾푡 + 1) MOT-SR: 휃 0 퐴+휃 1 퐴log(휃 2 푡 + 1) LLM-SR: 휃 0 퐴+ 휃 1 푒 −휃 2 푡 퐴 + 휃 3 퐴 휃 4 + 퐴 + 휃 5 퐴 1 + 휃 6 퐴 + 휃 7 푡 +휃 8 퐴+ 휃 9 푡퐴 These results clearly demonstrate that MOT-SR is capable of precisely capturing the exact components of the ground-truth equations, yielding structurally closer expressions than LLM-SR. We attribute this improvement to MOT-SR’s enhanced exploratory capability derived from multi-objective evaluation, as well as its strategy-guided generation process. 050100150200250 Iterations 10 12 10 10 10 8 10 6 10 4 Nmse Test All (log scale) NMSE Spread over Iterations Min-Max Range Min NMSE eq1. eq3. eq2. (a) Pareto Front Performance in CRK14. 050100150200 Iterations 10 12 10 10 10 8 10 6 10 4 10 2 10 0 Nmse Test All (log scale) NMSE Spread over Iterations Min-Max Range Min NMSE eq1. eq2. (b) Pareto Front Performance in CRK19. 050100150200 Iterations 10 11 10 9 10 7 10 5 10 3 10 1 Nmse Test All (log scale) NMSE Spread over Iterations Min-Max Range Min NMSE eq1. eq2. (c) Pareto Front Performance in CRK21. 050100150200250 Iterations 10 13 10 11 10 9 10 7 10 5 10 3 Nmse Test All (log scale) NMSE Spread over Iterations Min-Max Range Min NMSE eq.1 eq.3 eq.2 eq.4 (d) Pareto Front Performance in CRK36. Figure 12: Pareto Front Performance results across four benchmarks. J.2 Comparison of convergence dynamics. Figure13presents the wall-clock time evolution of NMSE (log scale) for MOT-SR and the baseline LLM-SR across four representative benchmarks. The results demonstrate that MOT-SR consistently achieves substantially faster error reduction, while LLM-SR converges more slowly and plateaus at significantly higher error levels. This advantage is particularly pronounced in the Oscillation tasks, where MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models 0481216 Time (hours) 10 13 10 10 10 7 10 4 10 1 NMSE (log scale) Oscillation 1 0481216 Time (hours) 10 11 10 9 10 7 10 5 10 3 10 1 NMSE (log scale) Oscillation 2 0481216 Time (hours) 10 2 10 1 10 0 NMSE (log scale) E. coli Growth 0481216 Time (hours) 10 2 10 1 NMSE (log scale) Stress-Strain MOT-SR (Llama-3.1)LLM-SR (Llama-3.1) Figure 13: Wall-clock time evolution of NMSE showing faster convergence of MOT-SR over LLM-SR. MOT-SR rapidly identifies accurate dynamical structures, but is also evident in more challenging domains such asE. coligrowth and stress– strain modeling, where the baseline stagnates prematurely. These findings highlight that the integration of tool-guided analysis and multi- objective optimization enables MOT-SR to explore the hypothesis space more efficiently, yielding superior equations under comparable computational budgets. Wang et al. J.3 Evolution of Discovered Equations 025050075010001250150017502000 Iteration 10 14 10 12 10 10 10 8 10 6 10 4 10 2 10 0 Normalized MSE c 1 x + c 0 v + c 2 ve c 3 x c 1 sin(c 2 x) + c 4 cos(c 5 x) + c 6 x + c 0 v + c 3 ve c 7 x + c 8 xv c 0 x + c 1 vc 2 x 2 + c 3 v 2 + c 4 x 3 + c 5 v 3 c 6 xv + c 7 sin(c 8 x) + c 9 v c 0 x + c 1 v + c 2 x 3 + c 3 v 3 c 4 xv + c 5 sin(c 6 x) Figure 14: Iterative evolution of equation structures under MOT-SR (Oscillation1). MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models 02505007501000 Iteration 10 14 10 12 10 10 10 8 10 6 10 4 10 2 10 0 Normalized MSE c 0 t 2 A + c 1 tA + c 2 A c 3 A 2 + c 4 t + c 5 t 2 + c 6 c 0 t 2 + c 1 A + c 2 e c 3 t + c 4 sin(c 5 t)A + c 6 tA A ( c 0 + c 1 e c 2 t + c 3 t c 4 + c 5 t 2 ) Ac 0 + Ac 1 log(c 2 t + 1) Figure 15: Iterative evolution of equation structures under MOT-SR (CRK36). Evolution of Equation Structures for Oscillation 1.Figure14presents the iterative evolution of equation structures for the problem os- cillation 1. In the early iterations, the generated expressions contained a wide range of nonlinear and polynomial terms, many of which represented exploratory hypotheses rather than meaningful components. Guided by the TOOL-AUGMENTED VARIABLE ANALYSIS mod- ule, the framework progressively refined these structures. When the analytical tools within this module (e.g., wavelet corr, FFT, and PCA mapping) revealed periodic features, sinusoidal terms such as sine and cosine functions were introduced into the candidate equations, thereby complementing the polynomial basis and improving alignment with the observed data. During the iteration, structurally inconsistent or statistically unsupported terms were systematically removed, while the retained non- linear and oscillatory components captured the essential dynamics of the system. The process ultimately converged toward compact for- mulations that balance accuracy, generalization, and complexity. This refinement demonstrates how the iterative mechanism, augmented by TOOL-AUGMENTED VARIABLE ANALYSIS, effectively distills the equation space and produces models consistent with both empirical behavior and theoretical plausibility. Evolution of Equation Structures for CRK36.Figure15presents the iterative evolution of equation structures for the problem CRK36. At the beginning of the process, the candidate functions contained a variety of nonlinear combinations involving both the temporal variable and A. Under the guidance of two correlation coefficients, the framework progressively refined these structures by eliminating inconsistent nonlinearities and retaining only the dominant dependencies. As the iterations advanced, the generated functions converged toward a fully linear dependence on A. K Gaussian Noise Evaluation To evaluate the robustness of MOT-SR in noisy environments, we conduct experiments on theOscillator1benchmark by injecting Gaussian noise with standard deviations휎 = 0.001, 0.002into the training data. Model performance is assessed under both in-domain (ID) and out-of-domain (OOD) test conditions. Unlike conventional symbolic regression approaches that rely solely on fitting accuracy, MOT-SR integrates structural priors—derived from analytical tools—and performs multi-objective evaluation that jointly considers accuracy, complexity, and generalization. This allows it to maintain high-fidelity symbolic recovery even when data are corrupted. Table 4shows that MOT-SR consistently outperforms LLMSR across all metrics, achieving significantly lower NMSE and higher Acc avg-0.01 scores under both noise settings. These gains stem from MOT-SR’s integration of tool-informed structural priors via the Meta Strategy Gen- erator and its multi-objective Equation Generator, which jointly enable the model to identify robust, generalizable expressions even under noisy supervision. Wang et al. Table 4: Comparison of MOT-SR and LLMSR onOscillator1under Gaussian noise (휎 = 0.001, 0.002). Metrics include NMSE and accuracy (Acc avg-0.01 ) on ID and OOD test sets. All models use LLaMA-3.1 as backbone. Model 휎 = 0.001휎 = 0.002 IDOODIDOOD NMSE↓Acc avg-0.01 ↑NMSE↓Acc avg-0.01 ↑NMSE↓Acc avg-0.01 ↑NMSE↓Acc avg-0.01 ↑ LLMSR1.59푒-589.20%0.002110.46%2.68푒-582.08%0.002510.50% MOT-SR 5.43e-1199.99%5.49e-798.79%3.98e-899.74%8.57e-797.88% L Assessment of Python AST Length as a Complexity Proxy Classical symbolic regression methods operate directly on mathematical expression trees, where structural complexity is commonly mea- sured by counting expression-tree nodes. MOT-SR allows the Equation Generator to produce executable Python functions without a fixed operator dictionary or function whitelist. This program-level representation does not directly provide a canonical mathematical expression tree. Known syntactic bias.In our implementation, each candidate function is parsed withast.parse, and its Python AST length is defined as the total number of nodes returned byast.walk. This metric reflects the syntactic structure of Python code and introduces a known bias. For example, a composite mathematical call such asnp.sin(x)can produce more Python AST nodes than a basic arithmetic operation, even when the two have comparable mathematical expression-tree complexity. Python AST length is thus treated as a computational proxy for mathematical expression complexity. Computational trade-off.A stricter alternative would convert each generated Python function into a normalized mathematical expres- sion, construct the corresponding expression tree, and count its nodes. This procedure aligns more closely with conventional complexity measures in symbolic regression, but adds a conversion step for every candidate generated during the search. MOT-SR evaluates thousands of candidates, making repeated conversion a substantial computational burden. Python AST length can be computed deterministically and directly from the generated source code with negligible overhead. Empirical ranking consistency.Pareto-based selection depends primarily on the relative complexity ordering of candidate equations. We therefore assess whether Python AST length preserves the ordering induced by mathematical expression-tree complexity. Let 퐶 AST (푓 ) denote the Python AST node count of a candidate equation푓, and let 퐶 Expr (푓 ) denote the node count of its normalized mathematical expression tree. For a set풬of candidate pairs, we define the pairwise ranking- disagreement ratio as 푟 dis = 1 |풬| ∑ (푓 푖 ,푓 푗 )∈풬 핀 [ sign ( 퐶 AST (푓 푖 ) − 퐶 AST (푓 푗 ) ) ≠sign ( 퐶 Expr (푓 푖 ) − 퐶 Expr (푓 푗 ) ) ] ,(23) where comparisons tied under either metric are excluded. We retrospectively computed both complexity measures for valid candidate equations collected from the search logs and evaluated nearly 400,000candidate pairs. The resulting pairwise disagreement ratio was approximately5%, indicating that the two metrics produced the same relative ordering for approximately95%of the evaluated pairs. Within the candidate distributions observed in our experiments, Python AST length therefore preserves the expression-tree ordering in most comparisons. Because the complexity objective is used to determine relative dominance and construct the Pareto front, this level of agreement indicates that the proxy has a limited effect on the practical selection order. Composite-function robustness.We further examine whether the additional syntactic cost of composite functions systematically sup- presses candidates containing trigonometric or other non-polynomial terms. The analysis covers tasks whose target structures depend strongly on such functions, including the Oscillation 1 and Oscillation 2 benchmarks and trigonometry-heavy tasks in LSR-Synth–Chemistry. MOT-SR discovers compact equations containing terms such as sin(⋅)and cos(⋅), while achieving strong numerical and symbolic recovery performance on these tasks. These results indicate that the additional AST cost assigned to composite calls does not prevent the correspond- ing structures from entering or remaining on the Pareto front. Python AST length is a language-dependent and syntax-dependent measure, rather than an exact representation of mathematical com- plexity. We use it as a low-overhead engineering proxy that supports flexible Python-based equation generation. Within the candidate distributions examined in this work, the retrospective ranking analysis and the results on composite-function-heavy tasks show that this proxy provides a sufficiently consistent complexity ordering for Pareto-based search. MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models M Prompt J.1 Prompt Design.The prompts used to construct the meta strategy are shown in Fig.16and Fig.17. Specifically, Fig.16presents the meta prompt used during initialization, while Fig.17corresponds to the prompts designed for iterative optimization. Together, these constitute the prompt inputs to the Meta Strategy Generator. J.2 Examples.The prompts used by the Equation Generator for the Oscillators 1 task are illustrated in Fig.18to Fig.21. The prompts received by the Meta Strategy Generator for the same task are shown in Fig.22to Fig.27. Fig.28provides examples of data analysis tool usage, drawn from tool invocations in the CRK36 task. Notably, through the integration of external tools, the language model successfully identified a strong linear relationship between the variable푥 0 (that is,퐴(푡)) and the target variable푦. As demonstrated in the final output, MOT-SR successfully recovered the ground-truth governing equation for the CRK36 problem. Wang et al. SYSTEM ROLE: You are a **scientist-assistant LLM** supervising a symbolic-regression pipeline. You must analyze existing equations, select appropriate analysis tools, and produce guidance for a generator LLM. ## Allowed analysis tools allowed_list Choose **up to 3** tools—exclusively from the list above—based **solely on the diagnostic needs revealed by the residual data you receive**.If a desired tool is not listed, omit it. ## RESPONSE FORMAT (STRICT) Return exactly one JSON object (no markdown fences). It must contain: 1. "structure_insight_prompt_for_generator_LLM" –either a string, or an object with `guidance.prompt_for_generator_LLM` (string). 2. "analysis_tools" –a JSON array of 1–5 tool names from the allowed list. Escape newlines inside strings as ` `. After the closing brace ``, output nothing else. Meta Prompt Figure 16: Meta Prompt. ## Context context ## Mathematical Pattern Analysis Task Analyze the following equations to identify **common mathematical patterns and shared structural elements**: .join(blocks) ## Available Analysis Tools tool_list ## Pattern Analysis Instructions 1. **Examine Mathematical Structures**: Look across all equations for: -Recurring mathematical functions (trigonometric, exponential, polynomial) -Similar variable interaction patterns (x*y, x+y, x^n, etc.) -Consistent coefficient magnitudes or sign patterns -Common denominators or fraction structures -Similar nesting or grouping of terms 2. **Identify Shared Elements**: Find mathematical components that appear in multiple high-performing equations 3. **Interpret Mathematical Significance**: Explain what these patterns might indicate about: -Underlying physical relationships -Mathematical constraints or principles -Optimal solution characteristics 4. **Generate Guidance**: Create a detailed prompt that will guide future equation generation based on successful patterns 5. **Select Analysis Tools**: Choose up to 5 data analysis tools that could help validate the patterns you identified ## RESPONSE FORMAT Return ONLY JSON with keys: structure_insight_prompt_for_generator_LLM, analysis_tools. -structure_insight_prompt_for_generator_LLM: Detailed guidance based on your pattern analysis for generating better equations -analysis_tools: Array of up to 5 tool names from the available list Focus on mathematical commonalities and their significance in equation generation. Refinement Prompt MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models ## Context context ## Mathematical Pattern Analysis Task Analyze the following equations to identify **common mathematical patterns and shared structural elements**: .join(blocks) ## Available Analysis Tools tool_list ## Pattern Analysis Instructions 1. **Examine Mathematical Structures**: Look across all equations for: -Recurring mathematical functions (trigonometric, exponential, polynomial) -Similar variable interaction patterns (x*y, x+y, x^n, etc.) -Consistent coefficient magnitudes or sign patterns -Common denominators or fraction structures -Similar nesting or grouping of terms 2. **Identify Shared Elements**: Find mathematical components that appear in multiple high-performing equations 3. **Interpret Mathematical Significance**: Explain what these patterns might indicate about: -Underlying physical relationships -Mathematical constraints or principles -Optimal solution characteristics 4. **Generate Guidance**: Create a detailed prompt that will guide future equation generation based on successful patterns 5. **Select Analysis Tools**: Choose up to 5 data analysis tools that could help validate the patterns you identified ## RESPONSE FORMAT Return ONLY JSON with keys: structure_insight_prompt_for_generator_LLM, analysis_tools. -structure_insight_prompt_for_generator_LLM: Detailed guidance based on your pattern analysis for generating better equations -analysis_tools: Array of up to 5 tool names from the available list Focus on mathematical commonalities and their significance in equation generation. Refinement Prompt Figure 17: Prompts Designed for Iterative Optimizationt. You are a helpful assistant tasked with discovering mathematical function structures for scientific systems. Complete the 'equation' function below, considering the physical meaning and relationships of inputs. # Scientific Analysis Based on Pareto Front Selection ## Core Generation Guidance ## Core Generation Guidance Generate equations that incorporate the identified mathematical patterns and shared structural elements. Emphasize leveraging variable interactions, consistent coefficient magnitudes, and common denominators. Prioritize equations that respect physical constraints and optimization characteristics. ## Generation Instructions Based on the above analysis, focus on: 1. Incorporating the identified recurring mathematical functions 2. Using the successful variable interaction patterns 3. Maintaining coefficient relationships that work well 4. Respecting the physical constraints and relationships 5. Building upon the shared mathematical elements Oscillator1 Figure 18: Example Prompt in Oscillator1 - For Equation Generator. Wang et al. ## Data Analysis Results Applied 3 analysis tools to all variable pairs: **Tools used**: spearman_corr, residual_var, corr **Dataset overview**: -2 independent variables analyzed -10000 total data points ## Correlation Analysis: x0_y **Strong correlation detected** (r=-0.920) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be negative -as one increases, the other decreases. ## Residual Variance Analysis: x0_y **Residual variance**: 0.000419 Low residual variance suggests the linear relationship captures most variation. ## Correlation Analysis: x0_y **Strong correlation detected** (r=-0.951) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be negative -as one increases, the other decreases. ## Correlation Analysis: x1_y **Weak correlation detected** (r=0.103) The relationship between x1_y and output appears weak or non-linear. ## Residual Variance Analysis: x1_y **Residual variance**: 0.004359 Low residual variance suggests the linear relationship captures most variation. ## Correlation Analysis: x1_y **Weak correlation detected** (r=-0.032) The relationship between x1_y and output appears weak or non-linear. Figure 19: Example Prompt in Oscillator1 - For Equation Generator. ## Generation Strategy Based on the analysis of Pareto front samples and data relationships: 1. **Build upon success**: Use the selected high-quality solution as foundation 2. **Targeted improvement**: Focus on reducing residuals in identified problem areas 4. **Data-driven insights**: Incorporate findings from variable relationship analysis 5. **Variable interactions**: Consider the significant correlations and patterns identified 6. **Maintain balance**: Preserve the multi-objective balance that led to Pareto selection 7. **Consider trade-offs**: Balance accuracy, generalization, and computational efficiency """ Find the mathematical function skeleton that represents acceleration in a damped nonlinear oscillator system with driving force, given data on position, and velocity. """ MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models ## Generation Strategy Based on the analysis of Pareto front samples and data relationships: 1. **Build upon success**: Use the selected high-quality solution as foundation 2. **Targeted improvement**: Focus on reducing residuals in identified problem areas 4. **Data-driven insights**: Incorporate findings from variable relationship analysis 5. **Variable interactions**: Consider the significant correlations and patterns identified 6. **Maintain balance**: Preserve the multi-objective balance that led to Pareto selection 7. **Consider trade-offs**: Balance accuracy, generalization, and computational efficiency """ Find the mathematical function skeleton that represents acceleration in a damped nonlinear oscillator system with driving force, given data on position, and velocity. """ Figure 20: Example Prompt in Oscillator1 - For Equation Generator. def equation(x: np.ndarray, v: np.ndarray, params: np.ndarray) -> np.ndarray: """ Mathematical function for acceleration in a damped nonlinear oscillator Args: x: A numpyarray representing observations of current position. v: A numpyarray representing observations of velocity. params: Array of numeric constants or parameters to be optimized Return: A numpyarray representing acceleration as the result of applying the mathematical function to the inputs. """ # Driving force considering correlation with x0_y and output f = params[0] * np.exp(-params[1] * v**2) # Driving force with damping # Incorporating strong linear relationship between x0_y and output linear_term= params[2] * x # Linear coefficient for damping # Damping effect considering correlation with x1_y and output damp = params[3] * (1 + 0.5 * np.sin(params[4] * x)) # Damping coefficient with trigonometric term # Nonlinear oscillator equation considering correlations a = -(damp + 0.5 * linear_term) * v # Additional non-linear terms for better model fit b = params[5] * np.sin(params[6] * x + params[7] * v) # Nonlinear term with hybrid sine function c = params[8] * v**2 + f * np.sin(params[9] * x) # Nonlinear term with exponenialtrig function # Regularization technique to prevent overfitting return a + b + c def equation_v1(x: np.ndarray, v: np.ndarray, params: np.ndarray) -> np.ndarray: """Improved version of `equation_v0`.""" Figure 21: Example Prompt in Oscillator1 - For Equation Generator. Wang et al. SYSTEM ROLE: You are a **scientist-assistant LLM** supervising a symbolic-regression pipeline. You must analyze existing equations, select appropriate analysis tools, and produce guidance for a generator LLM. ## Allowed analysis tools ccm_causality, corr, corr_dim_relation, dtw_distance, fft_cross_freq, granger_causality, ks_test_diff, lin_reg, lyapunov_relation, mutual_info, mutual_info_regression_score, pca_mapping, pearson_corr, residual_var, spearman_corr, wavelet_corr Choose **up to 3** tools—exclusively from the list above—based **solely on the diagnostic needs revealed by the residual data you receive**.If a desired tool is not listed, omit it. ## RESPONSE FORMAT (STRICT) Return exactly one JSON object (no markdown fences). It must contain: 1. "structure_insight_prompt_for_generator_LLM" –either a string, or an object with `guidance.prompt_for_generator_LLM` (string). 2. "analysis_tools" –a JSON array of 1–5 tool names from the allowed list. Escape newlines inside strings as ` `. After the closing brace ``, output nothing else. ## Context Symbolic regression problem: oscillator1. Dataset has 2 input variables and 1 target variable. Training data shape: (10000, 2) -> (10000,) Analyzing equations selected from Pareto front with residual data. ## Mathematical Pattern Analysis Task Analyze the following equations to identify **common mathematical patterns and shared structural elements**: Oscillator1 Figure 22: Example Prompt in Oscillator1 - For Meta Strategy Generator. ### Equation 1 Expression : """Mathematical function for acceleration in a damped nonlinear oscillator Args: x: A numpyarray representing observations of current position. v: A numpyarray representing observations of velocity. params: Array of numeric constants or parameters to be optimized params[0]: Spring constant, effectively the driving force coefficient params[1]: Damping coefficient params[2]: Amplitude, the maximum value of driving force params[3]: Angular frequency, c=driving force frequency params[4]: Mass of object, m params[5]: Friction or driving force damping coefficient params[6]: Non-linear damping exponent params[7]: Non-linear driving force frequency millers params[8]: Driving force input frequency peak frequency millers params[9]: Non-linear displacement spring constant magnetic field coefficients Returns: A numpyarray representing acceleration as the result of applying the mathematical function to the inputs. """ # Driving force c = params[0] # Original driving force coefficient A = params[2] # Original amplitude # Frictionless surface parameters (assuming frictionless) m = params[4] # mass of object # Non-linear ODE (On-Off Oscillators) # x-dependent, non-linear spring constant k = params[9] + params[1] * np.exp(-np.abs(v)) # Original spring constant # Incorporating driving force non-linearity # exponential riding f = c * (A + np.tanh(params[6] * (v + params[7]))) * np.sin(params[8] * x) # Non-linear ODE application return -m * (k * v + f) MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models ### Equation 1 Expression : """Mathematical function for acceleration in a damped nonlinear oscillator Args: x: A numpyarray representing observations of current position. v: A numpyarray representing observations of velocity. params: Array of numeric constants or parameters to be optimized params[0]: Spring constant, effectively the driving force coefficient params[1]: Damping coefficient params[2]: Amplitude, the maximum value of driving force params[3]: Angular frequency, c=driving force frequency params[4]: Mass of object, m params[5]: Friction or driving force damping coefficient params[6]: Non-linear damping exponent params[7]: Non-linear driving force frequency millers params[8]: Driving force input frequency peak frequency millers params[9]: Non-linear displacement spring constant magnetic field coefficients Returns: A numpyarray representing acceleration as the result of applying the mathematical function to the inputs. """ # Driving force c = params[0] # Original driving force coefficient A = params[2] # Original amplitude # Frictionless surface parameters (assuming frictionless) m = params[4] # mass of object # Non-linear ODE (On-Off Oscillators) # x-dependent, non-linear spring constant k = params[9] + params[1] * np.exp(-np.abs(v)) # Original spring constant # Incorporating driving force non-linearity # exponential riding f = c * (A + np.tanh(params[6] * (v + params[7]))) * np.sin(params[8] * x) # Non-linear ODE application return -m * (k * v + f) Figure 23: Example Prompt in Oscillator1 - For Meta Strategy Generator. Wang et al. ID Score : -0.00013475 (in-distribution performance) OOD Score : -5.38816e-05 (out-of-distribution performance) Eval Time : 10.6924s (computational efficiency) Avg Accuracy : -9.43156e-05 (overall accuracy) Gen Gap : 8.08681e-05 (generalization stability) Efficiency : 0.0855254 (time-normalized performance) **In-Distribution (ID) Dataset Analysis:** Worst Region Statistics: -Sample Count: 520 samples -Threshold: 0.000958914 -Mean Residual: -0.0008611 -Mean Abs Residual: 0.00112713 -Max Abs Residual: 0.00148016 -Std Residual: 0.000750307 Input Variables in Worst Region: x0: [-0.481813, -0.419474] x1: [-0.124159, 0.113473] Target Variable (y) in Worst Region: y: [0.048175, 0.0951525] **Out-of-Distribution (OOD) Dataset Analysis:** Worst Region Statistics: -Sample Count: 480 samples -Threshold: 0.000740516 -Mean Residual: -0.000603334 -Mean Abs Residual: 0.000866861 -Max Abs Residual: 0.000997573 -Std Residual: 0.000629749 Input Variables in Worst Region: x0: [-0.520795, 0.489177] x1: [-0.0900758, 0.282136] Target Variable (y) in Worst Region: y: [-0.0614007, 0.060305] Figure 24: Example Prompt in Oscillator1 - For Meta Strategy Generator. MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models ### Equation 2 Expression : """ Mathematical function for acceleration in a damped nonlinear oscillator. Args: x: A numpyarray representing observations of current position. v: A numpyarray representing observations of velocity. params: Array of numeric constants or parameters to be optimized params[0]: Spring constant (k), driving force coefficient (omega) params[1]: Damping coefficient (gamma) params[2]: Amplitude (A) params[3]: Angular frequency (omega) Returns: A numpyarray representing acceleration as the result of applying the mathematical function to the inputs. """ # Parameters for driving force omega = params[0] # Angular frequency omega_driving_force= omega # Parameters for friction or driving force damping gamma = params[1] # Damping coefficient gamma_friction_driving_damping= gamma # Mass and amplitude for non-linear oscillator m = params[4] # Mass A_driving_damping= params[2] # Amplitude # Calculate velocity induced force force = omega_driving_force* np.sin(params[3] * np.tanh(params[5] * v)) # Acceleration equation # Consider strong correlations between position and velocity acceleration = -m * (force + gamma_friction_driving_damping* v) return acceleration Figure 25: Example Prompt in Oscillator1 - For Meta Strategy Generator. ID Score : -1.16205 (in-distribution performance) OOD Score : -1.00194 (out-of-distribution performance) Eval Time : 0.698804s (computational efficiency) Avg Accuracy : -1.082 (overall accuracy) Gen Gap : 0.160103 (generalization stability) Efficiency : 0.58865 (time-normalized performance) **In-Distribution (ID) Dataset Analysis:** Worst Region Statistics: -Sample Count: 520 samples -Threshold: 0.101862 -Mean Residual: 0.0907549 -Mean Abs Residual: 0.109957 -Max Abs Residual: 0.113557 -Std Residual: 0.0622102 Input Variables in Worst Region: x0: [-0.512046, 0.444764] x1: [0.11177, 0.212819] Target Variable (y) in Worst Region: y: [-0.104672, 0.110646] **Out-of-Distribution (OOD) Dataset Analysis:** Worst Region Statistics: -Sample Count: 480 samples -Threshold: 0.0978116 -Mean Residual: 0.0446015 -Mean Abs Residual: 0.103459 -Max Abs Residual: 0.111418 -Std Residual: 0.0934521 Input Variables in Worst Region: x0: [-0.562765, 0.495206] x1: [0.0432888, 0.221801] Target Variable (y) in Worst Region: y: [-0.103036, 0.10915] Wang et al. ID Score : -1.16205 (in-distribution performance) OOD Score : -1.00194 (out-of-distribution performance) Eval Time : 0.698804s (computational efficiency) Avg Accuracy : -1.082 (overall accuracy) Gen Gap : 0.160103 (generalization stability) Efficiency : 0.58865 (time-normalized performance) **In-Distribution (ID) Dataset Analysis:** Worst Region Statistics: -Sample Count: 520 samples -Threshold: 0.101862 -Mean Residual: 0.0907549 -Mean Abs Residual: 0.109957 -Max Abs Residual: 0.113557 -Std Residual: 0.0622102 Input Variables in Worst Region: x0: [-0.512046, 0.444764] x1: [0.11177, 0.212819] Target Variable (y) in Worst Region: y: [-0.104672, 0.110646] **Out-of-Distribution (OOD) Dataset Analysis:** Worst Region Statistics: -Sample Count: 480 samples -Threshold: 0.0978116 -Mean Residual: 0.0446015 -Mean Abs Residual: 0.103459 -Max Abs Residual: 0.111418 -Std Residual: 0.0934521 Input Variables in Worst Region: x0: [-0.562765, 0.495206] x1: [0.0432888, 0.221801] Target Variable (y) in Worst Region: y: [-0.103036, 0.10915] Figure 26: Example Prompt in Oscillator1 - For Meta Strategy Generator. ## Available Analysis Tools ccm_causality, corr, corr_dim_relation, dtw_distance, fft_cross_freq, granger_causality, ks_test_diff, lin_reg, lyapunov_relation, mutual_info, mutual_info_regression_score, pca_mapping, pearson_corr, residual_var, spearman_corr, wavelet_corr ## Pattern Analysis Instructions 1. **Examine Mathematical Structures**: Look across all equations for: -Recurring mathematical functions (trigonometric, exponential, polynomial) -Similar variable interaction patterns (x*y, x+y, x^n, etc.) -Consistent coefficient magnitudes or sign patterns -Common denominators or fraction structures -Similar nesting or grouping of terms 2. **Identify Shared Elements**: Find mathematical components that appear in multiple high-performing equations 3. **Interpret Mathematical Significance**: Explain what these patterns might indicate about: -Underlying physical relationships -Mathematical constraints or principles -Optimal solution characteristics 4. **Generate Guidance**: Create a detailed prompt that will guide future equation generation based on successful patterns 5. **Select Analysis Tools**: Choose up to 5 data analysis tools that could help validate the patterns you identified ## RESPONSE FORMAT Return ONLY JSON with keys: structure_insight_prompt_for_generator_LLM, analysis_tools. -structure_insight_prompt_for_generator_LLM: Detailed guidance based on your pattern analysis for generating better equations -analysis_tools: Array of up to 5 tool names from the available list Focus on mathematical commonalities and their significance in equation generation. Figure 27: Example Prompt in Oscillator1 - For Meta Strategy Generator. MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models ## Correlation Analysis: x0_y **Strong correlation detected** (r=1.000) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be positive -as one increases, the other increases. pearson_corr: ## Correlation Analysis: x0_y **Strong correlation detected** (r=1.000) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be positive -as one increases, the other increases. spearman_corr: ## Mutual Information Regression Analysis: x1_y **Mutual information score**: 8.294050 Mutual information regression measures the dependency between variables in a regression context. mutual_info: ## Residual Variance Analysis: x0_y **Residual variance**: 0.057119 Low residual variance suggests the linear relationship captures most variation. residual_var: ## Linear Regression Analysis: x1_y **Strong linear relationship** (R²=0.996) Linear equation: output ≈ 0.000 ×x1_y + 0.000 Consider incorporating linear terms in your equation. lin_reg: Partial Data Analysis Tools Usage in CRK36: Figure 28: Examples of Data Analysis Tools Usage.