MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models Boxiao Wang
arXiv:2607.29561v1 [cs.LG] 31 Jul 2026
National Key Laboratory of Cognition and Decision Intelligence for Complex Systems Institute of Automation, Chinese Academy of Sciences Beijing, China [email protected]
Runxiang Wang
Chongming Li School of Astronomy and Space Science, University of the Chinese Academy of Sciences Beijing, China [email protected]
Zhiwei Chen
Yifan Zhang
Jian Cheng
Relativistic Astrophysics, Institute for Theoretical Physics Goethe University Frankfurt Frankfurt am Main, Germany [email protected]
National Key Laboratory of Cognition and Decision Intelligence for Complex Systems Institute of Automation, Chinese Academy of Sciences Beijing, China [email protected]
National Lab of Pattern Recognition Institute of Automation, Chinese Academy of Sciences Beijing, China [email protected]
Abstract Discovering compact and interpretable equations from observational data is essential for understanding complex scientific systems. 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 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-SR1 ), 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)
1
Kai Li
University of the Chinese National Key Laboratory of Academy of Sciences Cognition and Decision School of Advanced Intelligence for Complex Interdisciplinary Sciences Systems Beijing, China Institute of Automation, [email protected] Chinese Academy of Sciences Beijing, China [email protected]
Code is available at https://github.com/wswbx/MOT-SR.
orbital modeling, an important problem in space-based gravitationalwave 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.
Keywords Symbolic Regression, Extreme Mass-Ratio Inspirals, Gravitational Wave Modeling, Large Language Models
1 Introduction Symbolic Regression (SR) [34] aims to discover underlying mathematical 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 modeling 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], motivating 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
Wang et al.
Iterative Loop (Generation 𝑡 → 𝑡 + 1)
Observed Data 𝑋, 𝑦
Meta Strategy Generator
Equation Generator & Evaluation
Residual Calculation
Pareto Frontier Buffer ( 𝒫𝓉)
Scientific Toolbox
𝑒1 , 𝑒2 , … , 𝑒𝑁
Diversity-Aware Parent Selection
LLM Generation Correlation FFT/Wavelet Causal
Data Analysis Prompting
LLM Reasoning
Population of Equations
Search Strategy
Update Frontier ( 𝒫𝓉+1)
Multi-Objective Evaluation NMSE
(train-ID)
NMSE
(train-OOD)
AST Length (Complexity)
Non-Dominated Sorting & Pruning
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 strategy. Following this search strategy, an Equation Generator & Evaluation module performs diversity-aware parent selection, generates new equations, and evaluates them by NMSE ID and 𝐷 OOD (training-derived in-domain (ID)/out-ofon 𝐷train train domain (OOD) split) and by AST length for complexity. Nondominated 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 LLMSR [52] and LaSR [18] leverage LLM’s in-context learning capabilities and scientific priors to perform symbolic reasoning and equation generation. Despite encouraging progress, existing LLM-based SR methods face two key limitations. First, they typically lack systematic analysis 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 error, while overlooking other critical factors such as equation complexity, generalization, and diversity. This can lead to overfitting and premature convergence to locally optimal solutions. To address these challenges, we propose Multi-Objective Toolaugmented Symbolic Regression (MOT-SR), a unified framework (Figure 1) inspired by the human scientific modeling process. Scientists typically begin by analyzing data using a variety of tools, evaluate candidate hypotheses from multiple perspectives, and iteratively refine their modeling direction accordingly. MOT-SR operationalizes 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 variable relationships from multiple complementary dimensions, such as linearity, periodicity, and causality. These insights are converted into interpretable priors and structural constraints that guide equation generation. In parallel, MOT-SR adopts a multi-objective evaluation mechanism that jointly assesses candidate equations across three dimensions: accuracy on the in-domain (ID) training subset
ID , generalization on the out-of-domain (OOD) training subset 𝐷train 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 manner, 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 Equation Generator & Evaluation module performs diversity-aware parent selection and generates new candidate equations, which are evaluated and integrated into the Pareto-frontier buffer via nondominated sorting and pruning. We evaluate MOT-SR on five benchmarks spanning physics, chemistry, 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, generalization, and equation compactness. Beyond benchmark evaluation, we apply MOT-SR to extreme mass-ratio inspiral (EMRI) orbital evolution, a challenging problem in space-based gravitational-wave astronomy where small local errors accumulate over many orbital cycles. Using 58 configurations for equation discovery and 30 held-out configurations for a posteriori evaluation, MOT-SR identifies a compact symbolic correction that achieves the lowest trajectory-level integration error— approximately three orders of magnitude lower than a neural residual baseline and 26.8 times lower than LLM-SR. These results demonstrate that MOT-SR can recover transferable correction structures that remain reliable across unseen physical configurations and long horizons, highlighting the potential of interpretable equation discovery 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𝑖 ∈ ℝ𝑑 ,
𝑦𝑖 ∈ ℝ,
where x𝑖 denotes a 𝑑 -dimensional input vector and 𝑦𝑖 is the corresponding scalar output. The goal is to discover an explicit analytical equation 𝑓 (⋅) such that the predicted outputs 𝑦𝑖̂ = 𝑓 (x𝑖 ) closely approximate the ground truth 𝑦𝑖 . To assess the quality of a candidate equation, the normalized mean squared error (NMSE) is defined as 𝑛
NMSE(𝑓 , 𝐷) =
2
𝑓 (x𝑖 ) − 𝑦𝑖 1 ∑( ) . 𝑛 𝑖=1 𝜎𝑦
(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 accuracy, SR also values simplicity and generalization, seeking equations 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:
each variable conditioned on others. Together, these tools guide whether nonlinear or higher-order terms should be introduced. Time-frequency analysis tools help detect periodicity and transient dynamics, which frequently occur in oscillatory and mul• Equation skeleton generation: Structured prompts contain tiscale systems. MOT-SR employs two complementary methods: task-specific information (e.g., variable names, optimizaFast Fourier Transform [8] identifies dominant global frequency tion goals, example equations), guiding the LLM to procomponents, while the wavelet transform energy spectrum [11] capduce physically plausible equation skeletons. tures localized, non-stationary fluctuations. These insights support • Parameter optimization: The skeletons’ free parameters are the inclusion of periodic terms (e.g., sin) in candidate equations. optimized (e.g., via BFGS [14] ) and scored using NMSE. Causal inference tools are used to identify whether one vari• Feedback: High-quality equations are retained and reused able may influence another in a predictive or explanatory sense. as in-context examples, enabling iterative refinement through MOT-SR adopts Granger causality [16] for detecting temporal causal feedback-driven generation. influence in linear time-series, and Convergent Cross Mapping [55] While promising, LLM-based SR methods face key limitations: for identifying latent causality in nonlinear systems with potential they lack systematic modeling of variable dependencies, leading to delays. These methods provide structural signals that enhance the structurally under-informed equations; and they rely on a singleinterpretability and explanatory power of generated equations. objective evaluation focused solely on fitting error, neglecting facDynamic complexity tools help assess intrinsic system richtors such as complexity and generalization. These shortcomings ness and redundancy, guiding the pruning of over-specified comdiminish their effectiveness in solving more challenging tasks. ponents in equations. MOT-SR employs three methods: the Lyapunov exponent [61], which measures sensitivity to initial condi3 Method tions and indicates chaotic behavior; the correlation dimension [17], To address the above limitations, we propose MOT-SR (Multi-Objective which estimates the system’s effective degrees of freedom; and DyTool-augmented Symbolic Regression), a unified framework that namic Time Warping (DTW) [48], which evaluates time-shifted simintegrates external data analysis tools, multi-objective evaluation, ilarity between variable trajectories. Together, they offer structural and cooperative LLMs to enhance equation quality and search efcues for constructing compact and robust equations. ficiency. MOT-SR first extracts structural priors by analyzing variDistribution consistency tools evaluate whether input variable relationships with diverse analytical tools. It then introduces a ables behave uniformly across different input regions. MOT-SR multi-objective evaluation mechanism that jointly considers fitting uses the Kolmogorov–Smirnov (KS) test [35] to detect distributional error, equation complexity, and generalization, dynamically mainshifts between subdomains, informing the use of region-dependent taining a Pareto front through non-dominated sorting. Finally, two structures to reflect local variations in the data distribution. complementary LLMs work in tandem: one generates structural Rather than introducing new tools, the key innovation of MOTrefinement strategies, while the other synthesizes candidate equaSR lies in enabling LLMs to autonomously invoke and coordinate tions accordingly, forming a closed-loop process that continuously these tools to extract dependency patterns among variables. These guides and improves equation discovery. insights are distilled into concise guidance that informs variable selection and function composition, thereby improving responsive3.1 Tool-Augmented Variable Analysis ness to data characteristics and enhancing the scientific plausibilThis component constructs structural priors by quantifying diverse ity of generated equations. By adaptively combining outputs from variable relationships using a suite of carefully designed analytical heterogeneous analyses, for instance by linking correlation meatools. Details are provided in Appendix I. sures with periodicity detection, MOT-SR supports more targeted Linear correlation tools are essential for identifying domiequation discovery. The impact of this tool-augmented analysis on nant variables and constructing interpretable regression structures. equation generation is substantiated by the case studies in AppenMOT-SR integrates several complementary methods to assess lindix J.3. Looking ahead, MOT-SR opens the possibility for LLMs to ear dependencies from multiple statistical angles. The Pearson corsynthesize new tools, further expanding the scope of SR research. relation coefficient [42] quantifies pairwise linear associations, especially effective for Gaussian-like data. Simple linear regression 3.2 Multi-Objective Evaluation and residual variance analysis [38] evaluate predictive capacity and To enhance search efficiency and model quality, MOT-SR adopts error stability. PCA-based explained variance [19] identifies the key a multi-objective evaluation scheme that jointly considers predicdirections of structural variance. These tools jointly establish a tive accuracy, generalization, and structural simplicity. A Paretosolid basis for linear trend detection. based selection strategy maintains a diverse set of non-dominated Nonlinear dependency tools are employed to capture comcandidate equations, improving robustness and exploration of the plex interactions essential for modeling nonlinear systems. MOTsolution space. We define the ID/OOD regions within the training SR integrates three complementary methods: the Spearman rank set 𝐷train as a spatial split (Figure 2). The detailed configuration of correlation [54] measures monotonic associations based on rank, the ID-OOD partition for 𝐷train is provided in Appendix C. offering robustness to noise and non-Gaussian data; mutual information [50] measures the overall statistical dependency between variables without assuming any parametric form; and mutual inEvaluation Metrics. Unlike prior LLM-based SR methods that optiformation regression [43] quantifies the marginal contribution of mize only fitting error, MOT-SR employs a three-fold metric:
Wang et al.