From Score to Diagnosis: Search-Independent Post-Diagnosis for Symbolic Regression

Introduction

Symbolic regression (SR) discovers compact, interpretable mathematical expressions from observed data, enabling data-driven scientific discovery across physics, biology, chemistry, and engineering [44, 52, 18, 20, 43] (Udrescu & Tegmark 2020). Given observations $\{(x_i, y_i)\}_{i=1}^{N}$, SR searches for an expression $f : x \mapsto y$ that explains the underlying relationship. Modern SR solvers employ diverse search strategies to enable efficient and accurate SR: PySR uses multi-population evolutionary search [5]; LLM-SR iteratively refines equations using the code-generation and reasoning capabilities of LLMs [38]; LaSR co-evolves equations with a reusable concept library [10]; and SR-Scientist autonomously integrates data analysis, tool use, and equation refinement [45]. This diversity is valuable because different solvers offer complementary strengths.

Despite their diverse search strategies, most SR solvers rely on the same impoverished search feedback: a scalar fitness score, such as MSE or NMSE. By collapsing the residual landscape into a single value, these scores obscure error geometry and reveal little about where, why, or how a candidate fails. Indeed, similar scores may correspond to fundamentally distinct failure modes (Li et al.; Xiang et al.). This limited feedback leaves even sophisticated search mechanisms poorly informed, resulting in inefficient and misdirected exploration.

Recent methods RESTART and Deliberate Evolution (DE) design new SR solvers that embed structured diagnostics into their search loops, providing clearer directions for candidate refinement than scalar errors alone [32] (Li & Pan 2026). However, their feedback mechanism depend on solver-specific search information, such as search trajectories memory, parent-selection strategy, and operator policy. Thus, feedback mechanisms are tightly coupled to the method’s search process and remain an internal solver capability rather than reusable feedback. This coupling creates an undesirable trade-off: obtaining such enhanced feedback mechanism generally requires adopting the accompanying search framework, rather than retaining the distinctive mechanisms of existing solvers, such as LaSR’s concept–equation co-evolution, SR-Scientist’s autonomous tool-use loop, and PySR’s efficient genetic-programming engine.

This motivates a different question: can richer diagnostic feedback be made reusable across diverse solvers, so that their shared feedback bottleneck can be addressed without sacrificing their distinctive search mechanisms? Rather than introducing another new solver, we pursue a route orthogonal to solver design: constructing a diagnostic feedback mechanism that can be integrated with different SR solvers. Such a mechanism should be search-independent, neither depending on nor altering the solver’s underlying search mechanism; informative, revealing where and how a candidate fails; and solver-compatible, allowing diagnosis-guided revisions to be readily returned to the solver.

To instantiate this view, we propose PostDiag, a search-independent post-diagnostic module for SR. Unlike feedback mechanisms embedded in the search loop, PostDiag takes only candidate skeletons and observed data as input, without relying on solver-specific search mechanisms. It can therefore be integrated with solvers that employ different search mechanisms. For each candidate, PostDiag produces two complementary outputs: an evidence map that reveals where and how the candidate fails, and a risk localization that determines the appropriate intervention. PostDiag first introduces Patchwise Feature Profiling (PFP), which represents local candidate behavior using structured fingerprints of level, trend, and roughness. Rather than compressing diagnostic information into a single global error, discrepancies between these local fingerprints form an interpretable evidence map and support the estimation of Skeleton Expressiveness Risk and Skeleton Stability Risk. The former captures the structural mismatch remaining after coefficient fitting, while stability measures the candidate’s sensitivity to data resampling. The Diagnostic Risk Plane (DRP) uses these two quantities to assign the candidate to one of four intervention regions: preserve, repair, simplify, or rewrite. A conditional LLM then performs the corresponding skeleton intervention, and the revised candidate is returned to the original solver through a lightweight adapter. We integrate the PostDiag with different SR solvers spanning LLM-based and GP-based methods, and evaluate them on the LSR-Synth dataset of LLM-SRBench. The results show that PostDiag improves candidate refinement and accelerates convergence without replacing the solvers’ original search mechanisms. Ablation studies further validate the contribution of each diagnostic component.

Our contributions come in fourfold: (1) We introduce a solver-independent post-diagnostic perspective for SR that separates feedback design from search design, addressing a shared feedback bottleneck without replacing existing solvers. (2) We propose PostDiag, a reusable diagnostic module that takes only a candidate skeleton and observed data as input and returns evidence-supported diagnoses and candidate revisions through lightweight solver adapters. (3) We develop PFP to extract patch-wise diagnostic evidence and DRP to map skeleton adequacy and stability to targeted interventions. (4) We integrate the PostDiag with seven different LLM-based and GP-based SR solvers, demonstrating broad compatibility and consistent convergence improvements.

Related Work

Symbolic regression (SR) spans diverse search and learning paradigms, each contributing distinct inductive biases and capabilities [3, 36, 6, 16]. Genetic-programming methods evolve expression trees through population-based operators [19, 33, 31, 7, 40], with PySR providing efficient multi-population search, complexity control, and constant optimization [5]. Tree-search, reinforcement-learning, and neural methods instead leverage planning, learned generation policies, or pretrained representations [41, 34, 26, 30, 12, 2, 17, 37, 29 Li et al.]. More recently, LLM-based methods have incorporated scientific priors and reasoning: LLM-SR combines program proposal with evolutionary experience [38], LaSR co-evolves equations and a concept library [10], and many SR agent methods performs autonomous tool-assisted optimization [45, 15, 25]. Despite their diverse generation and search mechanisms, these methods share a common feedback bottleneck: candidate behavior is compressed into scalar fitness, which supports ranking but provides little information about failure modes, structural defects, or appropriate interventions.

Deliberate Evolution (DE) [32] augments an MSE-driven evolutionary process with data and residual diagnostics, adaptive search operators, and trajectory-level memory, demonstrating that structured feedback can improve symbolic search beyond reliance on scalar errors alone. However, DE addresses this shared feedback bottleneck by introducing a new end-to-end solver, in which diagnosis is tightly coupled with its parent-selection strategy, operator policy, and trajectory memory. Adopting DE therefore entails replacing the original solving framework, making it difficult to retain solver-specific strengths such as LaSR’s concept–equation co-evolution, SR-Scientist’s agentic closed-loop optimization, or PySR’s efficient genetic-programming engine. Its limitation lies in addressing a common challenge across SR methods through a dedicated solver rather than a capability that existing solvers can readily share.

PostDiag takes a different technical route. Instead of replacing existing solvers, it provides a shared diagnostic interface decoupled from their internal search processes. Given any candidate skeleton produced by the solver and the observation, PostDiag uses Patch-wise Feature Profiling to recover local behavioral information obscured by scalar errors and provide evidence for risk localization and candidate intervention. It can therefore enhance diagnostic feedback across heterogeneous SR solvers while preserving their original search mechanisms and distinctive strengths.

A Search-Independent Post-Diagnosis for SR

Problem Setup

Given a dataset $\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^{N}$, where $\mathbf{x}_i \in\mathbb{R}^{d}$ and $y_i \in\mathbb{R}$, we represent a candidate expression as $f_s(\mathbf{x}; \boldsymbol{\theta}_s)$. Here, $s$ denotes the symbolic skeleton and $\boldsymbol{\theta}_s$ is its numerical coefficients. PostDiag treats the skeleton $s$ as the object of diagnosis.

We randomly partition $\mathcal{D}$ into a disjoint fit set $\mathcal{D}_{\mathrm{fit}}$ and audit set $\mathcal{D}_{\mathrm{audit}}$. The fit set is used to estimate the numerical coefficients of the candidate skeleton, while the audit set is used to evaluate the fitted candidate and produce its diagnosis. This separation prevents coefficient fitting and candidate diagnosis from using the same observations.

Given a candidate skeleton $s$ and dataset $\mathcal{D}$, PostDiag outputs a structured diagnostic report $\mathcal{G}_s$ and a revised skeleton $\tilde{s}$:

Here, $\mathcal{G}_s=(c_s,\mathcal{E}_s)$ contains the candidate’s diagnostic coordinate $c_s$ on the diagnostic risk plane and supporting evidence $\mathcal{E}_s$, while $\tilde{s}$ is the candidate skeleton revised according to the diagnosis.

PostDiag neither generates initial candidates nor prescribes how a solver conducts its search. It diagnoses candidate skeletons supplied by the solver and returns the diagnostic results and revised skeletons through a unified interface. The same diagnostic core can therefore be integrated with different SR solvers through lightweight adapters.

Patch-wise Feature Profiling (PFP)

To support expressiveness and stability estimation without compressing evaluation into a scalar error, PostDiag uses Patch-wise Feature Profiling to construct structured fingerprints on the observation.

We first normalize each input dimension and denote the normalized coordinate of $x_i$ by $\tilde{x}_i$. We then construct $M$ local patches $\mathcal{P}=\{P_m\subseteq\mathcal{D}_{\mathrm{audit}}\}_{m=1}^{M}$. Such patches can be constructed in various ways. We adopt an anchor-$k$NN for its simplicity and robustness in multivariate settings. Specifically, we select spatially well-distributed anchors via farthest point sampling [9] and construct each patch from the $k$ nearest neighbors [4] of an anchor $\hat{x}_m$:

The patches may overlap, reducing sensitivity to neighborhood boundaries, while $M$ and $k$ are chosen to cover all audit samples.

Given a patch $P_m$, PFP summarizes the local behavior of $f$ using three complementary features.

Level. The level is the mean function value within the patch:

Equivalently, $\mu_m(f)$ is the best constant approximation of $P_m$ in the least-squares sense, i.e.,

Thus, $\mu_m(f)$ summarizes the zeroth-order behavior of $f$ on the patch $P_m$.

Trend. The trend captures the direction and magnitude in which $P_m$ varies across the normalized inputs. It is defined as the least-squares slope vector of $P_m$ on each direction:

where $\bar{x}_m=\frac{1}{k}\sum_{x_i\in P_m}x_i$ is the normalized patch centroid. $\beta_m(f)$ summarizes the first-order behavior of $f$ on the patch $P_m$. We retain $\beta_m(f)$ as a $d$-dimensional vector to capture trends in multivariate symbolic regression. Estimating local slopes via least squares, rather than endpoint differences, further improves robustness to endpoint anomalies and observation noise.

Roughness. Roughness captures the higher-order local behavior that cannot be explained by a constant level or a linear trend. For each $x_i\in P_m$, we can define the local affine approximation use Level and Trend as

We therefore denote the residual between $f$ and $\hat{f}$ on $P_m$ as an additional signal:

Accordingly, if $f$ is affine within $P_m$, its local behavior is fully explained by the Level and Trend terms, yielding $\rho_m(f)=0$. Local nonlinearity or curvature increases Roughness $\rho_m(f)$.

The local fingerprint of $f$ on $P_m$ is

and its complete fingerprint over the $\mathcal{D}_{\mathrm{audit}}$ is

By profiling different behaviors over the same patch partition, PFP provides a common basis for structured comparison. Each comparison is tied to a specific input-space patch and further decomposed into level, trend, and roughness features. It therefore reveals not only where behaviors differ, but also how they differ, while reducing reliance on individual pointwise residuals. In the next subsection, we use PFP as a comparison tool to assess skeleton adequacy and stability. Together, these comparisons yield the diagnostic coordinate $c_s$ and its supporting evidence map $\mathcal{E}_s$.

Diagnostic Risk Plane

PostDiag is designed to diagnose a candidate without depending on or altering the solver’s search mechanism. It therefore cannot rely on solver-specific signals such as search trajectories, visitation statistics, population states, or operator histories, and must infer the appropriate intervention from only the candidate skeleton and observed data. To meet this constraint, we introduce the Diagnostic Risk Plane (DRP), which characterizes each skeleton using two behavior-based risks: Skeleton Expressiveness Risk and Skeleton Stability Risk. The former asks whether the skeleton can represent the observed behavior after coefficient fitting, whereas the latter asks whether the skeleton exhibits consistent behavior when its coefficients are re-estimated from different subset of the fitting data.

Skeleton expressiveness risk. We quantify this risk by a fingerprint-based bias score $\operatorname{Bias}(s)$. Specifically, let $\tilde{\theta}_s$ denote the coefficients of skeleton $s$ fitted on $\mathcal{D}_{\mathrm{fit}}$, and let $\tilde{f}_s(x)=f_s(x;\tilde{\theta}_s)$. We evaluate $\tilde{f}_s$ on $\mathcal{D}_{\mathrm{audit}}$ and compare its fingerprint with that of the observed responses $\mathbf{y}_{\mathrm{audit}}$ over the same patches. For each patch $P_m$, we compute

The local expressiveness risk is

and the overall risk is

A low $\operatorname{Bias}(s)$ indicates that the fitted skeleton reproduces the observed behavior across patches, whereas a high value indicates insufficient expressiveness. The nonnegative weights $\omega_\mu$, $\omega_\beta$, and $\omega_\rho$ balance the three features. In practice, we assign a smaller weight to roughness because it is more sensitive to local fluctuations in finite patches.

Skeleton stability risk. Expressiveness alone cannot reveal whether the fitted behavior is reproducible or depends on the particular observations used for coefficient fitting. We characterize this dependence as skeleton stability risk. A stable skeleton should induce similar behavior when refitted on different resamples of the same observations, whereas a fragile skeleton relies on sample-specific coefficient configurations and therefore varies substantially across refits. We therefore quantify skeleton stability risk by the dispersion score $\operatorname{Var}(s)$ of the refitted function under different subset of the $\mathcal{D}_{\mathrm{fit}}$. Specifically, we generate $B$ fitting subsets $\{\mathcal{D}_b\}_{b=1}^{B}$ by bootstrap resampling [8] of the $\mathcal{D}_{\mathrm{fit}}$ and refit the coefficients of the same skeleton $s$ on each subset. This yields

. All refitted functions are evaluated on the same audit patches. For each patch $P_m$, we measure the dispersion of the fingerprints produced by the $\{f_s^{(b)}(x)\}_{b=1}^{B}$:

where $\overline{\beta}_m = \frac{1}{B} \sum_{b=1}^{B} \beta_m\left(f_s^{(b)}\right)$. The local stability risks of patch $P_m$ are

and the overall stability risks of skeleton $s$ can be defined as

A larger $\operatorname{Var}(s)$ indicates greater variation in the fitted behavior across resamples and thus higher skeleton stability risk.

Risk localization and diagnostic evidence. Each candidate $s$ is positioned in the Diagnostic Risk Plane (DRP) by the coordinate $c_s = \left(\operatorname{Bias}(f_s), \operatorname{Var}(f_s)\right)$. Meanwhile, the patch-level quantities are retained as the evidence map

Thus, $c_s$ provides a compact risk localization for selecting an intervention, while $\mathcal{E}_s$ identifies the patches and behavioral features responsible for that diagnosis.

Using thresholds $\tau_B$ and $\tau_V$, the DRP distinguishes four diagnostic states, each corresponding to a different failure mode and intervention strategy:

• SETTLED ($\operatorname{Bias} \leq\tau_B$, $\operatorname{Var} \leq\tau_V$): the skeleton is sufficiently expressive and stable. The intervention is therefore preservation-oriented: the skeleton is either preserved or only minimally refined according to the evidence map.

• UNDERFIT ($\operatorname{Bias} > \tau_B$, $\operatorname{Var} \leq\tau_V$): the skeleton exhibits stable but systematic behavioral discrepancies, indicating missing or incorrect local structure. The candidate should be repaired by adding or replacing local symbolic components guided by the evidence map.

• OVERFIT ($\operatorname{Bias} \leq\tau_B$, $\operatorname{Var} > \tau_V$): the skeleton matches the observations but is unstable under resampling, suggesting redundant or poorly identifiable structure. The candidate should be simplified.

• MISMATCHED ($\operatorname{Bias} > \tau_B$, $\operatorname{Var} > \tau_V$): the skeleton is neither expressive nor stable, indicating a more fundamental structural mismatch. The candidate should be rewritten.

The diagnostic coordinate determines what intervention is appropriate, while the evidence map indicates where and how the candidate fails. Both rely only on the candidate skeleton and observed data, enabling targeted diagnosis independent of the solver’s internal search process.

Diagnosis-Guided Intervention and Solver Integration

To reduce overhead, PostDiag first applies a low-cost NMSE-based gate to determine whether a candidate needs structured diagnosis. Candidates with either very low or very high errors are excluded from the PostDiag pipeline. Accurate candidates require no structural intervention, while poor candidates are left to the base solver’s native search mechanism for elimination or regeneration. PostDiag therefore focuses on intermediate candidates, for which more informative feedback is needed to determine an appropriate intervention. NMSE serves only as a coarse screening signal, whereas PFP and DRP provide the evidence and diagnostic decision that guide the intervention.

For each candidate entering the diagnosis pipeline, PostDiag produces $\mathcal{G}_s = (c_s, \mathcal{E}_s)$. The risk region determined by $c_s$ specifies the intervention type, while $\mathcal{E}_s$ identifies the patches and behavioral features that motivate the intervention. An LLM-based post-editor receives the original skeleton $s$, $c_s$, and $\mathcal{E}_s$, and proposes a diagnosis-guided revision $\tilde{s}$. The LLM acts only as a conditional skeleton editor and does not participate in the solver’s candidate generation or search control. Notably, a candidate entering the diagnostic process may still be classified as SETTLED. Although its scalar error is not yet low enough for direct acceptance, the DRP indicates that its skeleton is sufficiently expressive and stable. PostDiag still provides the LLM editor with the original skeleton, the preservation-oriented diagnosis, and the evidence map $\mathcal{E}_s$, but instructs it to either retain the skeleton or make only minor edits supported by the evidence.

Category	Method / setting	Overall Time (s)	Overall NMSE	Overall $A_{.01}$	Overall $A_{.001}$	Phys NMSE	Phys $A_{.01}$	Phys $A_{.001}$	MatSci NMSE	MatSci $A_{.01}$	MatSci $A_{.001}$	Chem NMSE	Chem $A_{.01}$	Chem $A_{.001}$	Bio NMSE	Bio $A_{.01}$	Bio $A_{.001}$
LLM-based	Direct LLM	30.5	7.72E-02	2.6	0.8	7.32E-02	3.0	0.8	6.24E-02	2.7	1.3	3.50E-02	2.8	0.9	1.63E-01	1.4	0.0
	Direct LLM+PostDiag	33.7	3.23E-02	4.4	2.3	5.22E-02	5.3	3.8	4.34E-03	4.0	2.7	8.76E-03	3.7	0.9	6.03E-02	4.2	1.4
	SGA	320.0	1.39E-02	11.9	5.9	9.43E-03	8.3	6.1	8.21E-03	12.0	5.3	4.92E-03	19.4	7.4	4.15E-02	6.9	4.2
	SGA+PostDiag	350.2	2.40E-03	18.6	9.8	2.66E-03	14.4	9.8	1.87E-04	25.3	6.7	8.65E-04	24.1	13.0	6.52E-03	11.1	8.3
	LLM-SR	167.5	1.45E-03	22.0	8.8	3.66E-03	14.4	8.3	3.12E-05	37.3	13.3	1.75E-04	27.8	8.3	8.15E-04	11.1	5.6
	LLM-SR+PostDiag	191.7	2.12E-04	31.0	15.2	3.48E-04	23.5	12.1	8.94E-06	60.0	30.7	5.84E-05	28.7	12.0	4.02E-04	18.1	9.7
	LaSR	559.2	1.06E-02	13.4	7.5	1.73E-02	11.4	6.8	2.48E-03	14.7	6.7	6.90E-03	18.5	11.1	1.22E-02	8.3	4.2
	LaSR+PostDiag	595.4	8.62E-04	21.4	12.7	1.52E-03	20.5	11.4	4.82E-05	29.3	14.7	6.45E-04	22.2	14.8	8.29E-04	13.9	9.7
	SR-Scientist	648.0	1.12E-03	25.6	11.9	2.74E-03	17.4	9.8	4.47E-05	46.7	21.3	1.49E-04	27.8	11.1	7.26E-04	15.3	6.9
	SR-Scientist+PostDiag	702.2	1.92E-04	35.7	20.4	3.11E-04	25.8	15.9	1.08E-05	62.7	37.3	6.47E-05	38.0	20.4	3.54E-04	22.2	11.1
GP-based	PySR	8.8	2.67E-03	23.8	11.4	4.11E-03	21.2	12.1	1.70E-04	44.0	18.7	8.93E-04	20.4	10.2	5.29E-03	12.5	4.2
	PySR+PostDiag	10.5	1.34E-03	28.7	13.2	1.82E-03	24.2	14.4	2.25E-04	50.7	17.3	5.29E-04	25.9	12.0	2.84E-03	18.1	8.3
	gplearn	95.5	1.34E-01	0.0	0.0	8.51E-02	0.0	0.0	2.36E-02	0.0	0.0	6.56E-02	0.0	0.0	4.42E-01	0.0	0.0
	gplearn+PostDiag	101.7	3.07E-02	1.0	0.0	4.29E-03	2.3	0.0	6.91E-03	0.0	0.0	4.29E-02	0.9	0.0	8.58E-02	0.0	0.0

We apply a lightweight validation before returning the proposed revision. Specifically, we refit the coefficients of $\tilde{s}$ on $\mathcal{D}_{\mathrm{fit}}$ and evaluate its prediction error on $\mathcal{D}_{\mathrm{audit}}$. The revision is accepted if it reduces the error; otherwise, it is rejected and $\tilde{s}=s$.

The revised skeleton can be integrated into different solvers via lightweight adapters. We consider two common solver families: GP-based and LLM-based solvers. For GP-based solvers, an accepted revision is inserted into the population as a repaired offspring. For LLM-based solvers, it replaces the original candidate, while $(s,\mathcal{G}_s,\tilde{s})$ is appended to the prompt as a diagnostic revision record. Additional adapter details are provided in Appendix F. PostDiag leaves candidate generation, selection, and exploration to the base solver, deriving diagnoses solely from the candidate skeleton and observed data without accessing solver-specific search states or histories.

Experiment

Experimental Setup

Dataset. Because PostDiag and most primary baselines use LLMs, we evaluate them on the LSR-Synth subset of LLM-SRBench (Shojaee et al., [39]). Unlike LSR-Transform, which may remain vulnerable to pretraining contamination (Xia et al., [45]), LSR-Synth combines canonical domain terms with newly synthesized terms and filters tasks for solvability, novelty, and scientific plausibility. It contains 129 tasks across physics (44), materials science (25), chemistry (36), and biology (24), each with in-domain and out-of-domain test sets.

Metrics. We evaluate the discovered equations from 3 perspectives: numerical accuracy, accuracy-to-tolerance, and search efficiency. Numerical accuracy is measured by the normalized mean squared error (NMSE) between the predictions of the discovered equation $f$ and the ground truth:

where $\hat{y}_i$ is the prediction target value and $y_i$ is the ground-truth value on $x_i$. NMSE normalizes the error by the target variance, making it directly comparable across tasks with different output scales. We further report $\mathrm{Acc}_{0.01}$ following Kamienny et al. [17]; Shojaee et al. [38]; Xia et al. [45] to show the task-level solving accuracy at stringency levels:

We report $\mathrm{Acc}_{\tau}$ for $\tau\in\{.01,.001\}$ Xia et al. [45]; tables use $A_\tau$. Search efficiency is measured by the average per-task search time (in seconds), which reflects the additional computational cost introduced by PostDiag.

Baselines. PostDiag is a search-independent diagnostic module rather than a standalone solver; we therefore compare each base solver with its PostDiag-enhanced version. We evaluate five LLM-based solvers: LLM-SR [38], LaSR [10], SR-Scientist [45], Direct LLM [39], and SGA [28]; and two GP-based solvers, PySR [5] and gplearn [40]. Integration details are in Appendix F.

Configuration of PostDiag. We describe the default configuration of PostDiag used in all experiments. For the data split, each training dataset $\mathcal{D}$ is randomly partitioned into a fitting set $\mathcal{D}_{\mathrm{tr}}$ and an audit set $\mathcal{D}_{\mathrm{aud}}$ with a 75%/25% ratio. On $\mathcal{D}_{\mathrm{aud}}$, PostDiag in the experiments of LSR-Synth constructs $M = 64$ overlapping anchor-$k$NN patches using farthest-point anchors, where $k = 64$. We suggest that the number of anchors and the neighborhood size should be scaled to the audit-set size to guarantee that every audit sample belongs to at least one patch. For fingerprint discrepancies, the three feature dimensions are combined with weights $(\omega_{\mu}, \omega_{\beta}, \omega_{\rho}) = (0.4, 0.4, 0.2)$; the smaller weight on roughness reflects its higher sensitivity to observation noise on small patches. Skeleton stability is estimated from $B = 3$ bootstrap refits, each of size $\lfloor\rho\lvert\mathcal{D}_{\mathrm{fit}}\rvert\rfloor$ points drawn with replacement, where $\rho= 0.7$. The Diagnostic Risk Plane thresholds are setted as $\tau_{B} = 0.25$ and $\tau_{V} = 0.15$ based on preliminary experiments on a small held-out subset of tasks. We use these thresholds transfers across all four benchmark domains without per-domain re-tuning.

Search budget. To evaluate the effectiveness of PostDiag under limited search resources, we adopt low-budget configurations in the main experiments. For each task, LLM-SR is limited to 100 LLM calls; LaSR runs for 20 outer iterations; SR-Scientist runs for 10 outer iterations, each with at most 25 interaction turns; PySR runs for 65 iterations; and gplearn evolves for 100 generations. Direct LLM generates 5 candidates in Low and 10 in High from a single prompt batch; SGA runs 5 rounds (80 generated candidates) in Low and 20 rounds (320 generated candidates) in High. For each pair of a base solver and its PostDiag-augmented variant, we keep the base search budget, internal hyperparameters, LLM backbone, and random seeds identical. The diagnostic and structural intervention costs introduced by PostDiag are not deducted from the base solver’s search budget, but are included in the reported end-to-end runtime. Measured results are aggregated over three random seeds. Table 1 uses Qwen3.5-27B for all LLM-based solvers and for PostDiag interventions on GP solvers; Appendix G reports Qwen3.5-122B-A10B comparison, symbolic accuracy, budget efficiency, and additional backbones, while Appendix E reports OOD results.

Results of Different SR Solvers Integrated with PostDiag

Table 1 compares seven SR solvers with and without PostDiag on LSR-Synth. Across LLM-based solvers, model scales, and domains, PostDiag consistently reduces NMSE and improves $\mathrm{Acc}_{0.01}$, indicating that stronger LLM backbones do not eliminate the need for structured diagnostic feedback. The improvements also extend to GP-based solvers, although it cannot fully compensate for the gplearn’s limited search capacity. These gains incur additional diagnostic cost from PostDiag, but require no modification to the underlying search mechanisms.

We evaluate all solvers under reduced search budgets to test whether PostDiag improves search efficiency. Larger budgets can mask differences in discovery speed as solvers approach convergence. Despite using far fewer search steps, PostDiag-enhanced solvers match or exceed results reported with much larger budgets. For example, in SR-Scientist [45], the authors report $\mathrm{Acc}_{0.01}$ of 30.49 for LLM-SR and 29.46 for their training-free variant with Qwen3-32B, using 1,000 LLM calls and 40 outer iterations. With Qwen3.5-27B, LLM-SR+PostDiag and SR-Scientist+PostDiag reach 31.0 and 35.7 using only 100 calls and 10 iterations, respectively, indicating improved search efficiency rather than gains from larger budgets.

Diagnostic Fidelity of the PostDiag

To assess whether PostDiag correctly identifies skeleton risks, we manually construct four labeled skeleton–data pair types. Settled pairs use the correct skeleton with varying Gaussian noise; Underfit pairs preserve the main structure and remove or replace one true term; Overfit pairs retain the correct skeleton and add redundant or unstable terms; and Mismatched pairs pair the data with an unrelated skeleton. Without access to these labels, PostDiag estimates Skeleton Bias and Variance via PFP and maps each candidate onto the DRP.

As shown in Figure 2a, the four pair types largely occupy their expected DRP regions. We observe a moderate overlap between underfit and mismatch near the stability boundary. This is expected, since both cases involve structural bias. When this bias is severe, coefficient refitting on resampled subsets may become unstable, making the boundary between an incomplete skeleton and a mismatched skeleton less sharp. Nevertheless, the clear quadrant structure shows that PFP and DRP provide a meaningful risk localization for different skeleton failure modes, yielding feedback with more explicit diagnostic information.

Ablation Study

To assess the contribution of different diagnostic signals, we compare three ablated variants with full PostDiag. “NMSE-only” provides only scalar candidate error to the LLM intervention module. “NMSE+DRP” infers skeleton adequacy and stability from NMSE-based signals and places the candidate on the DRP. “Evidence-only” supplies the PFP-based evidence map without DRP localization. “PostDiag” combines the evidence map, DRP localization, and explicit instructions to preserve, repair, simplify, or rewrite the skeleton. We evaluate all variants on LLM-SR, LaSR, SR-Scientist, and PySR. Results in Figure 2b are averaged across the four LSR-Synth domains. Qwen-27B serves as both the backbone for all LLM-based solvers and the intervention model for all four solvers.

As shown in Figure 2b, NMSE-only yields little improvement because the base solvers already rely on scalar fitness, so passing NMSE to the intervention model provides little additional information. NMSE+DRP improves the results, but its gains are limited. Estimating skeleton risk from a scalar NMSE is too coarse to reliably distinguish different failure modes, which can lead to inaccurate DRP localization. Evidence-only produces much larger improvements for LaSR and SR-Scientist, suggesting that these solvers benefit more from detailed information about where and how a candidate fails than from a coarse risk label. For PySR, NMSE+DRP is slightly more effective than Evidence-only, indicating that an explicit intervention operator is useful for guiding its population-based search. Full PostDiag performs best across all solvers, confirming that PFP, the evidence map, and DRP are complementary: PFP localizes risk, the evidence map reveals local defects, and DRP selects the appropriate structural intervention.

Conclusion and Limitations

This work addresses a feedback bottleneck shared by diverse SR solvers: scalar fitness can rank candidates but offers little guidance for structural refinement. Rather than developing a new solver, PostDiag separates feedback design from search design. Using only a candidate skeleton and observed data, it converts patch-wise behavior into expressiveness and stability risks, localized evidence, and targeted interventions. Experiments across LLM-based and GP-based solvers show that PostDiag improves accuracy and accelerates convergence without accessing or redesigning their internal search processes. These results establish post-diagnosis as an orthogonal way to strengthen existing SR solvers while preserving their distinctive search mechanisms.

Despite its search-independent diagnostic core, our current implementation is evaluated only on LLM-based and GP-based solvers, which provide natural interfaces for accepting revised candidates. Solvers without such interfaces, such as reinforcement-learning-based methods [34], may require more substantial adaptation. Future work will extend the diagnostic interface to broader search paradigms and improve the efficiency of diagnosis and intervention.

The Use of Large Language Models (LLMs)

Beyond the use of LLMs in the method described in this paper, LLMs served only as auxiliary tools for generating some of the cartoon characters in Figure 1 and for polishing and refining the manuscript’s language. Their contributions to research conception, experimental design, or manuscript preparation do not warrant authorship or contributor status. All substantive research ideas, methods, and results presented in this work were independently developed or produced by the authors.

Dataset Details

Task Composition

We conduct all experiments on LSR-Synth, the synthetic subset of LLM-SRBench (Shojaee et al. [39]). Following the evaluation setting used by SR-Scientist (Xia et al. [45]), the benchmark contains 129 symbolic regression tasks from four scientific domains: 44 in physics, 25 in material science, 36 in chemistry, and 24 in biology.

Representative tasks include damped oscillators, stress–strain relations, chemical kinetics, and population dynamics. Each task provides numerical observations together with the variable names and their scientific meanings. Its ground-truth expression combines established domain-specific terms with novel synthetic terms. Consequently, a solver cannot rely solely on memorized canonical equations and must instead use the observations to search for previously unseen structures. Compared with LSR-Transform, which may be more susceptible to memorization from pretraining, LSR-Synth is therefore better suited to evaluating data-driven equation discovery.

Data Generation and Splits

The construction of LSR-Synth consists of scientific-problem selection, known- and synthetic-term generation, parameter instantiation, numerical solvability and novelty checks, data generation, and expert validation. For each scientific problem, candidate known and synthetic terms are first generated and combined into complete expressions, whose coefficients are instantiated within scientifically meaningful ranges. For the dynamical systems in chemistry, biology, and physics, data are generated through numerical integration using SciPy’s solve_ivp function with the RK45 method. The static stress–strain relationships in material science are evaluated directly over predefined input ranges. Candidates that cannot be solved reliably, produce divergent or degenerate behavior, or are excessively similar to existing tasks are filtered out. The retained expressions and visualizations of their generated data are then cross-checked by domain experts.

Each task contains 5,000 observations. Chemistry, biology, and physics tasks are sampled uniformly over the time interval $t \in[0,60]$, whereas material-science tasks cover strain $\epsilon\in[0,0.6]$ and temperature $T \in[273,573]$ K. For time-dependent tasks, the 500 observations with the largest time values form the out-of-distribution (OOD) test set. For material-science tasks, the 500 observations with the largest temperature values form the OOD test set. Of the remaining 4,500 observations, 500 are randomly selected as the in-distribution (ID) test set, and the other 4,000 observations constitute the training set provided to the solver.

Data-Usage Protocol for PostDiag

All diagnosis and intervention operations performed by PostDiag are restricted to the training split and are strictly isolated from both test sets. Under the default configuration described in the main text, the training data for each task are further divided into a fitting set $\mathcal{D}_{\mathrm{tr}}$ and an audit set $\mathcal{D}_{\mathrm{aud}}$ using a 75%/25% split. Candidate coefficient fitting and bootstrap refitting use only $\mathcal{D}_{\mathrm{tr}}$, whereas the PPF fingerprints, skeleton adequacy and stability scores, and evidence maps are computed on $\mathcal{D}_{\mathrm{aud}}$. The acceptance check for an intervened candidate also uses only the training split.

The ID and OOD test sets are used solely to evaluate the final expression after search terminates. They are never used for candidate generation, diagnosis, structural intervention, rewrite acceptance, or candidate selection. Although an LLM-based base solver may consume the task and variable descriptions supplied by the benchmark, PostDiag itself receives only a candidate skeleton and the numerical observations from the training split.

Base Solver Configurations

Recent progress has led to increasing interest in symbolic regression (SR) ([3, 47, 29, 51, 23, 17, 24, 45, 38]). Unlike conventional predictive models, SR seeks to recover explicit mathematical relationships from observations, making the resulting equations directly interpretable and scientifically testable across domains such as physics, biology, and materials science [50, 22, 13, 11, 14, 49, 48]. A wide variety of symbolic regression methods have been developed in recent years. In this work, we integrate PostDiag with two of the most common SR paradigms: GP-based and LLM-based methods. Across these two families, we select seven widely used SR solvers as base solvers and incorporate PostDiag into each of them. We next describe the search budgets and integration strategies used for these solvers.

Search-Budget Protocol

The low-budget settings used in the main experiments are specified in § 4.1, while the high-budget reference settings are described in Appendix § 13.§ 3. To avoid repetition, this section does not restate the limits on LLM calls, outer iterations, or evolutionary generations. Instead, we report the internal configurations that remain fixed across budget levels and between each base solver and its PostDiag-augmented counterpart. Unless otherwise specified, we build on the official implementation of each method and keep its prompt templates, search operators, coefficient-fitting procedure, and candidate-evaluation protocol unchanged. The detailed integration of PostDiag with each solver is provided in Appendix § 12.

The Low and High settings retain each solver’s native search unit. Within each base/PostDiag pair, the native budget is identical, while PostDiag overhead is included in the reported end-to-end runtime.

Solver	Low	High	Native budget unit
Direct LLM	5	10	generated candidates
SGA	5 (80)	20 (320)	rounds (generated candidates)
LLM-SR	100	1,000	LLM calls
LaSR	20	40	outer iterations
SR-Scientist	$10 \times\leq25$	$40 \times\leq25$	outer iterations $\times$ turns
PySR	65	130	iterations
gplearn	100	200	generations

LLM-Based Solvers

Direct LLM. We follow the DataBlind Direct Prompting protocol of LLM-SRBench [39]. The LLM proposes equation skeletons from the scientific context without observing the numerical samples, while the numerical constants in each proposed skeleton are fitted by the benchmark evaluation pipeline. The PostDiag variant retains this direct generator and adds data-driven diagnosis and intervention after candidate evaluation; thus, DataBlind describes the base generator rather than the complete augmented pipeline. Low and High generate 5 and 10 candidates, respectively, from a single prompt batch.

SGA. SGA [28] uses a bilevel procedure in which an LLM proposes discrete scientific hypotheses and simulation-based optimization fits their continuous parameters. Its PostDiag variant attaches diagnostic feedback outside this bilevel loop, preserving the underlying hypothesis generation and continuous parameter optimization. Low uses 5 rounds and High uses 20 rounds, with 16 generated candidate attempts per round; the continuous-parameter optimizer is unchanged across budgets.

LLM-SR. LLM-SR [38] represents candidate equations as Python programs containing placeholders for numerical constants. An LLM proposes program skeletons conditioned on the task description and high-performing programs retrieved from an experience pool, while the numerical constants are fitted using BFGS. We retain the prompt templates, multi-island experience pool, candidate-ranking procedure, and coefficient-optimization pipeline of the official implementation. Apart from the search budgets specified elsewhere, all configurations are shared across LLM backbones and between the base and PostDiag-augmented variants.

LaSR. LaSR [10] maintains an LLM-induced library of natural-language concepts on top of a multi-population genetic-programming search. We use 10 populations with 33 individuals per population and run 550 GP cycles within each outer iteration. The maximum expression size is 30, and the concept library contains at most 20 concepts. The weights assigned to LLM-assisted crossover, mutation, and randomization are all set to 0.02. The symbolic operators and nested constraints follow the PySR configuration described below.

SR-Scientist. SR-Scientist [45] equips an LLM agent with data-analysis and equation-evaluation tools, allowing it to write and execute code while iteratively improving candidate equations. We retain its official tool interfaces, experience buffer, and candidate-evaluation procedure. Each search iteration uses a target MAPE of 0.1% and terminates early when the error reaches $10^{-4}\%$. The LLM backbone and generation settings are held fixed between each base/PostDiag pair.

GP-Based Solvers

PySR. PySR [5] performs island-based genetic programming using a Julia backend, together with expression-complexity regularization and numerical constant optimization. We use 15 populations with 33 individuals per population and run 550 GP cycles within each iteration. The maximum expression size is 50. The binary operator set is $\{+, -, \times, \div, \operatorname{pow}\}$, and the unary operator set is $\{\exp, \log, \sqrt{\cdot}, \sin, \cos\}$.

gplearn. gplearn [40] implements single-population tree-based genetic programming. We use a population size of 500 and a tournament size of 20. The function set contains addition, subtraction, multiplication, division, sine, cosine, square root, logarithm, absolute value, negation, and inverse. Other parameters, including the parsimony coefficient and the crossover and mutation probabilities, retain the defaults of the official implementation.

LLM Backbones

Model Selection

Our experiments use Qwen3.5-27B and Qwen3.5-122B-A10B as the principal LLM backbones [35]. Qwen3.5-27B is reported in the main table, while Qwen3.5-122B-A10B is reported in the appendix. Qwen3.5-27B is a 27B-parameter dense model, whereas Qwen3.5-122B-A10B adopts a mixture-of-experts architecture with 122B total parameters and approximately 10B active parameters. Although both models belong to the same family, they differ in capacity and architecture, allowing us to evaluate whether PostDiag remains effective across different backbone configurations.

We further report supplementary results with Qwen3.5-9B and gpt-oss-20b [1] in the appendix. Qwen3.5-9B evaluates PostDiag with a smaller model and, together with the two principal Qwen3.5 backbones, provides coverage across different model capacities. The gpt-oss-20b experiments provide complementary evidence from a different model family and help assess whether PostDiag depends on a particular family of LLMs.

Backbone Consistency

For the five LLM-based base solvers, LLM-SR, LaSR, SR-Scientist, Direct LLM, and SGA, each base method and its PostDiag-augmented counterpart use the same LLM backbone. For example, when a base solver uses Qwen3.5-27B, the corresponding PostDiag structural intervention is also performed by Qwen3.5-27B. This avoids confounding the effect of PostDiag with differences in backbone capability and keeps the primary distinction between each pair to whether PostDiag is applied.

For PySR and gplearn, whose original search procedures do not rely on an LLM, the base search remains unchanged. When PostDiag is added, we use Qwen3.5-27B consistently to perform the diagnosed structural interventions.

PostDiag Configuration

Each base solver retains its original configuration. Between a base method and its PostDiag-augmented counterpart, we keep the base search budget, internal hyperparameters, and random seeds fixed so that the comparison focuses on the diagnosis and structural intervention introduced by PostDiag.

For each PostDiag structural intervention, we set the sampling temperature to 0.3 and the maximum completion length to 8,192 tokens. Conditioned on the DRP state and Evidence Map, the model returns a structured preserve-or-revise decision, a revised expression skeleton, and the corresponding edit type. The LLM proposes only the structural revision and does not determine the final numerical coefficients. After parsing and validity checking, the coefficients of the resulting expression are refitted. We use the same PostDiag configuration for all base solvers and LLM backbones.

Compute Resources

Experiments are conducted on a compute cluster equipped with 16 NVIDIA A800 GPUs and 8 NVIDIA H800 GPUs. The LLM backbones are served locally using vLLM. PySR, gplearn, and the diagnostic computation of PostDiag run primarily on CPUs, while GPUs are used mainly for LLM inference. For each base solver, all reported runtime configurations—including Low, Low+PostDiag, and High—use the same GPU type, GPU count, vLLM serving configuration, and server-side concurrency limit. Each method retains its native request-scheduling pattern: LLM-SR issues up to four concurrent requests, whereas the interaction trajectory of SR-Scientist is executed sequentially; the other solvers follow their official implementations.

We report the mean end-to-end wall-clock time per task, measured from the start of the base solver’s search until it returns the final expression. The reported time includes base search, candidate evaluation, and coefficient fitting. For PostDiag-augmented variants, it also includes the overhead of diagnostic computation, coefficient refitting, and additional LLM-based structural intervention. All runtime measurements are averaged over three random seeds.

Out-of-Distribution Performance

LSR-Synth provides an independent out-of-distribution (OOD) test set for each task, with samples extending beyond the training range along a control variable such as time or temperature. The detailed split construction is provided in Appendix § 7. Following SR-Scientist (Xia et al., [45], we report OOD $\mathrm{Acc}_{0.01}$. After search terminates, the final expression returned by each method is evaluated directly on the OOD test set, without additional search, coefficient refitting, or candidate selection. The OOD data remain unseen throughout search, diagnosis, and structural intervention.

As shown in Table 2, PostDiag generally improves OOD accuracy across different base solvers and LLM backbones. With Qwen3.5-27B, the overall OOD $\mathrm{Acc}_{0.01}$ of LLM-SR increases from 21.2% to 28.9%, while that of SR-Scientist increases from 24.8% to 33.6%. With Qwen3.5-122B-A10B, the overall OOD accuracy of SR-Scientist increases from 40.6% to 43.8%. PySR also improves from 27.1% to 31.0%, suggesting that the extrapolation benefits of PostDiag are not limited to LLM-based solvers. In contrast, gplearn and its PostDiag-augmented variant both obtain zero OOD accuracy, indicating that diagnosis and structural intervention cannot fully compensate for insufficient search capability in the base solver.

Backbone/category	Method	Overall	Phys	MatSci	Chem	Bio
Qwen3.5-27B
LLM-based	Direct LLM [39]	1.0	0.8	2.7	0.9	0.0
LLM-based	Direct LLM [39]+PostDiag	2.3	2.3	2.7	2.8	1.4
LLM-based	SGA [28]	9.6	6.1	9.3	16.7	5.6
LLM-based	SGA [28]+PostDiag	16.3	12.9	21.3	21.3	9.7
LLM-based	LLM-SR	21.2	12.9	41.3	26.9	6.9
LLM-based	LLM-SR+PostDiag	28.9	20.5	61.3	27.8	12.5
LLM-based	LaSR	12.4	10.6	16.0	15.7	6.9
LLM-based	LaSR+PostDiag	19.1	17.4	30.7	18.5	11.1
LLM-based	SR-Scientist	24.8	16.7	50.7	25.0	12.5
LLM-based	SR-Scientist+PostDiag	33.6	22.7	64.0	36.1	18.1
Qwen3.5-122B-A10B
LLM-based	Direct LLM [39]	1.8	2.3	1.3	1.9	1.4
LLM-based	Direct LLM [39]+PostDiag	3.1	3.8	2.7	2.8	2.8
LLM-based	SGA [28]	11.9	8.3	12.0	18.5	8.3
LLM-based	SGA [28]+PostDiag	18.3	15.2	21.3	24.1	12.5
LLM-based	LLM-SR	26.6	18.9	53.3	25.9	13.9
LLM-based	LLM-SR+PostDiag	34.1	24.2	66.7	32.4	20.8
LLM-based	LaSR	14.0	13.6	18.7	15.7	6.9
LLM-based	LaSR+PostDiag	21.7	20.5	34.7	20.4	12.5
LLM-based	SR-Scientist	40.6	27.3	78.7	40.7	25.0
LLM-based	SR-Scientist+PostDiag	48.3	35.6	81.3	50.0	34.7
GP-based methods
GP-based	PySR	27.1	20.5	68.0	15.7	13.9
GP-based	PySR+PostDiag	31.0	22.7	74.7	20.4	16.7
GP-based	gplearn	0.0	0.0	0.0	0.0	0.0
GP-based	gplearn+PostDiag	0.0	0.0	0.0	0.0	0.0

The gains remain substantial under OOD evaluation, indicating that the benefits of PostDiag are largely preserved on out-of-distribution data. In material science, OOD accuracy is generally higher than the corresponding ID accuracy for both the base and PostDiag-augmented methods.

Details of Integration with Base Solvers

PostDiag is designed to be attached to existing SR solvers through lightweight adapters. We describe the integration recipe for each of the seven representative solvers in this section. The adapters handle representation conversion, diagnostic feedback, and candidate-pool insertion at the solver boundary, while the solver’s core candidate generation, selection, and search-update rules remain unchanged.

Direct LLM [39] generates equation skeletons from the scientific context without observing numerical samples; numerical constants are fitted by the benchmark pipeline. In our PostDiag, generated skeleton is diagnosed and revised. The accepted revisions enter the selected candidate pool, from which the final result is saved. The base generator remains DataBlind, while the augmented candidate evaluation is explicitly data-driven.

SGA [28] alternates LLM-based discrete hypothesis generation with simulation-based continuous parameter optimization. The proposed adapter parses each generated equation into SGA’s trainable PyTorch module, fits its parameters with the unchanged inner optimizer, and adds valid fitted candidates to the candidate pool. PostDiag diagnoses a fitted candidate, proposes a structural revision, and sends the parsed revision through the same module construction, fitting, and evaluation path before it becomes eligible for the selected pool.

LLM-SR [38] treats equation discovery as evolutionary program search, where an LLM proposes candidate skeletons conditioned on in-context examples drawn from a buffer of high-performing programs. PostDiag integrates through this buffer interface: for every newly admitted elite candidate $s$, PostDiag produces its diagnostic label and evidence map; the rewritten skeleton $\tilde{s}$ together with a brief textual annotation of the diagnosed failure mode is stored. Subsequent prompts draw both score-ranked and diagnosis-annotated examples, biasing the LLM’s next proposals toward the identified failure mode.

LaSR (Graveli et al., 2024) couples LLMs with genetic programming by periodically abstracting high-performing hypotheses into natural-language concepts, which then guide subsequent symbolic mutation. PostDiag integrates at the concept abstraction stage: when LaSR queries the LLM for concept induction, we additionally supply the diagnostic label and evidence map of the input hypotheses, so that the abstracted concepts carry structural failure-mode information. The resulting diagnosis-aware concepts guide the subsequent GP search toward diagnosed structural fixes, while LaSR’s underlying GP kernel remains untouched.

SR-Scientist (Xia et al. [45]) frames equation discovery as an autonomous agent that reads data, writes code to test hypotheses, and iterates on candidates through a tool-use loop that already accepts arbitrary numerical feedback signals. PostDiag integrates as an additional diagnose(candidate) tool that the agent can invoke after any candidate evaluation; the tool returns the structured diagnostic report and evidence map, which are then included in the agent’s next reasoning step. The integration only exposes a new tool endpoint and does not modify SR-Scientist’s control flow.

PySR (Cranmer [5]) performs island-based genetic programming with configurable custom operators and supports external individual injection through its population interface. PostDiag integrates through diagnosis-driven immigration: at the end of each generation, PostDiag diagnoses the best-performing individual. The candidate is revised through diagnosis-guided intervention and injected into the initial population of the next generation as an immigrant, thereby steering the subsequent evolutionary search. Since PySR natively supports the insertion of external individuals, the integration requires only a lightweight adapter that converts Python expressions into expression trees, without modifying its selection, mutation, or crossover operators.

gplearn (Stephens [40]) implements classical tree-based genetic programming and exposes init_method and warm_start for warm-starting the population from previously fitted regressors. PostDiag integrates by feeding the rewritten token $\tilde{s}$ as seeded individuals into the warm-started initial population of the next generation. Since gplearn’s tree grammar covers all operators used in PostDiag, the adapter is a thin conversion between our program representation and gplearn’s internal tree encoding, and does not alter gplearn’s selection or reproduction operators.

Across all seven solvers, PostDiag exposes the same standardized output, $(\mathcal{G}_s, \tilde{s})$. With only lightweight implementation changes, it equips each solver with post-diagnosis capabilities while leaving the core search algorithm intact. PostDiag therefore represents a complementary technical direction that is orthogonal to the development of new SR algorithms.

Additional Experimental Results

Unless otherwise specified, we use Qwen3.5-27B as the default experimental backbone, while Qwen3.5-122B-A10B is included to evaluate PostDiag with a stronger model. Since PySR and gplearn do not natively use an LLM, their PostDiag-augmented variants consistently use the default Qwen3.5-27B backbone.

Low-Budget Results with Qwen3.5-122B-A10B

Table 3 reports the counterpart of the main ID comparison using Qwen3.5-122B-A10B. Each base/PostDiag pair retains the same search configuration and random seeds, with the LLM backbone shared between the base solver and the PostDiag intervention.

Method / setting	Overall				Phys			MatSci			Chem			Bio
	Time (s)	NMSE	A.01	A.001	NMSE	A.01	A.001	NMSE	A.01	A.001	NMSE	A.01	A.001	NMSE	A.01	A.001
Direct LLM (Shojaee et al. [39])	35.3	4.62E-02	2.8	1.3	3.82E-02	3.0	0.8	5.78E-02	2.7	1.3	2.26E-02	3.7	1.9	8.42E-02	1.4	1.4
Direct LLM (Shojaee et al. [39])+PostDiag	39.3	9.54E-03	5.4	3.4	8.64E-03	6.1	3.8	2.86E-03	5.3	4.0	6.43E-03	4.6	2.8	2.28E-02	5.6	2.8
SGA (Ma et al. [28])	368.0	8.44E-03	14.7	8.8	5.88E-03	12.1	7.6	6.17E-03	13.3	8.0	5.10E-03	22.2	13.0	2.05E-02	9.7	5.6
SGA (Ma et al. [28])+PostDiag	406.2	1.44E-03	22.7	12.7	1.52E-03	20.5	12.9	2.06E-04	25.3	10.7	7.32E-04	28.7	15.7	3.66E-03	15.3	9.7
LLM-SR	201.1	9.68E-04	27.6	12.7	2.41E-03	21.2	9.8	1.91E-05	50.7	26.7	1.18E-04	27.8	9.3	5.86E-04	15.3	8.3
LLM-SR+PostDiag	231.7	1.31E-04	36.4	19.9	2.13E-04	28.0	15.9	5.82E-06	65.3	38.7	3.86E-05	36.1	15.7	2.47E-04	22.2	13.9
LaSR	716.3	8.71E-03	15.8	9.3	1.47E-02	14.4	7.6	1.92E-03	17.3	9.3	5.60E-03	19.4	13.0	9.45E-03	11.1	6.9
LaSR+PostDiag	762.1	6.76E-04	24.5	13.7	1.19E-03	23.5	10.6	3.76E-05	33.3	16.0	4.96E-04	25.0	17.6	6.68E-04	16.7	11.1
SR-Scientist	833.6	5.96E-04	41.3	24.3	1.43E-03	29.5	16.7	2.03E-05	74.7	50.7	8.12E-05	41.7	21.3	4.37E-04	27.8	15.3
SR-Scientist+PostDiag	902.2	1.04E-04	48.8	31.3	1.76E-04	36.4	22.7	6.08E-06	78.7	58.7	3.04E-05	50.9	29.6	1.82E-04	37.5	20.8

At the Overall level, PostDiag lowers NMSE and improves both accuracy thresholds for all solvers under Qwen3.5-122B-A10B. The stronger backbone therefore continues to benefit from diagnostic feedback, while the final-accuracy gap between Direct LLM and the iterative solvers remains, indicating that the upper of performance still depends on the underlying search procedure.

Symbolic Accuracy

We additionally evaluate whether the final expression is symbolically equivalent to the ground-truth expression. Following LLM-SRBench, we remove fitted parameters and constants before judging equivalence so that the metric measures recovery of the functional structure rather than numerical coefficient matching. We then adopt the more conservative two-stage adjudication procedure of SR-Scientist: GPT-OSS-120B evaluates each pair independently ten times, unanimous decisions are retained, and non-unanimous cases are resolved by human evaluators. Failed runs and invalid final expressions are counted as incorrect. Table 4 reports the pooled symbolic accuracy over all tasks and three seeds.

Base solver	Backbone	Base SA	+PostDiag SA	Δ SA
Direct LLM	Qwen3.5-27B	1.55	1.55	0.00
	Qwen3.5-122B-A10B	2.07	2.58	+0.52
SGA	Qwen3.5-27B	3.62	4.65	+1.03
	Qwen3.5-122B-A10B	3.88	5.17	+1.29
LLM-SR	Qwen3.5-27B	3.88	5.68	+1.81
	Qwen3.5-122B-A10B	4.39	6.46	+2.07
LaSR	Qwen3.5-27B	1.81	3.62	+1.81
	Qwen3.5-122B-A10B	2.58	4.91	+2.33
SR-Scientist	Qwen3.5-27B	4.91	7.24	+2.33
	Qwen3.5-122B-A10B	6.72	9.56	+2.84
PySR	N/A	4.13	5.68	+1.55
gplearn	N/A	0.00	0.00	0.00

PostDiag raises pooled SA for most solver–backbone pairs, with gains ranging from 0.52 points for Direct LLM under Qwen3.5-122B-A10B to 2.84 points for SR-Scientist under the same backbone. SGA gains 1.03 and 1.29 percentage points under Qwen3.5-27B and Qwen3.5-122B-A10B, respectively. Direct LLM under Qwen3.5-27B remains at 1.55%, and gplearn remains at 0.00%.

Budget-Efficiency Analysis

We compare PostDiag-augmented solvers against their unaugmented counterparts under the Low and High settings defined in Appendix B.1. The High settings are 10 generated candidates for Direct LLM [39]; 20 rounds, or 320 generated candidate attempts, for SGA [28]; 1,000 LLM calls for LLM-SR; 40 outer iterations for LaSR; 40 outer iterations with at most 25 interaction turns each for SR-Scientist; 130 iterations for PySR; and 200 generations for gplearn. Low+PostDiag retains the same native base-solver budget as Low. We report accuracy to tolerance as the primary criterion, arithmetic mean NMSE, and mean end-to-end wall-clock time per task run. Tables 5, 6, and 7 evaluate solution quality at matched native budgets and compare low-budget PostDiag variants with the corresponding high-budget base solvers.

The Low+PostDiag and High comparisons show that diagnosis and additional search are complementary rather than interchangeable. Direct LLM+PostDiag under Qwen3.5-27B exceeds High in Overall accuracy at similar runtime, whereas under Qwen3.5-122B-A10B it has higher accuracy

Method / setting	Overall Time (s)	Overall NMSE	Overall $A_{.01}$	Overall $A_{.001}$	Phys NMSE	Phys $A_{.01}$	Phys $A_{.001}$	MatSci NMSE	MatSci $A_{.01}$	MatSci $A_{.001}$	Chem NMSE	Chem $A_{.01}$	Chem $A_{.001}$	Bio NMSE	Bio $A_{.01}$	Bio $A_{.001}$
Direct LLM [39]	30.5	7.72E-02	2.6	0.8	7.32E-02	3.0	0.8	6.24E-02	2.7	1.3	3.50E-02	2.8	0.9	1.63E-01	1.4	0.0
Direct LLM [39]+PostDiag	33.7	3.23E-02	4.4	2.3	5.22E-02	5.3	3.8	4.34E-03	4.0	2.7	8.76E-03	3.7	0.9	6.03E-02	4.2	1.4
Direct LLM / High	45.5	5.78E-02	3.4	1.8	5.86E-02	3.0	1.5	3.98E-02	5.3	2.7	2.80E-02	3.7	1.9	1.20E-01	1.4	1.4
SGA [28]	320.0	1.39E-02	11.9	5.9	9.43E-03	8.3	6.1	8.21E-03	12.0	5.3	4.92E-03	19.4	7.4	4.15E-02	6.9	4.2
SGA [28]+PostDiag	350.2	2.40E-03	18.6	9.8	2.66E-03	14.4	9.8	1.87E-04	25.3	6.7	8.65E-04	24.1	13.0	6.52E-03	11.1	8.3
SGA / High	1250.0	1.70E-03	19.9	10.6	1.82E-03	15.9	12.1	2.27E-04	22.7	8.0	4.33E-04	26.9	13.0	4.89E-03	13.9	6.9
LLM-SR / Low	167.5	1.45E-03	22.0	8.8	3.66E-03	14.4	8.3	3.12E-05	37.3	13.3	1.75E-04	27.8	8.3	8.15E-04	11.1	5.6
LLM-SR+PD / Low	191.7	2.12E-04	31.0	15.2	3.48E-04	23.5	12.1	8.94E-06	60.0	30.7	5.84E-05	28.7	12.0	4.02E-04	18.1	9.7
LLM-SR / High	1603.4	3.65E-04	34.1	13.2	7.82E-04	28.0	12.9	1.20E-05	53.3	22.7	7.21E-05	37.0	11.1	4.08E-04	20.8	6.9
LaSR / Low	559.2	1.06E-02	13.4	7.5	1.73E-02	11.4	6.8	2.48E-03	14.7	6.7	6.90E-03	18.5	11.1	1.22E-02	8.3	4.2
LaSR+PD / Low	595.4	8.62E-04	21.4	12.7	1.52E-03	20.5	11.4	4.82E-05	29.3	14.7	6.45E-04	22.2	14.8	8.29E-04	13.9	9.7
LaSR / High	1118.4	8.13E-03	15.5	8.8	1.36E-02	14.4	7.6	1.73E-03	17.3	9.3	5.10E-03	19.4	12.0	9.30E-03	9.7	5.6
SR-Scientist / Low	648.0	1.12E-03	25.6	11.9	2.74E-03	17.4	9.8	4.47E-05	46.7	21.3	1.49E-04	27.8	11.1	7.26E-04	15.3	6.9
SR-Scientist+PD / Low	702.2	1.92E-04	35.7	20.4	3.11E-04	25.8	15.9	1.08E-05	62.7	37.3	6.47E-05	38.0	20.4	3.54E-04	22.2	11.1
SR-Scientist / High	2547.1	2.23E-04	37.7	21.2	4.47E-04	26.5	15.9	1.31E-05	70.7	48.0	5.03E-05	38.0	16.7	2.91E-04	23.6	9.7

Method / setting	Overall Time (s)	Overall NMSE	Overall $A_{.01}$	Overall $A_{.001}$	Phys NMSE	Phys $A_{.01}$	Phys $A_{.001}$	MatSci NMSE	MatSci $A_{.01}$	MatSci $A_{.001}$	Chem NMSE	Chem $A_{.01}$	Chem $A_{.001}$	Bio NMSE	Bio $A_{.01}$	Bio $A_{.001}$
Direct LLM [39]	35.3	4.62E-02	2.8	1.3	3.82E-02	3.0	0.8	5.78E-02	2.7	1.3	2.26E-02	3.7	1.9	8.42E-02	1.4	1.4
Direct LLM [39]+PostDiag	39.9	5.94E-03	5.4	3.4	8.64E-03	6.1	3.8	2.86E-03	5.3	4.0	6.43E-03	4.6	2.8	2.29E-02	5.6	2.8
Direct LLM / High	52.7	1.55E-03	3.9	2.1	3.86E-03	3.0	2.3	3.06E-05	6.7	2.7	1.89E-04	3.7	1.9	9.36E-04	2.8	1.4
SGA [28]	368.0	8.44E-03	14.7	8.8	5.88E-03	12.1	7.6	6.17E-03	13.3	8.0	5.10E-03	22.2	13.0	2.05E-02	9.7	5.6
SGA [28]+PostDiag	406.2	1.44E-03	22.7	12.7	1.52E-03	20.5	12.9	6.0E-04	25.3	10.7	7.32E-04	28.7	15.7	3.65E-03	15.3	9.7
SGA / High	1442.0	7.52E-04	26.9	15.0	1.43E-03	25.0	14.4	3.65E-05	28.0	18.7	3.89E-04	34.7	17.7	7.79E-04	20.8	9.7
LLM-SR / Low	201.1	9.68E-04	27.6	12.7	2.41E-03	21.2	9.8	1.91E-05	50.7	26.7	1.18E-04	27.8	9.3	5.86E-04	15.3	8.3
LLM-SR+PD / Low	231.7	1.31E-04	36.4	19.9	2.13E-04	28.0	15.9	5.82E-06	65.3	38.7	3.86E-05	36.1	15.7	2.47E-04	22.2	13.9
LLM-SR / High	1939.2	2.99E-04	39.3	16.8	5.15E-04	33.3	14.4	3.54E-06	66.7	34.7	4.42E-05	37.0	12.0	3.59E-04	25.0	9.7
LaSR / Low	716.3	8.71E-03	15.8	9.3	1.47E-02	14.4	7.6	1.92E-03	17.3	9.3	5.60E-03	19.4	13.0	9.45E-03	11.1	6.9
LaSR+PD / Low	762.1	6.76E-04	24.5	13.7	1.19E-03	23.5	10.6	3.76E-05	33.3	16.0	4.96E-04	25.0	17.6	6.68E-04	16.7	11.1
LaSR / High	1432.6	5.74E-03	19.1	10.9	4.16E-03	17.4	8.3	6.24E-04	22.7	12.0	3.14E-03	22.2	13.9	8.04E-03	13.9	9.7
SR-Scientist / Low	833.6	5.96E-04	41.3	24.3	1.43E-03	29.5	16.7	2.03E-05	74.7	50.7	8.12E-05	41.7	21.3	4.37E-04	27.8	15.3
SR-Scientist+PD / Low	902.2	1.04E-04	48.8	31.3	1.76E-04	36.4	22.7	6.08E-06	78.7	58.7	3.04E-05	50.9	29.6	1.82E-04	37.5	20.8
SR-Scientist / High	3289.4	9.21E-05	48.1	31.0	2.33E-04	34.1	21.2	2.95E-06	86.7	68.0	4.71E-06	47.2	25.9	5.75E-05	34.7	18.1

Method / setting	Overall Time (s)	Overall NMSE	Overall $A_{.01}$	Overall $A_{.001}$	Phys NMSE	Phys $A_{.01}$	Phys $A_{.001}$	MatSci NMSE	MatSci $A_{.01}$	MatSci $A_{.001}$	Chem NMSE	Chem $A_{.01}$	Chem $A_{.001}$	Bio NMSE	Bio $A_{.01}$	Bio $A_{.001}$
PySR / Low	8.8	2.67E-03	23.8	11.4	4.11E-03	21.2	12.1	1.70E-04	44.0	18.7	8.93E-04	20.4	10.2	5.29E-03	12.5	4.2
PySR+PD / Low	10.5	1.34E-03	28.7	13.2	1.82E-03	24.2	14.4	2.25E-05	50.7	17.3	5.29E-04	25.9	12.0	2.84E-03	18.1	8.3
PySR / High	15.5	1.76E-03	28.4	13.7	2.61E-03	25.8	15.2	8.12E-05	45.3	20.0	6.18E-04	25.0	12.0	3.66E-03	16.7	6.9
gplearn / Low	95.5	1.34E-01	0.0	0.0	8.51E-02	0.0	0.0	2.36E-02	0.0	0.0	6.56E-02	0.0	0.0	4.42E-01	0.0	0.0
gplearn+PD / Low	101.7	3.07E-02	1.0	0.0	4.29E-03	2.3	0.0	9.6E-03	0.0	0.0	4.29E-02	0.9	0.0	8.58E-02	0.0	0.0
gplearn / High	189.1	1.12E-01	0.0	0.0	6.73E-02	0.0	0.0	1.92E-02	0.0	0.0	5.31E-02	0.0	0.0	3.77E-01	0.0	0.0

but worse NMSE, exposing a metric tradeoff. For SGA under both backbones, Low+PostDiag approaches the High setting’s quality with substantially less end-to-end time, while High retains higher accuracy and lower NMSE. The remaining solvers show similarly method-dependent balances between diagnostic overhead and the benefit of a larger native search budget.

Supplementary Results with Additional LLM Backbones

To examine PostDiag under different backbones, we additionally evaluate Qwen3.5-9B and gpt-oss-20b under the same low-budget protocol. Table 8 reports the ID results for the LLM-based solvers and their PostDiag-augmented variants. All comparisons within a backbone use the same tasks, random seeds, search budgets, coefficient-fitting procedure, and validation protocol.

Category	Method / setting	Overall Time (s)	Overall NMSE	Overall A.01	Overall A.001	Phys NMSE	Phys A.01	Phys A.001	MatSci NMSE	MatSci A.01	MatSci A.001	Chem NMSE	Chem A.01	Chem A.001	Bio NMSE	Bio A.01	Bio A.001
Qwen3.5-9B
LLM-based	LLM-SR	129.8	2.36E-03	14.0	5.7	5.32E-03	9.8	3.0	2.38E-05	37.3	17.3	4.96E-04	6.5	2.8	2.14E-03	8.3	2.8
LLM-based	LLM-SR+PostDiag	146.8	3.93E-04	24.0	10.3	7.32E-04	18.2	6.8	1.20E-05	46.7	25.3	7.15E-05	22.2	8.3	6.52E-04	13.9	4.2
LLM-based	LaSR	382.5	2.16E-02	10.1	6.5	2.24E-02	6.8	4.5	4.15E-03	14.7	9.3	9.40E-03	13.0	9.3	5.64E-02	6.9	2.8
LLM-based	LaSR+PostDiag	407.9	3.27E-03	18.6	10.3	6.01E-03	13.6	6.1	6.26E-05	36.0	18.7	8.86E-04	18.5	12.0	5.16E-03	9.7	6.9
LLM-based	SR-Scientist	439.2	2.00E-03	13.7	4.1	4.48E-03	10.6	4.5	6.31E-05	33.3	10.7	4.12E-04	7.4	1.9	1.86E-03	8.3	0.0
LLM-based	SR-Scientist+PostDiag	477.2	5.39E-04	22.5	8.8	9.84E-04	15.9	6.1	1.69E-05	44.0	20.0	1.33E-04	21.3	7.4	8.76E-04	13.9	4.2
gpt-oss-20b
LLM-based	LLM-SR	244.7	1.31E-03	26.1	9.6	3.18E-03	18.9	4.5	2.62E-05	52.0	33.3	2.03E-04	22.2	3.7	8.91E-04	18.1	2.8
LLM-based	LLM-SR+PostDiag	260.9	2.39E-04	35.4	15.2	4.31E-04	25.8	9.1	8.52E-06	68.0	45.3	6.47E-05	31.5	8.3	3.87E-04	25.0	5.6
LLM-based	LaSR	395.4	1.21E-02	10.6	5.4	1.92E-02	9.1	5.3	3.18E-03	13.3	5.3	8.20E-03	13.0	7.4	1.43E-02	6.9	2.8
LLM-based	LaSR+PostDiag	419.6	9.53E-04	16.8	9.6	1.71E-03	14.4	8.3	5.04E-05	18.7	10.7	6.83E-04	19.4	11.1	9.13E-04	15.3	8.3
LLM-based	SR-Scientist	134.9	8.89E-04	33.3	15.5	2.12E-03	24.2	8.3	3.82E-05	53.3	29.3	1.31E-04	38.0	17.6	6.53E-04	22.2	11.1
LLM-based	SR-Scientist+PostDiag	159.1	1.60E-04	41.3	23.3	2.68E-04	30.3	17.4	9.47E-06	62.7	37.3	5.03E-05	45.4	24.1	2.81E-04	33.3	18.1

PostDiag Core Editor Prompts

System Prompt

You are a scientist working on symbolic regression. Use risk localization in the Diagnostic Risk Plane and the corresponding diagnostic evidence map to modify the symbolic structure of a candidate equation. Return a preserved or revised equation skeleton.

The input contains the current candidate, diagnostic coordinates, a diagnostic state, and local behavioral evidence. The numerical procedure has already determined the state. Follow its intervention:

- Settled (sufficiently expressive and stable)
  Prefer preserving the current skeleton. Make only small revisions clearly supported by local evidence; avoid unnecessary terms or replacing the main structure.

- Underfit (insufficiently expressive but stable)
  Preserve the main functional form. Check for missing terms, incorrect operators, or unsuitable variable combinations, and make local repairs targeting the main behavioral discrepancies.

- Overfit (sufficiently expressive but unstable)
  Prefer simplifying the skeleton. Check for redundant terms, repeated subexpressions, complex nesting, and structures that may depend on compensating coefficients.

- Mismatched (insufficiently expressive and unstable)

User Prompt Template

# Current candidate equation
{candidate_function}

# Diagnostic results
Skeleton Expressiveness Risk (Bias): {bias}
Skeleton Stability Risk (Var): {variance}
Diagnostic state: {diagnostic_state}

# Local behavioral evidence
The following evidence is computed on the same local patches of the
audit set:
- Level discrepancy: the absolute difference between the local means of
candidate predictions and observations.
- Trend discrepancy: the norm of the difference between their local
slope vectors.
- Roughness discrepancy: the absolute difference between their
deviations from their respective local linear approximations.
- Instability: local behavioral variation after bootstrap refitting of
the same skeleton. Larger values indicate greater instability.

Patch ID | Level | Trend | Roughness | Instability
{patch_evidence_rows}

# Expression constraints
Input variables: {variables}
Allowed operators and functions: {allowed_operators}
Fitted parameter slots: {parameter_slots}
Expression-size and nesting limits: {structural_constraints}
Input variable ranges: {input_ranges}

# Output requirements
Use the diagnostic state and local evidence to return one preserved or
revised equation. Represent every fitted constant as params[i] and use
only the supplied parameter slots.

Output a single Python code block with the following function signature.
The function body must contain only one return expression:
```python

def equation({signature}, params):
    return ...