Learning to Design SR Solvers: Harness Optimization for Symbolic Regression

Qian Li; Yuxiao Hu; Qihao Duan; Jinyue Yan; Yuntian Chen

cs.LGOct 12, 2026

Learning to Design SR Solvers: Harness Optimization for Symbolic Regression

Authors: Qian Li, Yuxiao Hu, Qihao Duan, Jinyue Yan, Yuntian Chen

Organizations: Shanghai Jiao Tong University, Shanghai, China · Eastern Institute of Technology, Ningbo, China · The Hong Kong Polytechnic University, Hong Kong, China · Epigenetics Laboratory, Max Planck Institute for Heart and Lung Research

Abstract

Symbolic regression (SR) seeks explicit equations that explain observed data, and recent methods increasingly employ large language models (LLMs) and agents for equation discovery. However, these systems typically rely on manually designed workflows and prompts. We argue that the SR agent harness should itself be an optimization target: its design shapes search priors and how numerical feedback guides subsequent equation refinement. We therefore propose HarnessSR, a framework for jointly optimizing prompts and workflows to construct reusable SR agents. During both harness optimization and agent execution, HarnessSR supports either PySR or an LLM-based equation generator through tool-specific interfaces, allowing the framework to optimize search guidance for different equation-discovery mechanisms. HarnessSR first optimizes role prompts through equation-discovery feedback in basic workflows, then jointly searches prompt assignments and workflows with brief adaptation to align the two. Finally, it refines role instructions through reflection on residual patterns, candidate equation structures, and search outcomes, improving how these signals guide subsequent equation search. On LSR-Synth in LLM-SRBench, HarnessSR achieves a mean ID $\mathrm{Acc}_{0.01}$ of 61.33% across three backbones, a 45.0% relative improvement over the next-best baseline.

Figures & tables

Figure 1: Overview of HarnessSR. It initializes role prompts, jointly optimizes workflows and prompt assignments, and refines instructions with SR-specific feedback. The top-performing agent harness is then frozen and reused on new tasks. The agent workflow shown below is illustrative.

Parameter	Values	Effect on the workflow
$u_{\mathrm{ana}}$	$\{0,1\}$	Activates Analyst to examines observational evidence before initial planning.
$b$	$\{1,2,3,4,5\}$	Number of parallel search branches, each executing Planner $\rightarrow$ PySR, or Planner $\rightarrow$ $\cdots$ $\rightarrow$ PySR
$t_{\mathrm{rep}}$	$\{0,1,2,3\}$	Number of search-plan refinement rounds in each branch, following Diagnoser $\rightarrow$ Repairer $\rightarrow$ PySR.
$t_{\mathrm{dis}}$	$\{0,1,2\}$	Number of cross-branch discussion rounds; in each round, each branch uses the shared results to update its own search plan and run PySR.
$u_{\mathrm{syn}}$	$\{0,1\}$	Activates Synthesizer to integrate branch findings into a unified plan, followed by one additional PySR run.

Table 1: Effects of workflow parameters. PySR as an example of base SR tools.

Target	Initialization workflow	$w$
Planner	Planner $\to$ PySR $\to$ Selection	(0,1,0,0,0)
Analyst	Analyst $\to$ Planner $\to$ PySR $\to$ Selection	(1,1,0,0,0)
Diag–Re	Planner $\to$ PySR $\to$ Diagnoser $\to$ Repairer $\to$ PySR $\to$ Selection	(0,1,1,0,0)
Synthesizer	Two parallel [Planner $\to$ PySR] branches $\to$ Synthesizer $\to$ PySR $\to$ Selection	(0,2,0,0,1)

Table 2: Local workflows for prompt initialization. PySR as an example of base SR tools.

Figure 2: 2-branch search with 1 cross-branch discussion round:

w

=(0,2,0,1,0). PySR as an example.

Non-LLM methods
Method	Overall		Physics		Chemistry		Biology		Materials		Time $\downarrow$ [-1pt]min/task
Method	$\mathrm{Acc}_{0.01}$	$\mathrm{Acc}_{0.001}$	$\mathrm{Acc}_{0.01}$	$\mathrm{Acc}_{0.001}$	$\mathrm{Acc}_{0.01}$	$\mathrm{Acc}_{0.001}$	$\mathrm{Acc}_{0.01}$	$\mathrm{Acc}_{0.001}$	$\mathrm{Acc}_{0.01}$	$\mathrm{Acc}_{0.001}$	Time $\downarrow$ [-1pt]min/task
GPLearn	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	3.15
E2E	0.78	0.00	0.00	0.00	0.00	0.00	0.00	0.00	4.00	0.00	0.43
NeSymReS	4.65	1.29	0.00	0.00	4.63	0.00	5.56	0.00	12.00	6.67	15.84
DSR	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	1.66

Table 3: ID accuracy (%) and mean solving time (min/task) on LSR-Synth. Bold indicates the best result in each column.

Method	Optimized component	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$
SR-Scientist	None	32.30	16.02
SR-Scientist+RL	Model weights	40.92	20.69
HarnessSR	Prompts and workflow	57.36	37.73

Table 4: Two optimization routes based on Qwen3-Coder-30B-A3B-Instruct and their ID accuracy (%).

Figure 3: Mean LLM calls per task versus ID accuracy on LSR-Synth (logarithmic horizontal axis). Baseline curves show call checkpoints; HarnessSR shows the full execution of its selected solver configuration. Curve construction is detailed in Appendix D.1 .

Stage	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$	Solver executions
Initial harness	5.43	1.55	0
After Stage 1	27.91	16.54	1,792
After Stage 2	43.41	32.82	2,752
After Stage 3 (full)	65.37	38.76	4,728

Table 5: ID accuracy (%) and cumulative solver executions after each optimization stage (Qwen3.6-35B-A3B, LSR-Synth). The initial harness uses initial prompts and a single-branch workflow.

Workflow search	Refinement feedback	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$	LLM calls/task $\downarrow$
Fixed-pool search	Execution traces	49.35	28.17	37
Full-space search	Execution traces	53.49	33.07	22
Fixed-pool search	SR diagnostics	56.59	32.30	37
Full-space search	SR diagnostics	65.37	38.76	23

Table 6: Workflow search and feedback on LSR-Synth (Qwen3.6-35B-A3B): ID accuracy (%) and mean LLM calls per task. Fixed-pool search selects presampled workflows; full-space search can propose new ones. Instructions are revised using execution traces directly or SR diagnostics of these traces.

Appendix figures & tables14 assets

Supplementary material from the paper’s appendix.

Appendix

Method	Training tasks	Solver executions	Optimization GPUs	Final product
SR-Scientist+RL	1,024	15,360	$32\times$ H200	Updated model weights
HarnessSR	160	4,728	$2\times$ H100	Role prompts and solver workflow

Appendix

Table 7: One-time optimization investment for two adaptation routes.

Method	SA $\uparrow$
uDSR	0.52
PySR	5.17

Appendix

Table 8: Symbolic accuracy (SA, %) on LSR-Synth.

Domain	Typical problem	Observed variables and equation target	Tasks
Physics	Damped oscillators and other dynamics	Time $t$ , position $x(t)$ , velocity $v(t)$ ; acceleration relation	44
Chemistry	Reaction kinetics	Time $t$ , reactant concentration $c(t)$ ; rate of concentration change	36
Biology	Population growth and dynamics	Time $t$ , population $P(t)$ ; population growth rate	24
Materials science	Stress–strain relations	Strain $\epsilon$ and temperature $T$ ; stress $\sigma$	25
Total			129

Appendix

Table 9: Four scientific domains and evaluation task counts in LSR-Synth.

Cumulative calls	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$
4	2.33	0.00
8	2.33	0.00
12	2.33	0.00
16	2.33	0.00
20	2.33	0.00
40	6.20	0.00

Appendix

Table 10: LLM-SR call checkpoints and ID accuracy over 129 tasks (%).

Iteration	Mean cumulative calls	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$
1	12.40	9.30	5.43
5	54.69	21.71	11.63
10	103.37	29.46	13.95
15	149.98	33.33	17.83
20	196.41	37.21	19.38
25	240.91	42.64	21.71

Appendix

Table 11: SR-Scientist iteration checkpoints and ID accuracy over 129 tasks (%).

Iteration	Mean cumulative calls	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$
1	40.63	0.78	0.00
5	135.33	4.65	2.33
10	253.52	6.20	2.33
15	371.05	5.43	3.10
20	499.16	8.53	3.88
25	627.83	10.08	4.65

Appendix

Table 12: LaSR iteration checkpoints, actual mean cumulative calls, and ID accuracy over 129 tasks (%).

SR backend	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$
LLM	59.43	34.11
LLM + PySR	66.67	40.05

Appendix

Table 13: ID accuracy (%) on LSR-Synth with LLM-based and hybrid SR backends for Qwen3.6-35B-A3B. Role prompts and workflows are independently optimized for each backend.

Figure 4: ID and OOD

\mathrm{Acc}_{0.01}

(%) for methods using Qwen3.6-35B-A3B.

Figure 5: ID and OOD

\mathrm{Acc}_{0.01}

(%) for methods using Qwen3-Coder-30B-A3B-Instruct.

Figure 6: ID and OOD

\mathrm{Acc}_{0.01}

(%) for methods using GPT-OSS-20B.

Non-LLM methods
Method	Overall	Physics	Chemistry	Biology	Materials
Method	$\mathrm{Acc}_{0.01}$ $\uparrow$	$\mathrm{Acc}_{0.01}$ $\uparrow$	$\mathrm{Acc}_{0.01}$ $\uparrow$	$\mathrm{Acc}_{0.01}$ $\uparrow$	$\mathrm{Acc}_{0.01}$ $\uparrow$
PySR	31.78	27.27	18.52	11.11	78.67
Qwen3.6-35B-A3B
LLM-SR	34.88	29.55	25.00	12.50	80.00
SR-Scientist	44.19	31.82	25.00	50.00	88.00

Appendix

Table 14: OOD accuracy (%) on LSR-Synth.

Figure 7: Final solver workflow for the PySR backend in the main experiments. Three branches search and repair independently; after synthesis, the final equation is selected from all candidates.

Case (domain)	ID NMSE	OOD NMSE	$\mathrm{Acc}_{0.01}$ ID / OOD	$\mathrm{Acc}_{0.001}$ ID / OOD
BPG10 (population growth)	$2.9237\times 10^{-15}$	$4.5305\times 10^{-14}$	Pass/Pass	Pass/Pass
PO18 (physical oscillator)	$1.4546\times 10^{-15}$	$1.6621\times 10^{-15}$	Pass/Pass	Pass/Pass
MatSci22 (materials science)	$1.2594\times 10^{-11}$	$2.4181\times 10^{-10}$	Pass/Pass	Pass/Pass

Appendix

Table 15: Numerical results for three successful equation-recovery cases.

Training tasks	ID $\mathrm{Acc}_{0.01}$ $\uparrow$	ID $\mathrm{Acc}_{0.001}$ $\uparrow$	Solver executions
80	52.71	28.94	2,648
160	65.37	38.76	4,728
320	68.48	41.34	9,016
640	69.51	42.12	16,312

Appendix

Table 16: ID accuracy (%) and cumulative solver executions across training-set sizes (Qwen3.6-35B-A3B, LSR-Synth). Each nested training set undergoes all three optimization stages; 160 tasks correspond to the main setting.

Explore similar work

Aug 5, 2026

A-SR: Self-Evolving Agentic LLMs for Symbolic Regression via Hierarchical Coordination

Symbolic regression aims to discover closed-form equations from data, but existing LLM-guided methods often rely on a unified proposal loop that compresses heterogeneous search failures into a scalar score and a single prompt. We propose A-SR, a self-evolving agentic framework that shifts the control unit from expression edits to role-conditioned evidence views. A-SR coordinates formula discovery through routing among coordination protocols, an online evaluator-reward role policy, and state-routed process memory. During search, evaluator feedback characterizes reliability and productivity, updates role-level utilities, and routes elite motifs, failure traces, and validity diagnostics to different agents. The framework self-evolves at two timescales: within a run, it adapts the search process without updating LLM parameters; across runs, recorded trajectories can be distilled into open-source LLMs as role-conditioned proposal priors. Averaged over the four LSR-Synth scientific domains in LLM-SRBench, A-SR improves [email protected] over baselines from 25.79% to 48.30% with Llama3.1-8B, while A-SR-LoRA improves the corresponding Qwen3-4B result from 24.58% to 38.29%. On four real-world scientific discovery tasks, A-SR obtains the best in-distribution or out-of-distribution normalized mean squared error on 7 of 8 reported metrics.

Wenxiao Zhao, Dong Liu, Kaiyi Xu +9

Shanghai AI Laboratory · University of California at Los Angeles · Tongji University

Multi-Agent LLM Systems LLM Agent Training Large Language Model-Guided Scientific Discovery Symbolic Regression

Ask PDF

Sep 28, 2026

SRHarness: A Harness for Agentic Symbolic Regression

Recent agentic symbolic regression approaches increasingly rely on large language models to analyze data, select scientific operations, and refine hypotheses over long search trajectories. In such systems, performance depends not only on the underlying model and search strategy, but also on the runtime infrastructure that supports scientific search. We introduce SRHarness, a domain-specific harness for agentic symbolic regression built around three mechanisms: composable scientific actions that provide a common interface over raw, transformed, and candidate-derived quantities; persistent scientific state that retains evaluated hypotheses and exposes compact model-facing views; and trajectory lifecycle management that coordinates continuation, branching, restart, and termination. On LLM-SRBench, SRHarness consistently improves both numerical generalization and symbolic recovery under matched LLM backbones. With DeepSeek-v4-flash-0731, it achieves 93.69% symbolic accuracy on LSR-Transform, compared with 62.16% for SR-Scientist, and retains 72.97% accuracy on an anonymized variant that removes scientific descriptions and variable semantics, versus 39.64% for SR-Scientist. Under the same DeepSeek-v4-flash-0731 backbone, SRHarness also substantially outperforms Codex (72.97% vs. 20.72%) and reaches performance comparable to Codex with GPT-5.5, while simply providing Codex with the same scientific tools does not reproduce this advantage. These results show that effective agentic symbolic regression depends not only on models or tools, but also on structured runtime support for organizing scientific actions, accumulated hypotheses, and long-horizon search.

Zihan Yu, Shixuan Zhou, Hao Huang +2

Department of Electronic Engineering, BNRist, Tsinghua University, Beijing, China · College of AI, Tsinghua University, Beijing, China · Department of Earth System Science, Tsinghua University, Beijing, China

Long-Horizon LLM Agents Large Language Model-Guided Scientific Discovery LLM Agent Harnesses Symbolic Regression

Ask PDF

Aug 26, 2026

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

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.

Yating Ling, Wenjing Cun, Zhitang Chen

Symbolic Regression Large Language Model-Guided Scientific Discovery Genetic Programming

Ask PDF