Source-linked AI summary

Learning to Design SR Solvers: Harness Optimization for Symbolic Regression

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

arXiv:2610.16966v1cs.LG

TL;DR

SR agents increasingly use LLMs and tools, but their workflows and prompts are commonly designed by hand even though harness choices shape equation search and feedback use. HarnessSR optimizes role prompts and workflows, then refines instructions using SR-specific feedback; on LSR-Synth it achieves a mean ID Acc0.01 of 61.33% across three backbones, 45.0% above the next-best baseline.

  • Problem

    SR solver harnesses are largely designed by hand, despite harness choices shaping how equation hypotheses are proposed and how evidence guides subsequent search.

    Source: §1 · pp. 1–2
  • Method

    HarnessSR jointly optimizes role instructions, demonstrations, and workflow organization in three stages, then uses SR execution feedback to refine role instructions while keeping model weights and underlying search algorithms fixed.

    Source: §1 · p. 2
  • Results

    Across three backbones on LSR-Synth, HarnessSR achieves a mean ID Acc0.01 of 61.33%, a 45.0% relative improvement over the next-best baseline.

    Source: §1 · p. 2
  • Takeaways & Limitations

    Results on LSR-Synth support harness optimization as a way to construct stronger, reusable SR agents and to use numerical feedback in developing search strategies reused across tasks.

    Source: §5 · p. 9

    The present formulation searches within a predefined workflow space, leaving open whether more open-ended structures could reveal additional reusable strategies.

    Source: §5 · p. 10

Abstract

from arXiv · show

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.

1 INTRODUCTION

Symbolic regression turns observations into explicit equations that can serve as testable scientific hypotheses. HarnessSR treats the agent harness as an optimization target, jointly optimizing prompts and workflows before refining instructions with SR-specific feedback.

1 INTRODUCTION

  • Symbolic regression produces inspectable equations that can be related to conservation laws, dimensional constraints, and domain theories, making them testable scientific hypotheses beyond numerical prediction.The paper situates this use across physics, biology, materials science, and other domains.Source: §1 · p. 1
  • Treating the harness as a design variable matters because its prompts and workflow shape search priors, information flow, feedback use, and allocation of finite computation.The authors argue that these choices affect the quality and efficiency of equation discovery.Source: §1 · p. 2
  • Across three backbones on LSR-Synth, HarnessSR achieves a mean ID Acc0.01 of 61.33%, a 45.0% relative improvement over the next-best baseline.The framework jointly optimizes prompts and workflows and refines role instructions using SR execution feedback.Source: §1 · p. 2
  • HarnessSR optimizes role instructions, demonstrations, and workflow organization through three stages: role-prompt optimization, joint workflow–prompt search, and SR-feedback-guided instruction refinement.The resulting top-performing harness is frozen and reused on new SR tasks while model weights and underlying search algorithms remain fixed.Source: §1 · p. 2

2 RELATED WORK

Prior SR methods range from evolutionary and neural-guided search to LLM-based equation generation and tool-using agents. The paper identifies a remaining focus: optimizing the overarching solver harness itself, rather than only adapting search within a largely predefined harness.

Related SR solvers

  • Classical SR includes genetic programming, while PySR combines evolutionary search with expression simplification and constant optimization; neural-guided methods learn distributions over expressions.The passage names DSR as a recurrent-neural-network example guided by risk-seeking policy gradients.Source: §2 · p. 2
  • LLM-based SR methods generate structural hypotheses and combine them with numerical optimization, residual feedback, concept libraries, adaptive operators, or reflective memory.The cited examples include LLM-SR, SGA, LaSR, RESTART, and Deliberate Evolution.Source: §2 · p. 2
  • Agentic SR systems use tools and feedback-driven interaction, including autonomous analysis and evaluation, physical-property extraction, and hierarchical coordination.The examples are SR-Scientist, KeplerAgent, and A-SR.Source: §2 · p. 2
  • Representative SR systems retain largely hand-designed overarching harnesses, even when agents autonomously select tools within predefined interfaces, templates, and control rules.The paper therefore focuses on optimizing the solver harness under a limited budget to obtain reusable SR solvers.Source: §2 · p. 3

Self-Evolving Agents

  • Methods for self-evolving agents optimize instructions, demonstrations, component composition, or executable code, while the cited work notes limited attention to harness optimization for scientific discovery.The passage describes MIPRO, MASS, and Meta-Harness.Source: §2 · p. 3

3 METHOD

HarnessSR optimizes role prompts and workflow structure while holding the base language model and underlying search algorithms fixed. Its three-stage procedure moves from role-level prompt initialization to joint workflow–prompt search and SR-guided instruction refinement.

3.1 PROBLEM FORMULATION

  • The harness comprises role prompts, workflow structure, and shared fixed tools, and is optimized on training tasks for average equation quality under workflow-validity and execution-budget constraints.Role prompts include instructions and demonstrations; workflows control role activation, parallel search, and iterative refinement.Source: §3.1 · pp. 3–4

3.2 EXECUTION WITHIN AN SR AGENT HARNESS

  • At execution, an optional Analyst examines the data, branches run Planner-generated searches with Diagnoser–Repairer rounds, and optional cross-branch replanning and synthesis add further searches.Tool-specific adapters translate role outputs into PySR search parameters or LLM-generator prompts and parameterized program skeletons, while leaving the base algorithms unchanged.Source: §3.2 · pp. 4–5

3.3 STAGE 1: ROLE-LEVEL PROMPT OPTIMIZATION

  • Stage 1 optimizes each role’s instructions and demonstrations within simple workflows, retaining three high-scoring prompt packages per role for later search.Planner prompts are optimized first; the other roles are optimized with the Planner fixed, and Diagnoser–Repairer prompts are retained as paired configurations.Source: §3.3 · p. 5

3.4 STAGE 2: JOINT WORKFLOW AND PROMPT OPTIMIZATION

  • Stage 2 jointly searches valid workflow configurations and assignments from the Stage 1 prompt libraries, using TPE to propose candidates based on evaluation scores.Promising candidates then receive brief instruction adaptation with workflows and demonstrations fixed, addressing prompt–workflow compatibility before the larger Stage 3 budget.Source: §3.4 · pp. 5–6

3.5 STAGE 3: SR-GUIDED INSTRUCTION REFINEMENT

  • Stage 3 fixes each candidate’s workflow, demonstrations, and shared components, and refines only role instructions using SR execution traces and diagnostic feedback.Diagnostics assess evidence-to-action alignment, structural exploration, constraint justification, and expression quality; revised child harnesses are evaluated against parents, and the top harness is selected across workflows.Source: §3.5 · pp. 6–7

4 EXPERIMENTS

HarnessSR is evaluated on LSR-Synth across three backbones, with comparisons covering accuracy, solving time, model-weight post-training, call efficiency, optimization stages, and workflow search. Its strongest reported results are higher accuracy than the next-best baseline and gains from combining workflow search with SR diagnostics.

4.2 MAIN RESULTS

  • 45.0% and 63.8% higher mean ID Acc0.01 and Acc0.001, respectively, than SR-Scientist; mean solving time is reduced by 58.0–85.7% versus the fastest LLM-based baseline in each model group.HarnessSR leads at both accuracy thresholds across all four domains with the two Qwen models, and in physics, biology, and materials science with GPT-OSS-20B.Source: Table 3 · p. 7, §4.2 · p. 8

4.1 EXPERIMENTAL SETUP

  • The evaluation uses 129 LSR-Synth tasks across four scientific domains, three backbones, and non-LLM and LLM-based baselines.Harness optimization uses 1,024 public synthetic training tasks; accuracy-to-tolerance and mean solving time are reported, while optimization effort is measured by solver executions.Source: §4.1 · pp. 7–8

4.3 COMPARISON WITH REINFORCEMENT-LEARNING POST-TRAINING

  • HarnessSR reaches 57.36% ID Acc0.01 versus 40.92% for SR-Scientist+RL on the same Qwen3-Coder base model, and also leads at ID Acc0.001.HarnessSR keeps model weights fixed while optimizing role search instructions and their interactions to improve guidance for SR tools.Source: Table 4 · p. 8, §4.3 · p. 8

4.4 LLM CALLS VERSUS ACCURACY

  • HarnessSR achieves the highest accuracy at both thresholds with an average of 23 LLM calls per task, about 6.2% of SR-Scientist’s final-point call count.Baseline call checkpoints are compared with HarnessSR’s full execution of its selected solver configuration; competing methods remain less accurate despite hundreds or up to a thousand calls.Source: Figure 3 · p. 8, §4.4 · pp. 8–9

4.5 CONTRIBUTION OF OPTIMIZATION STAGES

  • ID Acc0.001 rises from 16.54% to 32.82% through the first two optimization stages, which optimize prompts in basic workflows and then jointly search workflows and prompt assignments.Stage 2 briefly adapts role instructions to fit the selected workflows’ role combinations and interactions.Source: Table 5 · p. 9, §4.5 · p. 9
  • Stage 3 raises ID Acc0.01 from 43.41% to 65.37% and ID Acc0.001 from 32.82% to 38.76% through reflection on SR diagnostics while keeping the workflow fixed.The diagnostics connect candidate-equation structures and residual patterns with role actions to guide subsequent equation search.Source: §4.5 · p. 9

4.6 WORKFLOW SEARCH AND DIAGNOSTIC FEEDBACK

  • Combining full-space search with SR diagnostics raises ID Acc0.01 from 49.35% to 65.37% and reduces mean LLM calls per task from 37 to 23, compared with fixed-pool search using execution traces.SR diagnostics improve accuracy under both search strategies, with larger gains in full-space search.Source: Table 6 · p. 9, §4.6 · p. 9

5 CONCLUSION AND LIMITATION

HarnessSR treats the SR agent harness as an optimization target and jointly optimizes role instructions, demonstrations, and workflow organization while keeping model weights and search algorithms fixed. The authors identify the predefined workflow space as a boundary of the present formulation.

5 CONCLUSION AND LIMITATION

  • HarnessSR jointly optimizes role instructions, demonstrations, and workflow organization while keeping model weights and underlying search algorithms fixed.The conclusion presents reusable SR agents as the supported outcome of this approach on LSR-Synth.Source: §5 · p. 9
  • Solver organization is proposed as an optimization dimension alongside equation representations and search algorithms.The authors connect numerical feedback to both within-task equation refinement and search strategies reused across tasks.Source: §5 · pp. 9–10
  • The present formulation searches within a predefined workflow space; whether more open-ended structures reveal additional reusable strategies remains open.Source: §5 · p. 10

A ONE-TIME OPTIMIZATION COST

The section compares one-time optimization investment for fixed-weight harness search and reinforcement-learning adaptation, then describes the symbolic-accuracy evaluation.

A ONE-TIME OPTIMIZATION COST

  • HarnessSR keeps model weights fixed while searching reusable role prompts and workflows; SR-Scientist+RL updates weights through reinforcement learning.Table 7 reports the one-time investment for obtaining each final solver.Source: A ONE-TIME OPTIMIZATION COST · p. 13, Table 7 · p. 13
  • The comparison uses the same Qwen3-Coder-30B-A3B-Instruct base model for both adaptation routes.SR-Scientist+RL training task counts and solver execution counts refer to optimization-stage activity.Source: A ONE-TIME OPTIMIZATION COST · p. 13

B SYMBOLIC ACCURACY

C.1 TRAINING TASKS AND SCIENTIFIC DOMAINS

The setup uses LSR-Synth tasks for harness optimization and evaluates the frozen solver on official tasks spanning four scientific domains.

C.1 TRAINING TASKS AND SCIENTIFIC DOMAINS

  • LSR-Synth combines established scientific-equation terms with newly constructed terms; the study uses it because transformed equations may retain pretraining clues.LSR-Transform instead modifies existing equations.Source: §C.1 · p. 13
  • Harness optimization draws from 1,024 public synthetic training tasks, with the training-scale comparison using 80, 160, 320, and 640 tasks.Training tasks provide role demonstrations, assess candidate prompts and workflows, and supply feedback for instruction revision.Source: §C.1 · p. 13
  • The final solver harness is frozen before evaluation on the 129 official LSR-Synth tasks.Each task provides observed data and separate in-distribution and out-of-distribution test points.Source: §C.1 · p. 13
  • LSR-Synth spans physics, chemistry, biology, and materials science, with Table 9 summarizing representative tasks and evaluation counts.The domain examples follow LLM-SRBench (Shojaee et al., 2025b).Source: §C.1 · p. 13, Table 9 · p. 14

C.2 LOCKED-EQUATION EVALUATION AND AGGREGATION

Evaluation locks one equation per task, tests it on separate ID and OOD points, and aggregates accuracy using a trimmed relative-error criterion.

C.2 LOCKED-EQUATION EVALUATION AND AGGREGATION

  • Each method locks one equation per task, and an external evaluator computes predictions on separate ID and OOD test points.The solver receives observed data, variable descriptions, and scientific context.Source: §C.2 · p. 14
  • The evaluation discards the largest 5% of sorted relative errors and tests task-level accuracy at thresholds 0.01 and 0.001.The trimming count is floor(0.05n_q) for a task with n_q test points.Source: §C.2 · p. 14
  • Accuracy in percent is aggregated over a set of N tasks using the remaining task-level indicators.The supplied definition introduces the aggregation formula but does not provide its displayed expression.Source: §C.2 · p. 14
  • ID and OOD evaluation use the same locked equation and numerical criterion, but draw test points from different input distributions..Source: §C.2 · p. 14

C.3 BASE MODELS AND PER-TASK TIME

The evaluation fixes model weights during harness optimization, defines per-task timing, and reports resource controls and optimization settings.

C.3 BASE MODELS AND PER-TASK TIME

  • Harness optimization modifies role prompts and workflows while holding base-model weights fixed, with at most 40 role-model calls and 24 SR backend calls per task.The three backbones are Qwen3.6-35B-A3B, Qwen3-Coder-30B-A3B-Instruct, and GPT-OSS-20B, serving different comparison roles.Source: §C.3 · p. 14
  • Mean per-task time runs from task start to locking the final equation and includes model calls and SR backend solving.It is averaged across evaluation tasks; one-time harness optimization cost is reported separately.Source: §C.3 · p. 14
  • Timing runs methods sequentially with exclusive use of shared GPU inference resources, while each method retains its own parallelism.Resource monitoring found no sustained queues, GPU request counts below limits, KV-cache occupancy at or below 10%, and normalized CPU utilization at or below 60%.Source: §C.3 · p. 15

C.4 HARNESS OPTIMIZATION SETTINGS

  • The 160-task main optimization set is split into three 32-task stage sets and a disjoint 64-task evaluation set for candidate selection.Each stage uses its corresponding D set for optimization, while Deva is used for selection.Source: §C.4 · p. 15
  • Stage 2 initializes 8 workflows and 16 prompt-assigned harnesses, then evaluates 8 more configurations for 24 total before brief adaptation of the top 6.Brief adaptation uses an eight-task subset of D2.Source: §C.4 · p. 15
  • Stage 3 performs at most four instruction-refinement rounds for each of the six selected candidates.Prompt-bundle construction, candidate acceptance, and final selection follow Section 3 procedures.Source: §C.4 · p. 15

D COMPARISON METHODS

The comparison spans evolutionary, neural, LLM-assisted, and tool-using SR solvers, while backend experiments test HarnessSR with PySR, LLM equation generation, and a hybrid. Call-accuracy trajectories track how selected equations perform as inference calls or iterations accumulate.

D COMPARISON METHODS

  • The evaluated baselines cover evolutionary, neural-symbolic, LLM-assisted, and autonomous tool-interaction approaches.The set includes GPLearn, E2E, NeSymReS, DSR, uDSR, PySR, LaSR, LLM-SR, SR-Scientist, and SR-Scientist+RL.Source: D COMPARISON METHODS · pp. 15–16
  • The comparison includes SR-Scientist+RL to contrast model-weight post-training with harness optimization in final solving performance.The variant retains SR-Scientist’s inference-time tool workflow.Source: D COMPARISON METHODS · p. 16

D.1 CALL-ACCURACY TRAJECTORIES

  • Call-accuracy comparisons use actual cumulative LLM calls for LLM-SR and SR-Scientist iterations, with accuracy evaluated across the same 129 tasks.For SR-Scientist, tasks ending early retain their final equation and cumulative calls at later checkpoints.Source: §D.1 · pp. 16–17

E SYMBOLIC REGRESSION BACKENDS AND CANDIDATE SELECTION

  • The study independently optimizes prompts and workflows for PySR, LLM equation generation, and a hybrid, then freezes each harness for evaluation.PySR role guidance maps to search settings; the LLM adapter inserts structural hypotheses into prompts and initializes parameterized program skeletons.Source: E SYMBOLIC REGRESSION BACKENDS AND… · p. 18, §E.1 · p. 18, §E.2 · p. 18

E.1 PYSR BACKEND

  • The PySR selector retains candidates close to the lowest observed-point error and chooses the expression with the fewest nodes.Both HarnessSR and the independent PySR baseline use this selection rule.Source: §E.1 · p. 18

E.3 HYBRID BACKEND

  • The hybrid completes separate PySR and LLM solves, pools their candidates, and selects one equation without further LLM calls or constant fitting.LLM-only and hybrid candidates are evaluated using the 95th percentile of observed-point relative error, expression complexity, and numerical validity.Source: §E.3 · p. 18
  • HarnessSR’s hybrid backend achieves higher accuracy than LLM-only and slightly exceeds the main-experiment PySR configuration.Prompts and workflows are independently optimized for each backend.Source: §E.3 · p. 19, Table 13 · p. 19

F OUT-OF-DISTRIBUTION RESULTS

HarnessSR reports OOD accuracy on LSR-Synth and illustrates its final PySR workflow and three successful equation-recovery cases. The OOD results compare accuracy across models and domains, while the cases show recovered structures and ID/OOD predictions.

F OUT-OF-DISTRIBUTION RESULTS

G.1 INITIAL INSTRUCTIONS FOR FIVE ROLES

  • The Analyst organizes context and observations, the Planner proposes a search plan, and the Diagnoser uses residuals to choose exploration or repair.The Repairer revises the plan, while the Synthesizer combines candidates across branches; the Analyst does not see ground truth or held-out rows.Source: §G.1 · p. 21

G.2 FINAL PYSR SOLVER WORKFLOW AND ROLE INSTRUCTIONS

  • The final PySR workflow runs three independent search branches, each with an initial plan and three diagnose–repair rounds, then synthesizes candidates for a final search.All candidates from initial search, repair, and synthesis are selected together; full execution uses 23 role-model calls and 13 SR-backend calls.Source: §G.2 · p. 21, Figure 7 · p. 22

H.1 POPULATION GROWTH: BPG10

  • The population-growth prediction retains fractional-power and linear terms with a shared coefficient, and passes both accuracy thresholds on ID and OOD test points.Its ID and OOD NMSE values are 2.9237 × 10^-15 and 4.5305 × 10^-14, respectively.Source: §H.1 · p. 26

H.2 NONLINEAR OSCILLATOR: PO18

  • For the nonlinear oscillator, the prediction retains the ground truth’s state-variable combinations and passes both thresholds, with ID and OOD NMSE below 2 × 10^-15.The terms couple a linear displacement-restoring component with velocity multiplied by sin(x).Source: §H.2 · p. 26

H.3 MATERIAL STRESS: MATSCI22

  • For material stress, a cosine term expresses the same periodic family as the ground-truth sine term with a phase shift, and both accuracy thresholds are passed.ID and OOD NMSE are 1.2594 × 10^-11 and 2.4181 × 10^-10, respectively.Source: §H.3 · pp. 26–27

I EFFECT OF TRAINING-SET SIZE

Larger harness-optimization training sets improve accuracy, with diminishing gains per doubling and greater one-time search effort at the largest scale.

I EFFECT OF TRAINING-SET SIZE

  • 52.71% to 65.37% Acc0.01 is the increase from 80 to 160 training tasks, with accuracy continuing to improve at 320 and 640 tasks but smaller gains per doubling.Table 16 reports ID accuracy and cumulative solver executions across training-set sizes; 160 tasks is the main setting.Source: I EFFECT OF TRAINING-SET SIZE · p. 27, Table 16 · p. 27
  • About 3.45 times as many solver executions raise Acc0.01 from 65.37% at 160 tasks to 69.51% at 640 tasks.The 640-task setting also improves accuracy at the stricter threshold.Source: I EFFECT OF TRAINING-SET SIZE · p. 27, Table 16 · p. 27
Loading 2610.16966v1…