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

A-SR: Self-Evolving Agentic LLMs for Symbolic Regression via Hierarchical Coordination
Aug 5, 2026Wenxiao Zhao, Dong Liu, Kaiyi Xu +9Multi-Agent LLM Systems LLM Agent Training
Shanghai AI Laboratory · University of California at Los Angeles · Tongji University
SRHarness: A Harness for Agentic Symbolic Regression
Sep 28, 2026Zihan Yu, Shixuan Zhou, Hao Huang +2Long-Horizon LLM Agents Large Language Model-Guided Scientific Discovery
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
InsightSR: Refining Symbolic Regression Search Spaces via Parallel Semantic and Structural LLM Guidance
Aug 26, 2026Yating Ling, Wenjing Cun, Zhitang ChenSymbolic Regression Large Language Model-Guided Scientific Discovery

Aug 5, 2026J/K move · Enter open · S save

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

Wenxiao Zhao, Dong Liu, Kaiyi Xu, Feng Liu+8

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

Open paper AskPDF arXiv

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.

Abstract

Explore similar work

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

SRHarness: A Harness for Agentic Symbolic Regression

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