InsightSR: Refining Symbolic Regression Search Spaces via Parallel Semantic and Structural LLM Guidance

# InsightSR: Refining Symbolic Regression Search Spaces via Parallel Semantic and Structural LLM Guidance
Source: [https://arxiv.org/html/2608.25291](https://arxiv.org/html/2608.25291)
###### Abstract

> Symbolic regression \(SR\) seeks to discover parsimonious mathematical laws from observational data, yet conventional approaches often struggle with the vast combinatorial search space of physically meaningful expressions\. We present InsightSR, a framework that embeds Large Language Models \(LLMs\) as a guiding layer around the PySR genetic programming engine\. Rather than relying on LLMs to generate expressions directly, InsightSR uses LLMs to progressively transform the search space itself through two complementary pathways: a Semantic Seed Pathway that proposes dimensionally consistent functional skeletons, and a Structural Feature Pathway that recommends nonlinear feature transformations\. These transformations accumulate over iterations, broadening the input space and shifting the symbolic search from constructing deep expression trees over raw variables to assembling shallow trees over a rich, semantically informed feature set\. A post\-generation feedback loop evaluates candidates, categorizes features by their empirical utility, and refines the guidance for the next iteration, transforming the discovery process from open\-ended generation into iterative, self\-correcting refinement\. Across three benchmarks, InsightSR achieves a 95% exact recovery rate on the Feynman benchmark and 80\.18% accuracy on the LLM\-SRBench LSR\-Transform task, substantially outperforming state\-of\-the\-art genetic programming and neural\-symbolic methods while maintaining strong out\-of\-distribution generalization on real\-world datasets\.

## Introduction

Symbolic regression \(SR\) aims to distill underlying mathematical laws directly from observed data, serving as a cornerstone for automated scientific discovery\[schmidt2009distilling,brunton2016discovering\]\. While SR produces interpretable closed\-form expressions, its optimization is inherently NP\-hard\[virgolin2022symbolic\]and the combinatorial search space of candidate expressions grows exponentially with expression complexity\. A further challenge is the difficulty of incorporating physical constraints, such as dimensional consistency and known conservation laws, into the search process\. As a result, purely data\-driven SR methods frequently yield expressions that achieve high numerical accuracy but violate basic physical principles\.

Existing approaches predominantly rely on stochastic search over the discrete space of mathematical expressions\. The most widely used are Genetic Programming \(GP\) variants\[koza1994genetic,cranmer2023pysr,virgolin2021improving\], which evolve populations of expression trees through mutation and crossover, and more recent methods based on Monte Carlo Tree Search \(MCTS\)\[sun2022symbolic,shojaee2023transformer\], which frame discovery as a sequential decision\-making task\. Although these methods can recover approximate symbolic forms, they suffer from low sample efficiency and tend to overfit under observational noise\. More critically, without mechanisms to enforce physical constraints, they frequently produce expressions that are mathematically flexible but physically meaningless, lacking the parsimony and domain consistency necessary for scientific interpretation\.

The integration of neural networks has further expanded the SR toolkit\. Deep reinforcement learning methods such as DSR\[petersen2021deep\]and uDSR\[landajuela2022unified\]treat equation discovery as a policy optimization problem, using learned gradients to guide the search through expression space\. Transformer\-based models such as NeSymReS\[biggio2021neural\]and E2E\[kamienny2022end\]leverage large\-scale pre\-training to map data patterns directly to mathematical expressions\. Other specialized approaches, such as PhySO\[landajuela2021discovering\], incorporate dimensional constraints into the search to enforce physical consistency\. Although these neural methods accelerate discovery, they are fundamentally data\-driven: their performance depends heavily on the training distribution, and they often fail to generalize to out\-of\-distribution regimes\. Moreover, the physical constraints they enforce are typically hard\-coded during training rather than dynamically reasoned about, limiting their adaptability to novel physical contexts\.

Recent work has explored integrating Large Language Models \(LLMs\) into SR to leverage scientific priors encoded\. LLM\-SR\[shojaee2025llmsr\]generates equation skeletons via an LLM and optimizes their parameters through BFGS, storing successful patterns in an experience buffer to guide future iterations\. SR\-LLM\[guo2025srllm\]casts SR as a reinforcement learning problem with an LLM\-based policy network, coupled with retrieval\-augmented generation for incremental knowledge accumulation\. LaSR\[li2024lasr\]focuses on LLM\-guided feature engineering, recommending transformations to expand the input space\. PiSR via LLM\[taskin2026pisr\]incorporates LLM\-derived evaluations into the loss function as a physics\-informed regularization term\. While promising, these methods are often constrained by two limitations\. First, they incur substantial inference overhead, as they require frequent LLM queries throughout the search process\. Second, and more importantly, each method employs LLM guidance along a single dimension—either structural hypothesis generation or input space transformation—without exploiting the potential synergy of integrating both\. This separation overlooks a key opportunity: structural priors can inform which features to engineer, and engineered features can simplify the structural forms needed to fit the data\.

To address these limitations, we propose InsightSR, a framework that embeds LLMs as a guiding layer around the PySR genetic programming engine\[cranmer2023pysr\]\. Rather than relying on the LLM to generate candidate expressions directly, InsightSR uses the LLM to progressively transform the search space through two complementary pathways\. A Semantic Seed Pathway proposes dimensionally consistent functional skeletons that provide a physics\-informed warm\-start, while a Structural Feature Pathway recommends nonlinear feature transformations that accumulate over generations, progressively broadening the input space\. This shifts the burden of the search from constructing deep expression trees over raw variables to assembling shallow combinations over a semantically enriched feature set\. A learning\-from\-results feedback loop closes each generation: the LLM evaluates the Pareto\-optimal candidates, categorizes features by their empirical utility, and updates a knowledge base that informs both pathways in the subsequent iteration, transforming symbolic discovery from open\-ended generation into iterative, self\-correcting refinement\.

Our work makes three key contributions:

1. 1\.Semantic Seed Pathway\.We introduce a mechanism that leverages domain\-aware LLMs to propose dimensionally consistent functional skeletons that serve as physics\-informed seeds for the evolutionary search\. By auditing and reconfiguring candidate topologies to align with target physical units, this pathway prunes the search space of physically inconsistent expressions before numerical fitting begins\. The LLM also provides per\-operator complexity biases that steer the genetic programming engine toward physically plausible operator combinations\.
2. 2\.Structural Feature Pathway\.We design a complementary pathway in which the LLM recommends nonlinear feature transformations based on the problem context and the performance history of prior generations\. These transformations accumulate across iterations, progressively broadening the input space\. This shifts the evolutionary search from constructing deep expression trees over raw variables to assembling shallow combinations over an enriched feature basis, substantially reducing the structural complexity that the symbolic engine must resolve\.
3. 3\.Closed\-Loop Iterative Refinement\.We implement a learning\-from\-results mechanism that closes the discovery loop\. After each generation, the LLM evaluates candidates along multiple dimensions including numerical accuracy, physical interpretability, and feature utility\. These assessments are accumulated in a dynamic knowledge base that informs both the semantic and structural pathways in subsequent iterations, enabling the search to self\-correct and progressively converge toward the optimal symbolic form\.

## Method

![Refer to caption](https://arxiv.org/html/2608.25291v1/figures/aa.png)Figure 1:The system architecture of InsightSR\. The framework follows a four\-phase iterative loop per generation: \(1\)Contextual Initialization, where the LLM performs automated unit synthesis and domain analysis; \(2\)Parallel Guidance Strategy, which concurrently generates seed expressions \(Strategy A\) and recommends feature enhancements \(Strategy B\); \(3\)Search with PySR, where the integrated heuristics steer the PySR engine; and \(4\)Strategic Analysis, where the LLM evaluates the top\-KKcandidates and formulates a plan for the next iteration\.The goal of SR is to discover an analytical expressionf:ℝd→ℝf:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}that accurately describes the relationship in a dataset𝒟=\{\(𝐱i,yi\)\}i=1N\\mathcal\{D\}=\\\{\(\\mathbf\{x\}\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{N\}\. Following the standard GP formulation, we treat this as a multi\-objective optimization that balances numerical accuracy against structural parsimony:

minf∈ℱ⁡ℒ⁡\(f,𝒟\)=MSE​\(y,f⁡\(𝐗\)\)\+λ⋅𝒞⁡\(f\),\\min\_\{f\\in\\mathcal\{F\}\}\\mathcal\{L\}\(f,\\mathcal\{D\}\)=\\text\{MSE\}\(y,f\(\\mathbf\{X\}\)\)\+\\lambda\\cdot\\mathcal\{C\}\(f\),\(1\)
where𝒞⁡\(f\)\\mathcal\{C\}\(f\)denotes the structural complexity offf\(measured by the number of nodes in its expression tree\) andλ\\lambdacontrols the parsimony penalty\.

InsightSR organizes the discovery process into a four\-phase iterative loop, summarized in Figure[1](https://arxiv.org/html/2608.25291#Sx2.F1)and formalized in Algorithm[1](https://arxiv.org/html/2608.25291#alg1)\. Each generation embeds LLM\-derived domain knowledge at three intervention points—before, during, and after the evolutionary search—shifting PySR from undirected stochastic exploration toward purposefully guided optimization\. We detail each phase below\.

Algorithm 1InsightSR Search Procedure0:Dataset

\(𝐗,y\)\(\\mathbf\{X\},y\), problem context

𝒞\\mathcal\{C\}, units

𝒰\\mathcal\{U\}
0:Optimal symbolic expression

f∗f^\{\*\}
Initialize knowledge base

𝒦←∅\\mathcal\{K\}\\leftarrow\\emptyset
for

g=1g=1to

GGdo

// Dual\-Path Heuristic Guidance

ℱs​e​e​d←LLM\.SemanticSeeds​\(𝒞,𝒰,𝒦\)\\mathcal\{F\}\_\{seed\}\\leftarrow\\text\{\\sc LLM\.SemanticSeeds\}\(\\mathcal\{C\},\\mathcal\{U\},\\mathcal\{K\}\)⊳\\trianglerightSeed Pathway

𝐗a​u​g←LLM\.StructuralFeatures​\(𝒞,𝒦\)\\mathbf\{X\}\_\{aug\}\\leftarrow\\text\{\\sc LLM\.StructuralFeatures\}\(\\mathcal\{C\},\\mathcal\{K\}\)⊳\\trianglerightFeature Pathway

𝒢←MergeGuidance​\(ℱs​e​e​d,𝐗a​u​g\)\\mathcal\{G\}\\leftarrow\\text\{MergeGuidance\}\(\\mathcal\{F\}\_\{seed\},\\mathbf\{X\}\_\{aug\}\)

candidates←PySR\.GuidedSearch​\(𝐗a​u​g,y,𝒢\)\\text\{candidates\}\\leftarrow\\text\{\\sc PySR\.GuidedSearch\}\(\\mathbf\{X\}\_\{aug\},y,\\mathcal\{G\}\)

eval,insights←LLM\.StrategicAnalysis​\(candidates\)\\text\{eval\},\\text\{insights\}\\leftarrow\\text\{\\sc LLM\.StrategicAnalysis\}\(\\text\{candidates\}\)

𝒦←KB\.Accumulate​\(insights,eval\)\\mathcal\{K\}\\leftarrow\\text\{\\sc KB\.Accumulate\}\(\\text\{insights\},\\text\{eval\}\)

if

Loss​\(f∗\)<10−10\\text\{Loss\}\(f^\{\*\}\)<10^\{\-10\}then

break

endif

endfor

returnBest expression

f∗f^\{\*\}from candidates

### Contextual Initialization and Unit Synthesis

The framework first resolves the physical dimensions of each variable\. It extracts available dimensional information from the dataset metadata and problem description\. For variables whose units remain ambiguous or unspecified, an LLM infers probable dimensions based on the semantic context of the problem and the physical role of the target variable\. The resulting unit assignments persist as dimensional constraints throughout the discovery loop\. By enforcing dimensional homogeneity, the framework prunes physically inconsistent candidates from the search space, ensuring that the symbolic search remains grounded in the underlying physics\.

### Parallel Guidance Strategy

The core of InsightSR is a parallel guidance strategy in which the LLM simultaneously contributes along two complementary dimensions: top\-down semantic seeding and bottom\-up feature engineering\.

Semantic Seed Pathway\.This pathway leverages domain\-aware LLMs to propose dimensionally consistent functional skeletons that serve as physics\-informed seeds for the evolutionary search\. Given the problem metadata and the unit constraints resolved in Section 3\.1, the LLM generates candidate expression topologies whose dimensional consistency is audited before they enter the population\. Each skeleton contains symbolic constants that are subsequently optimized numerically\. By seeding the population with physically plausible structures, this pathway prunes a vast region of dimensionally inconsistent expressions from the search space before the evolutionary optimization begins\.

Structural Feature Pathway\.In parallel, this pathway enriches the input space through LLM\-guided feature engineering\. At each generation, the LLM analyzes the performance history and identifies high\-utility nonlinear transformationsϕ⁡\(𝐱\)\\phi\(\\mathbf\{x\}\)based on patterns observed in surviving candidates\. These transformations are appended to the input matrix:

𝐗a​u​g=\[𝐱1,…,𝐱d,ϕ1​\(𝐱\),…,ϕk​\(𝐱\)\],\\mathbf\{X\}\_\{aug\}=\[\\mathbf\{x\}\_\{1\},\\dots,\\mathbf\{x\}\_\{d\},\\phi\_\{1\}\(\\mathbf\{x\}\),\\dots,\\phi\_\{k\}\(\\mathbf\{x\}\)\],\(2\)whereϕj​\(𝐱\)\\phi\_\{j\}\(\\mathbf\{x\}\)includes operations such as power\-law termsxinx\_\{i\}^\{n\}and transcendental mappingssin⁡\(xi\)\\sin\(x\_\{i\}\)\. Crucially, these features accumulate across generations, progressively broadening the input space\. This shifts the search from constructing deep expression trees to discover complex nonlinear relationships to assembling shallow combinations over an already\-enriched feature set, substantially reducing the structural depth that the evolutionary engine must resolve\.

The two pathways converge within the PySR engine\. Rather than searching the unconstrained spaceℱ\\mathcal\{F\}over raw inputs𝐱\\mathbf\{x\}, the optimizer operates over an informed manifold:

f∗=arg⁡minf∈ℱs​e​e​d⁡ℒ⁡\(f⁡\(𝐗a​u​g\),𝒟\)\.f^\{\*\}=\\arg\\min\_\{f\\in\\mathcal\{F\}\_\{seed\}\}\\mathcal\{L\}\(f\(\\mathbf\{X\}\_\{aug\}\),\\mathcal\{D\}\)\.\(3\)Semantic seeds provide macro\-level structural constraints grounded in physics, while engineered features supply micro\-level building blocks derived from empirical patterns\. Together, they confine the search to a physically plausible region of the expression space without sacrificing the numerical flexibility needed for high\-precision fitting\.

### Search with PySR

The synthesized guidance consists of three components: symbolic seeds, operator preferences, and the augmented feature matrix𝐗a​u​g\\mathbf\{X\}\_\{aug\}\. These are passed to a modified PySR engine for evolutionary search\. To steer the search toward LLM\-informed functional structures, we introduce a complexity\-biasing mechanism that adjusts the structural penalty of each expression tree\. The adjusted complexity𝒞′​\(f\)\\mathcal\{C\}^\{\\prime\}\(f\)is computed as the sum of baseline operator costscoc\_\{o\}modified by an LLM\-provided biasωo\\omega\_\{o\}:

𝒞′​\(f\)=∑o∈f\(co\+ωo\)\.\\mathcal\{C\}^\{\\prime\}\(f\)=\\sum\_\{o\\in f\}\(c\_\{o\}\+\\omega\_\{o\}\)\.\(4\)Operators recommended by the LLM receive a negative biasωo\\omega\_\{o\}, granting them a complexity discount that promotes their survival on the Pareto front\. Conversely, operators deemed physically implausible are penalized with a positive bias, increasing their effective cost and suppressing their propagation across generations\. In addition, we employ a warm\-start mechanism that retains the best candidate from the previous generation as an elite member of the new population\. This ensures that prior refinements are preserved and built upon, enabling progressive convergence toward the optimal symbolic form\.

### Multi\-Dimensional Evaluation and Strategic Feedback

After each evolutionary search completes, the LLM serves as a post\-generation critic to evaluate the Pareto\-optimal candidates and guide the next iteration\. The evaluation synthesizes multiple criteria: numerical accuracy, physical interpretability, and the completeness with which the expression incorporates the relevant variables\. In parallel, the system performs a diagnostic analysis of the feature space, categorizing each engineered transformation as effective or ineffective based on its contribution to reducing the residual error\. These assessments are stored in a dynamic knowledge base, enabling the framework to learn from both successful patterns and failure modes across generations\. The knowledge base then informs the next iteration in three ways: it refines the seed expressions proposed by the Semantic Seed Pathway, adjusts the operator complexity biases used by the PySR engine, and prioritizes feature transformations for the Structural Feature Pathway\. By closing the loop between empirical search and semantic reasoning, the framework progressively narrows the search toward physically meaningful expressions until convergence\.

## Experimental Setup

### Benchmarks

We evaluate on three benchmarks summarized in Table[1](https://arxiv.org/html/2608.25291#Sx3.T1)\.

Table 1:Benchmark datasets overviewBenchmarkProblemsExampleFeynman\[udrescu2020ai\]100m01−v2/c2\\frac\{m\_\{0\}\}\{\\sqrt\{1\-v^\{2\}/c^\{2\}\}\}LLM\-SRBench\[shojaee2025llm\]240−k​A​\(t\)2\+kz​A​\(t\)2β​A​\(t\)4\+1\-kA\(t\)^\{2\}\+\\frac\{k\_\{z\}A\(t\)^\{2\}\}\{\\beta A\(t\)^\{4\}\+1\}Real\-World\[shojaee2025llmsr\]40\.8sin⁡\(x\)−0\.5​v3−0\.2​x3−0\.5​x​v−x​cos⁡\(x\)\\sin\(x\)\-0\.5v^\{3\}\-0\.2x^\{3\}\-0\.5xv\-x\\cos\(x\)TheFeynman 100benchmark\[udrescu2020ai\]consists of 100 physics equations derived from the Feynman Lectures on Physics, covering classical mechanics, electromagnetism, quantum mechanics, and thermodynamics\. Each problem has a known ground\-truth expression, enabling exact recovery evaluation\.

TheLLM\-SRBench\[shojaee2025llm\]comprises five categories: LSR\-Transform and LSR\-Synth across Physics, Chemistry, Biology, and Material Science\. Each category features different data\-generating mechanisms and varying complexity, assessing cross\-domain generalization\.

TheReal\-Worldbenchmark\[shojaee2025llmsr\]includes four real\-world datasets \(Oscillator 1, Oscillator 2, E\. coli, Stress\-Strain\)\. Each dataset provides both in\-distribution \(ID\) and out\-of\-distribution \(OOD\) test splits, enabling evaluation of out\-of\-distribution generalization\.

### Implementation Details

Table[2](https://arxiv.org/html/2608.25291#Sx3.T2)lists the key hyperparameters used in our experiments\.

Table 2:HyperparametersParameterValueLLMQwen\-3\.5 27BLLM Temperature0\.3PySR Iterations100–500LLM Generations3–30Populations30Population Size50Parsimony0\.003Max Expression Size20–25Early Stop Threshold10−1010^\{\-10\}Loss FunctionMSE
### Evaluation Metrics

The Coefficient of Determination \(R2R^\{2\}\) measures the proportion of variance captured by the model, defined as:

R2=1−∑i=1n\(yi−y^i\)2∑i=1n\(yi−y¯\)2R^\{2\}=1\-\\frac\{\\sum\_\{i=1\}^\{n\}\(y\_\{i\}\-\\hat\{y\}\_\{i\}\)^\{2\}\}\{\\sum\_\{i=1\}^\{n\}\(y\_\{i\}\-\\bar\{y\}\)^\{2\}\}\(5\)whereyiy\_\{i\}represents the ground truth,y^i\\hat\{y\}\_\{i\}is the predicted value, andy¯\\bar\{y\}is the empirical mean\. Consistent with our optimization objective, the Normalized Mean Squared Error \(NMSE\) is defined asNMSE=1−R2\\text\{NMSE\}=1\-R^\{2\}\.

The Accuracy under Tolerance \(accτ\\text\{acc\}\_\{\\tau\}\) is defined as a binary success metric for each problem\. A discovery is considered successful \(11\) only if the maximum relative error across allnntest samples is within the thresholdτ\\tau; otherwise, it is marked as failure \(00\):

accτ=𝕀⁡\(max1≤i≤n⁡\|yi−y^iyi\|≤τ\)\\text\{acc\}\_\{\\tau\}=\\mathbb\{I\}\\left\(\\max\_\{1\\leq i\\leq n\}\\left\|\\frac\{y\_\{i\}\-\\hat\{y\}\_\{i\}\}\{y\_\{i\}\}\\right\|\\leq\\tau\\right\)\(6\)where𝕀⁡\(⋅\)\\mathbb\{I\}\(\\cdot\)is the indicator function\. The Exact Recovery Rate is defined as a special case where the symbolic form is mathematically equivalent to the ground truth, typically resulting inaccτ→1\\text\{acc\}\_\{\\tau\}\\to 1forτ→0\\tau\\to 0\.

## Results

### Feynman Benchmark

Figure[2](https://arxiv.org/html/2608.25291#Sx4.F2)summarizes the performance on the Feynman benchmark\. InsightSR achieves a 95% exact recovery rate and an averageR2R^\{2\}of 0\.9999, outperforming the previous state of the art\. Figure[3](https://arxiv.org/html/2608.25291#Sx4.F3)shows the Pareto front of loss versus expression complexity\. The median loss declines consistently as complexity increases, and the narrowing interquartile range at higher complexity levels indicates enhanced numerical stability\. We further evaluated robustness to data perturbations\. As shown in Table[3](https://arxiv.org/html/2608.25291#Sx4.T3), accuracy remains high under moderate noise, with Acc0\.1reaching 79% at a noise level of 0\.001\.

Figure 2:Feynman benchmark: exact recovery rate \(bars\) with average R2overlay \(diamonds\)\.![Refer to caption](https://arxiv.org/html/2608.25291v1/figures/loss_convergence.png)Figure 3:Loss vs\. expression complexity Pareto front\. The solid line shows the median; shaded region shows the 25th–75th percentile\.Table 3:Noise vs Accuracy PerformancenoiseAcc0\.1\(%\)0\.05370\.01680\.005720\.00179097
### LLM\-SRBench

We further evaluate our method on the LLM\-SRBench across five scientific domains\. As summarized in Table[4](https://arxiv.org/html/2608.25291#Sx4.T4), our method achieves a substantial gain in the LSR\-Transform category, reaching an Acc0\.1of 80\.18% compared to the previous best of 50\.45% \(LaSR\)\.

Performance is consistent across the four domain\-specific synthetic datasets\. Specifically, our method achieves an accuracy of 66\.67% in Chemistry, 50\.00% in Biology, 40\.91% in Physics, and 92\.00% in Material Science\. In terms of numerical precision, our method yields NMSE values that are several orders of magnitude lower than other LLM\-based solvers\. These results are further illustrated in Figure[4](https://arxiv.org/html/2608.25291#Sx4.F4), which shows the relative performance across all methods\.

Table 4:LLM\-SRBench results\.MethodLSR\-TransformChemistryBiologyPhysicsMaterial Sci\.Acc0\.1↑\\uparrowNMSE↓\\downarrowAcc0\.1↑\\uparrowNMSE↓\\downarrowAcc0\.1↑\\uparrowNMSE↓\\downarrowAcc0\.1↑\\uparrowNMSE↓\\downarrowAcc0\.1↑\\uparrowNMSE↓\\downarrowDirect Prompting6\.310\.26313\.882\.2e\-24\.160\.4659\.090\.0650\.000\.048SGA8\.110\.23216\.665\.5e\-412\.510\.0139\.090\.05136\.116\.0e\-4LaSR50\.450\.00138\.929\.1e\-520\.831\.5e\-431\.819\.9e\-472\.049\.2e\-6LLM\-SR39\.640\.00952\.774\.1e\-629\.163\.1e\-636\.367\.6e\-588\.283\.2e\-9Ours80\.181\.3e\-1466\.678\.0e\-750\.007\.3e\-740\.911\.6e\-592\.004\.7e\-8Figure 4:LLM\-SRBench Acc0\.1comparison across five categories\.
### Real\-World Datasets

Table[5](https://arxiv.org/html/2608.25291#Sx4.T5)compares the NMSE results of our method against several baselines on four real\-world datasets, covering both in\-distribution \(ID\) and out\-of\-distribution \(OOD\) scenarios\. Overall, our method obtains the lowest NMSE in most test cases\.

On the Oscillator 1 and Oscillator 2 datasets, our method achieves ID NMSE values of9\.55×10−119\.55\\times 10^\{\-11\}and2\.45×10−92\.45\\times 10^\{\-9\}, respectively, representing an improvement of several orders of magnitude over the previous best LLM\-based method \(LLM\-SR\)\. For the E\. coli dataset, our method yields an ID NMSE of1\.76×10−31\.76\\times 10^\{\-3\}and an OOD NMSE of1\.41×10−21\.41\\times 10^\{\-2\}\. In the Stress\-Strain task, our results are2\.00×10−22\.00\\times 10^\{\-2\}\(ID\) and5\.32×10−25\.32\\times 10^\{\-2\}\(OOD\), which are comparable to or better than LLM\-SR\. These results are further illustrated in Figure[5](https://arxiv.org/html/2608.25291#Sx4.F5)\.

Table 5:Real\-world datasets results\. ID: in\-distribution, OOD: out\-of\-distribution\.MethodOscillator 1Oscillator 2E\. coliStress\-StrainIDOODIDOODIDOODIDOODGPlearn1\.55e\-025\.57e\-017\.55e\-013\.19e\+001\.08e\+001\.04e\+001\.06e\-014\.09e\-01NeSymReS4\.70e\-035\.38e\-012\.49e\-016\.47e\-01—1\.45e\+007\.93e\-016\.38e\-01E2E8\.20e\-033\.72e\-011\.40e\-011\.91e\-016\.32e\-011\.45e\+002\.26e\-015\.87e\-01DSR8\.70e\-032\.45e\-015\.80e\-021\.95e\-019\.45e\-012\.43e\+003\.33e\-011\.11e\+00uDSR3\.00e\-047\.00e\-043\.20e\-031\.50e\-033\.32e\-015\.46e\+005\.02e\-021\.76e\-01PySR9\.00e\-043\.11e\-012\.00e\-049\.80e\-033\.76e\-021\.01e\+003\.31e\-021\.30e\-01LLM\-SR4\.65e\-075\.00e\-042\.12e\-073\.81e\-052\.14e\-022\.64e\-022\.10e\-025\.16e\-02Ours9\.55e\-111\.15e\-042\.45e\-091\.03e\-051\.76e\-031\.41e\-022\.00e\-025\.32e\-02Figure 5:LLM\-SR benchmark NMSE comparison\.
### Qualitative Case Study

To provide an intuitive understanding of our method’s performance, we compare the ground\-truth equations against the discovered expressions across all three benchmarks\. Table[Qualitative Case Study](https://arxiv.org/html/2608.25291#Sx4.SSx4)highlights how our method recovers the core functional forms, with matching terms bolded in the discovered column\.

Table 6:Ground truth vs\. discovered equations\. Matching terms are bolded\.BenchmarkProblemGround TruthDiscoveredFeynmanI\.6\.2ae−θ2/22​π\\displaystyle\\frac\{e^\{\-\\theta^\{2\}/2\}\}\{\\sqrt\{2\\pi\}\}0\.3989⋅𝐞−θ𝟐/𝟐0\.3989\\cdot\\mathbf\{e^\{\-\\theta^\{2\}/2\}\}I\.8\.14\(x2−x1\)2\+\(y2−y1\)2\\sqrt\{\(x\_\{2\}\-x\_\{1\}\)^\{2\}\+\(y\_\{2\}\-y\_\{1\}\)^\{2\}\}\(𝐱𝟏−𝐱𝟐\)𝟐\+\(𝐲𝟏−𝐲𝟐\)𝟐\\mathbf\{\\sqrt\{\(x\_\{1\}\-x\_\{2\}\)^\{2\}\+\(y\_\{1\}\-y\_\{2\}\)^\{2\}\}\}LLM\-SRBenchII\.6\.15b8​π​Ef​ϵ​r33​sin⁡\(2​θ\)\\displaystyle\\frac\{8\\pi E\_\{f\}\\epsilon r^\{3\}\}\{3\\sin\(2\\theta\)\}4​π3⋅𝐄𝐟​ϵ​𝐫𝟑sin⁡θ​cos⁡θ\\displaystyle\\frac\{4\\pi\}\{3\}\\cdot\\frac\{\\mathbf\{E\_\{f\}\\epsilon r^\{3\}\}\}\{\\mathbf\{\\sin\\theta\\cos\\theta\}\}I\.48\.2−c​1−c4​m2/En2\-c\\sqrt\{1\-c^\{4\}m^\{2\}/E\_\{n\}^\{2\}\}−𝐜​1−𝐦𝟐​𝐜𝟒/𝐄𝐧𝟐\-\\mathbf\{c\}\\sqrt\{1\-\\mathbf\{m^\{2\}c^\{4\}/E\_\{n\}^\{2\}\}\}Real\-WorldOscillator 10\.8​sin⁡x−0\.5​v30\.8\\sin x\-0\.5v^\{3\}−0\.2​x3−0\.5​x​v−x​cos⁡x\-0\.2x^\{3\}\-0\.5xv\-x\\cos x−0\.5​\(𝐯𝟑\+𝐱𝐯\)\-0\.5\(\\mathbf\{v^\{3\}\}\+\\mathbf\{xv\}\)\+sin⁡𝐱⁡\(−0\.275​cos⁡𝐱\+0\.075\)\+\\mathbf\{\\sin x\}\(\-0\.275\\mathbf\{\\cos x\}\+0\.075\)Oscillator 2