Paper deep dive
RISR: Residual-Informed Scientific Equation Discovery with Large Language Models
Haobo Li, Wenshuo Zhang, Wenxiao Zhao, Eunseo Jung, Rui Sheng, Yushi Sun, Peiqin Zhuang, Hao Chen, Fenghua Ling
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 10/10/2026, 1:26:22 AM
Summary
The paper introduces RISR (Residual-Informed Scientific Equation Discovery), a method that uses large language models (LLMs) guided by residual error patterns to discover scientific equations. RISR employs a residual encoder to compress numerical error states into tokens that condition an LLM for formula proposal, and a dual-view relational encoder to predict the utility of candidate corrections for refinement. Evaluated on the LLM-SRBench, RISR achieves high in-domain and out-of-domain accuracy, outperforming baselines with the same backbone.
Entities (6)
Relation Signals (5)
RISR → evaluatedon → LLM-SRBench
confidence 95% · We evaluate RISR on scientific tasks from the LLM-SRBench.
RISR → solves → Symbolic Regression
confidence 95% · We introduce RISR, a residual-informed method that uses these error patterns to guide formula discovery...
RISR → uses → Residual Encoder
confidence 95% · A residual encoder compresses aligned inputs, targets, current predictions, and residuals into continuous tokens that condition a language model to propose formulas.
RISR → uses → Dual-View Relational Encoder
confidence 95% · For subsequent refinement, a dual-view relational encoder uses additive and regularized multiplicative residuals to predict the post-fit utility of candidate corrections.
RISR → backedby → Qwen3-Coder-30B-A3B-Instruct
confidence 90% · With Qwen3-Coder-30B-A3B-Instruct as the LLM backbone, RISR achieves 63.57% and 38.50% accuracy...
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Symbolic regression combines structural search with numerical fitting, but aggregate fit scores do not describe how the remaining error varies across inputs. We introduce RISR, a residual-informed method that uses these error patterns to guide formula discovery and learn which corrections are worth fitting. A residual encoder compresses aligned inputs, targets, current predictions, and residuals into continuous tokens that condition a language model to propose formulas. For subsequent refinement, a dual-view relational encoder uses additive and regularized multiplicative residuals to predict the post-fit utility of candidate corrections. We evaluate RISR on scientific tasks from the LLM-SRBench. RISR achieves 63.57% and 38.50% ID accuracy at the 1% and 0.1% pointwise relative-error tolerances, respectively. The corresponding OOD accuracies are 56.07% and 38.24%. RISR outperforms the reported baselines using the same backbone. The results show that our residual-informed approach can improve numerical equation recovery.
Tags
Links
- Source: https://arxiv.org/abs/2610.11387v1
- Canonical: https://arxiv.org/abs/2610.11387v1
Trouble viewing inline? Open PDF directly →
Full Text
77,830 characters extracted from source content.
Expand or collapse full text
RISR: Residual-Informed Scientific Equation Discovery with Large Language Models Haobo Li Affiliation: Shanghai AI Laboratory Wenshuo Zhang Affiliation: The Hong Kong University of Science and Technology Wenxiao Zhao Affiliation: University of California, Los Angeles Eunseo Jung Affiliation: The Hong Kong University of Science and Technology Rui Sheng Affiliation: The Hong Kong University of Science and Technology Yushi Sun Affiliation: The Hong Kong University of Science and Technology Peiqin Zhuang Affiliation: Shanghai AI Laboratory Hao Chen Affiliation: Shanghai AI Laboratory Fenghua Ling Affiliation: Shanghai AI Laboratory Abstract Symbolic regression combines structural search with numerical fitting, but aggregate fit scores do not describe how the remaining error varies across inputs. We introduce RISR, a residual-informed method that uses these error patterns to guide formula discovery and learn which corrections are worth fitting. A residual encoder compresses aligned inputs, targets, current predictions, and residuals into continuous tokens that condition a language model to propose formulas. For subsequent refinement, a dual-view relational encoder uses additive and regularized multiplicative residuals to predict the post-fit utility of candidate corrections. We evaluate RISR on scientific tasks from the LLM-SRBench. RISR achieves 63.57% and 38.50% ID accuracy at the 1% and 0.1% pointwise relative-error tolerances, respectively. The corresponding OOD accuracies are 56.07% and 38.24%. RISR outperforms the reported baselines using the same backbone. The results show that our residual-informed approach can improve numerical equation recovery. 1 Introduction Symbolic regression (SR) discovers explicit mathematical relationships from observations, providing interpretable models for scientific analysis (Schmidt and Lipson, 2009; Udrescu and Tegmark, 2020; Cranmer, 2023). These expressions make dependencies between variables explicit and allow proposed relationships to be examined and compared. Discovering an equation requires both discrete search over expression structures and continuous estimation of numerical coefficients. The two are closely coupled: the usefulness of a structural term depends on its fitted coefficients and its combination with the rest of the equation. A discovery procedure must therefore explore candidate structures while accounting for their numerical behavior on the observations. Existing methods address this coupled problem through different representations and search strategies. Sparse identification selects a small set of terms from a candidate function library (Brunton et al., 2016). Neural expression search learns a sampling policy through risk-seeking policy gradients (Petersen et al., 2021), while hybrid frameworks integrate neural search, genetic programming, and other complementary strategies (Landajuela et al., 2022). Another line of work trains expression generators on synthetic equation–observation pairs, learning to map numerical observations to symbolic expressions (Biggio et al., 2021; Kamienny et al., 2022). Together, these approaches provide mechanisms for proposing structures and organizing exploration of the expression space. Recent work increasingly uses numerical evaluation to guide large language models (LLMs) during search. TPSR incorporates fitting accuracy and expression complexity into Transformer decoding through Monte Carlo tree search (Shojaee et al., 2023). For LLM-based program generation, FunSearch combines a pretrained language model with an evaluator in an evolutionary search procedure (Romera-Paredes et al., 2024). LLM-based equation-discovery systems use scientific knowledge and executable programs to construct and revise candidate relationships (Grayeli et al., 2024; Shojaee et al., 2025a; Xia et al., 2026). LLM-SR scores fitted candidates by their mean squared error and uses a score-organized experience buffer to construct subsequent prompts (Shojaee et al., 2025a). SR-Scientist further allows the model to analyze data and residual errors through code execution (Xia et al., 2026). The form of numerical feedback determines what guidance is available for the next proposal. Aggregate fit scores rank candidate equations, but do not preserve the location or variable dependence of their errors. Equations with comparable overall losses can therefore require different corrections. Pointwise residuals expose how the remaining error varies across inputs, including nonlinear patterns, variable interactions, and localized discrepancies. Turning this information into useful search guidance requires accounting for the current equation and the effects of numerical fitting. A correction suggested by a residual pattern may produce different improvements depending on its parameterization and how it is combined with the current equation. We study how to encode these numerical states for formula generation and how to learn which corrections are worth fitting. We introduce Residual-Informed Symbolic Regression (RISR), a method for scientific equation discovery. Its search encoder compresses aligned inputs, targets, predictions, and residuals into continuous tokens that condition LLM proposals (Li and Liang, 2021). These proposals are complemented by a learned structural-edit generator. Candidate formulas are fitted and evaluated to obtain an initial equation. For subsequent refinement, a dual-view encoder captures variable-specific patterns and interactions from additive and regularized multiplicative residuals. It models relationships between variables through attention, following the principle of attention-based set encoding (Lee et al., 2019). A scoring head uses these representations to predict the utility of candidate corrections after coefficient fitting. Corrections are fitted numerically and accepted when the training loss improves, after which the residuals and action ranking are refreshed. The encoders are trained separately using synthetic data constructed from randomly transformed and composed elementary functions, independently of benchmark equations and observations. The search encoder and its projection are trained from scratch on residual-to-formula examples while the LLM remains frozen. The refinement encoder and scoring head learn from paired current equations with shared observations and candidate actions. Corrections are fitted on one subset of observations, and their outcomes on a disjoint subset provide utility supervision. We evaluate RISR on scientific problems in the released LLM-SRBench data, covering materials science, chemistry, biology, and physics (Shojaee et al., 2025b). With Qwen3-Coder-30B-A3B-Instruct as the LLM backbone, RISR achieves 63.57% and 38.50% accuracy at the 1% and 0.1% relative-error tolerances. Its performance exceeds that of existing non-LLM or LLM-based methods with the same LLM as backbone and remains competitive with substantially larger backbones. Our contributions are summarized as follows: • We introduce RISR, a residual-informed method that guides both formula search and subsequent refinement. A search encoder conditions LLM proposals on numerical error patterns, while a dual-view refinement encoder supports correction selection. • We construct complementary synthetic datasets for training RISR. Random transformations and compositions generate residual-to-formula examples, while paired current equations with shared observations and candidate actions provide fitted-utility supervision on observations held out from coefficient fitting. • We evaluate RISR on scientific problems. Our method improves overall accuracy over the reported 30B baselines and achieves the highest overall ID accuracy at the 1% tolerance among the baselines, including those using larger backbones. 2 Related work Symbolic equation discovery Symbolic regression combines expression search with numerical evaluation to recover mathematical relationships from data (Schmidt and Lipson, 2009). To make this search more efficient, physics-inspired methods such as AI Feynman exploit properties including symmetry and separability to simplify the search (Udrescu and Tegmark, 2020). PySR uses a multi-population evolutionary procedure that alternates expression evolution, simplification, and coefficient optimization (Cranmer, 2023). Moving beyond explicit search procedures, neural approaches learn from generated equations and observations: NeSymReS pretrains a Transformer to predict symbolic expressions, while end-to-end symbolic regression predicts expressions and constants jointly before numerical refinement (Biggio et al., 2021; Kamienny et al., 2022). TPSR incorporates fitting accuracy and expression complexity into Transformer decoding through Monte Carlo tree search (Shojaee et al., 2023). The score of a fitted equation does not determine which structural correction will improve it. The outcome also depends on the current residual pattern, the correction’s parameterization, and its composition with the equation. RISR addresses this decision by learning correction utility from fitted outcomes. Its dual-view encoder represents additive and regularized multiplicative correction signals, while a shared scorer uses these representations to rank corrections before fitting. LLM-based Symbolic Regression LLM-based methods combine program generation with executable evaluation. FunSearch places a pretrained LLM and an evaluator within an evolutionary program-search loop (Romera-Paredes et al., 2024). LLM-SR combines language-model proposals, parameter fitting, and evolutionary search for SR (Shojaee et al., 2025a). LaSR learns a textual concept library from high-performing hypotheses and uses it to guide subsequent search (Grayeli et al., 2024). Alternatively, SR-Scientist leverages agents for data analysis, equation evaluation, and iterative optimization, including residual inspection through code execution (Xia et al., 2026). Beyond making residuals available as agents’ input, a central question is how to represent the current numerical state for subsequent formula proposals. RISR learns this representation from synthetic residual-to-formula examples. It jointly processes aligned inputs, targets, current predictions, and residuals, compressing them into fixed-size continuous tokens that condition a frozen LLM alongside the textual search context. 3 RISR Given a set of training observations =(xi,yi)i=1nD=\(x_i,y_i)\_i=1^n, where xi∈ℝdx_i ^d contains the input variables of observation i and yi∈ℝy_i is its observed response, RISR seeks an explicit equation that predicts the response from the inputs. As shown in Figure 1, RISR proceeds in two stages. In stage A, residual-informed formula search uses a search residual encoder to supply pointwise numerical feedback to a frozen LLM through continuous tokens. LLM proposals and learned structural edits are numerically fitted and evaluated to obtain an initial equation p0p_0. In stage B, residual-guided refinement then uses a dual-view refinement encoder and utility head to rank correction actions for the current equation. Selected corrections are fitted and accepted when they improve the training loss. Figure 1: Overview of RISR. (A) Residual-informed formula search conditions a frozen LLM on textual context and numerical tokens from the search residual encoder. LLM proposals and learned structural edits are fitted and evaluated to select an initial equation p0p_0. (B) Residual-guided refinement encodes additive and regularized multiplicative residual views and predicts the fitted utility of unary and product corrections. 3.1 Residual-Informed Formula Search Inspired by recent work on LLM-guided program search and symbolic regression (Romera-Paredes et al., 2024; Grayeli et al., 2024; Shojaee et al., 2025a), we iteratively generate and evaluate equation programs. The LLM receives two complementary inputs: textual context describing previously evaluated equations and numerical tokens representing a fitted parent’s residual state. Textual context. Each prompt contains equation programs sampled from the candidate pool. The highest-scoring equation included in the prompt is the parent equation used for numerical conditioning. The textual context provides the equation-completion instruction, variable descriptions, and sampled programs. Numerical residual. The search residual encoder processes aligned tuples (xi,yi,p(xi),yi−p(xi))(x_i,y_i,p(x_i),y_i-p(x_i)). A point-attention path aggregates these observations, while a spectral-field path represents coordinate–residual interactions using frequency features and residual-region pooling. Learned gates fuse the two paths, and a projection maps their outputs to 16 continuous tokens in the LLM embedding space. These tokens accompany the textual context, conditioning executable formula proposals on the parent’s input-dependent error pattern. A learnable structural-edit generator also supplies candidates. Given the observations and current formula, the generator combines its own residual encoder with an encoding of the formula’s abstract syntax tree (AST) to predict an addition, replacement, deletion, or stop action. The resulting equation p0p_0, with the lowest mean squared error (MSE), initializes downstream refinement. Appendix A gives more details of the encoder architecture and proposal procedures. 3.2 Residual-Guided Refinement Stage B constructs residual views, scores candidate actions, and fits the selected correction. 3.2.1 Dual-View Residual Encoding For a current equation p, let pi=p(xi)p_i=p(x_i) be its prediction and ri=yi−pir_i=y_i-p_i its residual. We consider two ways to correct its predictions. Let h(x)h(x) denote a candidate correction. 1) An additive update, p+hp+h, adds an output increment, so matching the observations requires h(xi)≈rih(x_i)≈ r_i. 2) A multiplicative update, p(1+h)p(1+h), applies a relative adjustment to the current prediction. Because its output increment is pih(xi)p_ih(x_i), matching the observations requires pih(xi)≈rip_ih(x_i)≈ r_i. The two update modes therefore motivate different correction signals for the same residual state. The additive view uses rir_i directly. For the multiplicative view, direct division by pip_i would be unstable near zero, so we use the regularized signal: ui=piripi2+λ2.u_i= p_ir_ip_i^2+λ^2. (1) Here λ is 5% of the prediction root mean square (RMS), with a positive floor. The dimensionless quantity uiu_i approximates the fractional adjustment ri/pir_i/p_i when |pi||p_i| is large relative to λ, and suppresses unstable ratios near zero. Appendix B gives more details. The two residual views shown combine the same input coordinates, targets, and predictions with their respective correction signals, r and u. Each view summarizes correlations between transformed coordinates and observation signals. Pairwise statistics capture relationships between variables and the alignment of product corrections with each view’s signal. The two views share an embedding network and two relation-attention blocks that encode the variables. Pair features contribute to both attention biases and values. For each variable, we fuse its representations from the two views with a pooled global representation, producing jh_j for variable j. This set-structured design shares parameters across variables and their relationships (Zaheer et al., 2017; Lee et al., 2019). Padding and masks support inputs with up to eight variables. 3.2.2 Correction Utility Prediction Each correction action selects an additive or multiplicative update, one variable or two distinct variables, and an operator for each selected variable. We use eight operator families: Φ=z,z2,z3,sinz,tanhz,expz,log(1+z2),(1+z2)−1 = \z,\ z^2,\ z^3,\ z,\ z,\ z,\ (1+z^2),\ (1+z^2)^-1 \ (2) where z denotes a dimensionless scalar argument. A unary action applies one operator to one variable, while a product action combines two transformed variables. Masks exclude actions involving inactive variables. A shared query network combines each variable representation jh_j with an operator embedding ke_k to form a variable–operator query. For product actions, the two queries are combined symmetrically and augmented with their product correlations from both residual views. A shared scoring head combines this representation with the update mode and action type to predict the fitted utility Qθ(,p,a)Q_θ(D,p,a). Appendix B specifies the action encoding and network architecture. 3.2.3 Coefficient Fitting and Adaptive Updates Coefficient fitting. For the selected action, numerical fitting determines the coefficients of its correction function. For a subset S of observation indices, we define the relative squared-error objective LS(f)=1|S|∑i∈S(f(xi)−yisi,S)2,L_S(f)= 1|S| _i∈ S ( f(x_i)-y_is_i,S )^2, (3) where si,Ss_i,S is the target magnitude with a positive floor scaled to the target RMS on S. During refinement, fitting and acceptance use all training observations. An action with selected variables J and operators ϕkj _k_j defines the correction ha(x)=α∏j∈Jϕkj(bjzj+cj)+β,h_a(x)=α _j∈ J _k_j(b_jz_j+c_j)+β, (4) where zjz_j is the standardized value of input variable j and |J|∈1,2|J|∈\1,2\. The parameters bj,cjb_j,c_j control each operator’s input scale and shift, while α,βα,β control the correction’s amplitude and offset. The selected update mode gives either p+hap+h_a or p(1+ha)p(1+h_a). We fit the correction by variable projection: for each setting of the nonlinear parameters bj,cjb_j,c_j, weighted least squares determines the amplitude α and offset β (Golub and Pereyra, 1973; O’Leary and Rust, 2013). Adaptive action selection. At each attempt, the frozen scorer selects the highest-scoring untried action. The fitted correction is accepted if it reduces the training relative squared error by more than 0.1%0.1\%. After acceptance, the updated equation determines new residual views and a refreshed action ranking. Rejected actions leave the equation unchanged, and every attempted action is marked as tried. Refinement uses at most four action-fitting attempts, including rejected attempts, and makes no additional LLM calls. Appendix B gives the point-selection protocol, numerical fitting settings, and early-stopping criterion. 3.3 Synthetic Data and Model Training We train the search residual encoder, structural-edit generator, and refinement model separately using synthetic data generated from elementary functions with sampled transformations, coefficients, and input coordinates, independently of the benchmark equations and observations. All three components are trained offline, while the pretrained LLM parameters remain frozen. Appendix C provides more details. 3.3.1 Search residual Encoder Data and Training Target formulas combine randomly transformed elementary functions into compositional products, smooth regime transitions, or condition-gated responses. Current equations omit components or use affine approximations. Each example pairs a residual state from 256 observations over one to six variables with its target formula. Validation rejects nonfinite, nearly constant, or excessively large targets and checks numerical alignment and residual consistency. All target programs pass parsing and compilation. Audits confirm no cross-split overlap in example identities or exact normalized numerical inputs. The dataset contains 6,000 training and 768 validation examples. We train the search residual encoder and projector from scratch by minimizing the LLM loss on target-formula tokens, with the pretrained LLM parameters frozen. 3.3.2 Structural-Edit Data and Training For structural-edit generator training, we generate 20,000 training and 2,000 validation examples, each with 256 observations over one to six variables. Current formulas omit a target component, contain an extra or incorrect component, or already match the target, providing supervision for addition, deletion, replacement, and stopping. Numerical evaluation of candidate edits supplies outcome regression and ranking targets. We train the program encoder and editing modules using supervised edit losses, parameter and outcome regression, and candidate ranking. 3.3.3 Refinement-Encoder Data and Training Refinement-encoder dataset. Targets combine three to six unary or product components. For each target, we construct paired current equations with different residual patterns, covering both coarse and near-target states. Each pair shares 256 observations and eight candidate actions selected by frozen policies and random sampling. Sharing observations and actions allows supervision to capture how an action’s utility changes with the current equation. The splits contain 8,192 training pairs, 1,024 validation pairs, and three test sets of 1,024 pairs covering new observation clouds, held-out product actions, and eight-variable inputs. Together, they provide 24,576 states and 196,608 action-utility labels. Validation checks action alignment, variable eligibility, and acceptance decisions. Fitted-action utility. For each synthetic problem, we divide the 256 observations into 192 fitting rows F and 64 utility rows H. Rows in F supply the encoder inputs, coefficient fits, and acceptance decisions; rows in H measure the resulting utility. The losses LFL_F and LHL_H apply Eq. 3 to these respective subsets. For each current equation and candidate action, we fit the correction on F and accept it if it reduces LFL_F by more than 0.1%0.1\%. Let pa′p _a denote the accepted equation, or the unchanged parent when the action is rejected. Its supervision target is the clipped log reduction in error on H: U(p,a)=clip[−3,3][log10LH(p)+εLH(pa′)+ε],U(p,a)=clip_[-3,3] [ _10 L_H(p)+ L_H(p _a)+ ], (5) where ε>0 >0 stabilizes the logarithm. Rejected actions have zero utility, while accepted actions receive positive or negative utility according to their effect on H. The target therefore measures the outcome of fitting a correction beyond the rows used to estimate its coefficients. Training objective. We jointly train the relational encoder and action scorer. The loss is: ℒ=ℒvalue+0.5ℒrank+0.5ℒΔ.L=L_value+0.5L_rank+0.5L_ . (6) The value term ℒvalueL_value regresses measured utilities. The ranking term ℒrankL_rank learns their ordering within each state (Burges et al., 2005). The paired-state term ℒΔL_ matches utility differences for the same action across the two current equations. We train with AdamW (Loshchilov and Hutter, 2019) for 16 epochs, corresponding to 2,048 updates, at a learning rate of 3×10−43× 10^-4. 4 Evaluation 4.1 Evaluation Setting Dataset As many datasets may have been included in LLM pretraining corpora, LLM-SRBench (Shojaee et al., 2025b) combines known terms in the underlying equations with novel, synthetically generated terms to create problems that go beyond memorization. To ensure scientific rigor, two subject-matter experts verified all equations and visualizations of the corresponding generated data. We evaluate our approach on LSR-Synth, following Xia et al. (Xia et al., 2026). It contains 25 problems in materials science, 36 in chemistry, 24 in biology, and 44 in physics. Each problem contains a training set for SR methods, along with separate in-domain (ID) and out-of-domain (OOD) test sets. Table 1: Scientific ID accuracy (%), measured by Acc1%Acc_1\% and Acc0.1%Acc_0.1\%. The best and second-best results are highlighted in bold and underlined, respectively. Method Overall Materials Chemistry Biology Physics Acc1%Acc_1\% Acc0.1%Acc_0.1\% Acc1%Acc_1\% Acc0.1%Acc_0.1\% Acc1%Acc_1\% Acc0.1%Acc_0.1\% Acc1%Acc_1\% Acc0.1%Acc_0.1\% Acc1%Acc_1\% Acc0.1%Acc_0.1\% Without LLMs E2E 0.26 0.00 1.33 0.00 0.00 0.00 0.00 0.00 0.00 0.00 NeSymReS 3.10 0.78 8.00 4.00 2.78 0.00 4.17 0.00 0.00 0.00 DSR 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 uDSR 29.46 12.40 36.00 8.00 50.00 25.00 29.17 8.33 9.09 6.82 PySR 29.46 14.47 53.33 22.67 25.93 11.11 16.67 6.95 25.76 16.67 GLM-4.5-Air (Zeng et al., 2025) LaSR 14.21 8.53 18.67 12.00 16.67 10.18 11.11 5.56 11.36 6.82 LLM-SR 35.92 14.47 61.33 38.67 30.56 3.71 18.06 4.17 35.61 15.15 SR-Scientist 48.32 25.06 81.33 70.67 45.37 11.11 40.28 16.66 36.37 15.15 Qwen3-Coder-480B-A35B-Instruct (Yang et al., 2025) LaSR 11.89 7.49 13.33 6.67 16.67 9.26 6.95 6.95 9.85 6.82 LLM-SR 41.08 18.09 80.00 52.00 36.11 9.26 30.56 18.06 28.79 6.06 SR-Scientist 49.09 24.55 86.67 69.33 40.74 5.56 50.00 26.39 34.09 13.64 GPT-OSS-120B (Agarwal et al., 2025) LaSR 16.02 10.08 20.00 12.00 18.52 12.96 9.72 5.56 15.15 9.09 LLM-SR 28.16 11.37 64.00 38.67 22.22 3.71 11.11 2.78 21.97 6.82 SR-Scientist 63.57 49.35 74.67 60.00 81.48 64.81 66.67 43.05 40.91 34.09 Qwen3-Coder-30B-A3B-Instruct (Yang et al., 2025) LaSR 12.66 8.27 20.00 10.67 17.59 12.04 8.33 6.95 6.82 4.55 LLM-SR 24.55 7.24 29.33 6.67 30.55 10.18 15.28 5.55 21.97 6.06 SR-Scientist 32.30 16.02 81.33 52.00 22.22 5.56 22.22 8.33 18.18 8.33 SR-Scientist + RL 40.92 20.69 85.33 65.33 37.38 7.46 29.17 11.11 25.00 11.37 RISR 63.57 38.50 89.33 58.67 62.04 30.56 50.00 25.00 57.58 40.91 Evaluation Metrics We report accuracy-to-tolerance following the released SR-Scientist (Xia et al., 2026; Kamienny et al., 2022; Biggio et al., 2021). Point i passes when |y^i−yi|/|yi|≤τ| y_i-y_i|/|y_i|≤τ. Writing this pointwise indicator as bi(τ)b_i(τ), task success is Passτ(f)=[1N∑i=1Nbi(τ)≥0.95].Pass_τ(f)=1\! [ 1N _i=1^Nb_i(τ)≥ 0.95 ]. (7) Accuracy is the percentage of tasks that pass, reported at τ=0.01τ=0.01 and 0.0010.001. Figure 2: ID and OOD accuracy at a 1%1\% relative-error tolerance. Colors denote methods. Filled and hatched bars represent ID and OOD results, respectively. All LLM-based methods use Qwen3-Coder-30B-A3B-Instruct. Baselines We compare RISR against both non-LLM and LLM-based symbolic regression methods. The non-LLM baselines include E2E (Kamienny et al., 2022), NeSymReS (Biggio et al., 2021), DSR (Petersen et al., 2021), uDSR (Landajuela et al., 2022), and PySR (Cranmer, 2023), while the LLM-based baselines comprise LaSR (Grayeli et al., 2024), LLM-SR (Shojaee et al., 2025a), and SR-Scientist (Xia et al., 2026). We use the published results reported by (Xia et al., 2026). Implementation and Protocol We use Qwen3-Coder-30B-A3B-Instruct as the search backbone. Its pretrained parameters remain frozen during encoder training and evaluation. Each search has 64 outer candidate slots, including initialization. LLM proposals and learned structural edits share this budget. Stage B first fits the selected equation’s coefficients and then allows at most four new correction-fitting attempts, including rejections, without further LLM calls. Within-run equation selection and fitting use training observations, leaving the test set for evaluation. 4.2 Main Results ID Performance Table 1 reports the ID performance of RISR, averaged across results of three seeds. More details can be found in Appendix D. Among methods using the same backbone, RISR achieves 63.57% and 38.50% accuracy at the 1% and 0.1% tolerances, respectively, exceeding SR-Scientist+RL by 22.65 and 17.81 percentage points. The gains span all four domains at the 1% tolerance, with particularly large improvements in chemistry and biology: 24.66 and 20.83 percentage points, respectively. This breadth indicates that the overall improvement is shared across scientific domains rather than concentrated in a single problem category. The 30B model results are also competitive with systems built on substantially larger backbones. RISR exceeds the strongest reported Qwen3-Coder-480B result in seven of the ten dimensions. At the 1% tolerance, across all baselines, it matches the highest reported overall accuracy, 63.57%63.57\% for SR-Scientist with GPT-OSS-120B. At the stricter 0.1% tolerance, its 38.50% accuracy ranks second overall, ahead of all reported GLM-4.5-Air and Qwen3-Coder-480B configurations. These results demonstrate the effectiveness of the complete RISR pipeline at the 30B scale and show that leading accuracy at the 1% tolerance can be achieved without adopting a larger LLM backbone. As for domains, the strongest cross-backbone results occur in materials science and physics. RISR achieves the highest materials accuracy at the 1% tolerance, reaching 89.33%. In physics, it leads at both tolerances with 57.58% and 40.91%, exceeding the strongest competing results of 120B model by 16.67 and 6.82 percentage points, respectively. In chemistry, RISR ranks second at both tolerances. In biology, it ties for second place at the 1% tolerance. OOD performance RISR achieves the highest overall OOD accuracy among the compared methods at the same 30B model scale (Figure 2). At the 1% tolerance, its accuracy reaches 56.07%, compared with 39.50% for SR-Scientist + RL, a gain of 16.57 percentage points. At 1%1\%, OOD accuracy is 84.00%84.00\% in materials science, 48.15%48.15\% in chemistry, 48.61%48.61\% in biology, and 50.76%50.76\% in physics. Physics retains 40.91%40.91\% accuracy at the stricter 0.1%0.1\% tolerance. These results show that the advantage of RISR extends beyond ID fitting to prediction outside the training domain using the same discovered equations. 4.3 Ablation Study Table 2: The ablation results with one seed setting. Accuracy is reported at the 1%1\% / 0.1%0.1\% tolerances. For Stage B, all rows share the same Stage-A initial equations and numerical fitting protocol. Configuration ID accuracy ↑ OOD accuracy ↑ LLM-only formula search 7.75 / 0.00 9.30 / 3.10 + residual-informed formula search (Stage A) 36.43 / 23.26 40.31 / 30.23 + coefficient fitting 59.69 / 34.11 53.49 / 38.76 + random corrections 60.47 / 34.88 52.71 / 36.43 + learned-scored corrections (RISR) 63.57 / 37.21 55.04 / 38.76 Table 2 isolates the contributions of coefficient fitting and the learned correction-selection policy. The complete residual-informed search raises accuracy at the 1% tolerance from 7.75% to 36.43% on ID data and from 9.30% to 40.31% on OOD data compared with the LLM-only method. Residual-guided refinement provides additional gains, reaching 63.57% ID and 55.04% OOD accuracy at the 1% tolerance. Compared with random corrections, learned-scored corrections improve ID and OOD accuracy at the 1% tolerance by 3.10 and 2.33 percentage points. The results show a positive contribution from residual-informed formula search and learned action selection under the fixed correction budget. 4.4 Symbolic Accuracy Table 3: SA (%) on one seed setting. Method LLM backbone SA ↑ uDSR — 0.77 PySR — 4.65 LLM-SR GLM-4.5-Air 5.43 LLM-SR GPT-OSS-120B 4.65 LLM-SR Qwen3-Coder-480B 5.43 SR-Scientist GLM-4.5-Air 7.75 SR-Scientist GPT-OSS-120B 7.00 SR-Scientist Qwen3-Coder-480B 7.00 RISR Qwen3-Coder-30B 13.18 Symbolic accuracy (SA) evaluates whether a discovered equation can represent the target relationship. We follow the symbolic-assessment criterion of SR-Scientist (Xia et al., 2026): A candidate is counted as positive when an assignment of its adjustable constants can represent the target expression. GPT-6-Astra serves as our judge, and non-unanimous cases are reserved for human review. 17 problems receive ten Yes votes, and 112 receive ten No votes, leaving no voting disagreements. We additionally manually reviewed all 17 cases classified as positive by the LLM judge. RISR obtains an SA of 13.18%13.18\% (17/12917/129), exceeding existing methods (Table 3). SA and numerical accuracy are complementary: predictive agreement on held-out observations does not necessarily imply structural recovery under the constant-assignment criterion. Appendix E documents more details. 4.5 Case Study Figure 3: Case study on PO30. RISR retains the target components, removes an extraneous term through structural editing, and calibrates the coefficients during refinement. PO30 illustrates how RISR contributes to equation recovery (Figure 3). The target describes a nonlinear oscillator, with position x and velocity v. During Stage A, the formula proposal contains all four target components, but also includes an extra term c5∗vc_5*v. Then the learned structural-edit generator replaces this extra term with a constant and submits the edited formula for coefficient fitting and scoring. At Stage B, coefficient fitting calibrates the selected equation without changing its structure. The remaining constant becomes negligible, and the two independently fitted sinusoidal coefficients closely match the shared coefficient in the target. We examine more cases (Appendix F), including those in which numerical error is substantially reduced despite incomplete recovery of the underlying formula. These cases may inform future research aimed not only at reducing numerical error but also at recovering the scientific formulas. 5 Conclusion We present RISR, a two-stage method that uses the pointwise residual information of a current equation to guide scientific equation discovery. Residual-conditioned LLM proposals and structural-edit generators support formula search, while a dual-view encoder predicts which corrections are likely to improve the equation after fitting. For the evaluation, RISR reaches 63.57% and 38.50% ID accuracy at the 1% and 0.1% tolerances. The case study illustrates how structural editing and coefficient calibration contribute at different points in the discovery process. The evaluation results show the effectiveness of our method in the task of scientific equation discovery. AI Use Statement In this work, we used generative AI tools to assist in the writing of proofs, implement methods, assist with translation, and formulate mathematical claims. We have not used generative AI tools to generate synthetic data sets, propose or refine hypotheses, design or provide feedback on research methodology or experiments, clean and reformat datasets, support qualitative and thematic data analysis, or interpret results. Providing critical ingredients for proving mathematical claims and helping develop theoretical models or conceptual frameworks is not applicable to this work. We have reviewed all AI-assisted work. We take responsibility for the final content of this work, including text, claims, or artifacts produced with the aid of generative AI. Ethics Statement RISR is developed and evaluated using publicly available scientific benchmark data and synthetic numerical datasets generated through transformations and compositions of elementary mathematical functions. This work focuses on methodological research in scientific equation discovery. We have not identified specific ethical concerns associated with the data and experiments presented in this study. Reproducibility Statement The main paper and appendices describe the model architectures, synthetic data generation procedures, training settings, and evaluation protocols used in our experiments. To facilitate reproducibility, we will publicly release the synthetic datasets and the code for data generation, model training, equation search and refinement, and evaluation. References Agarwal et al. (2025) S. Agarwal, L. Ahmad, J. Ai, S. Altman, A. Applebaum, E. Arbus, R. K. Arora, Y. Bai, B. Baker, H. Bao, et al. Gpt-oss-120b & gpt-oss-20b model card. arXiv preprint arXiv:2508.10925. Cited by: Table 1. Biggio et al. (2021) L. Biggio, T. Bendinelli, A. Neitz, A. Lucchi, and G. Parascandolo Neural symbolic regression that scales. In International conference on machine learning, p. 936–945. Cited by: §1, §2, §4.1, §4.1. Brunton et al. (2016) S. L. Brunton, J. L. Proctor, and J. N. Kutz Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences 113 (15), p. 3932–3937. External Links: Document, Link, https://w.pnas.org/doi/pdf/10.1073/pnas.1517384113 Cited by: §1. Burges et al. (2005) C. Burges, T. Shaked, E. Renshaw, A. Lazier, M. Deeds, N. Hamilton, and G. Hullender Learning to rank using gradient descent. In Proceedings of the 22nd International Conference on Machine Learning, ICML ’05, New York, NY, USA, p. 89–96. External Links: ISBN 1595931805, Link, Document Cited by: §3.3.3. Cranmer (2023) M. Cranmer Interpretable machine learning for science with pysr and symbolicregression.jl. arXiv preprint arXiv:2305.01582. External Links: Document Cited by: §1, §2, §4.1. Golub and Pereyra (1973) G. H. Golub and V. Pereyra The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate. SIAM Journal on Numerical Analysis 10 (2), p. 413–432. External Links: Document, Link, https://doi.org/10.1137/0710036 Cited by: §3.2.3. Grayeli et al. (2024) A. Grayeli, A. Sehgal, O. Costilla-Reyes, M. Cranmer, and S. Chaudhuri Symbolic regression with a learned concept library. In Proceedings of the 38th International Conference on Neural Information Processing Systems, NIPS ’24, Red Hook, NY, USA. External Links: ISBN 9798331314385 Cited by: §1, §2, §3.1, §4.1. Kamienny et al. (2022) P. Kamienny, S. d’Ascoli, G. Lample, and F. Charton End-to-end symbolic regression with transformers. In Proceedings of the 36th International Conference on Neural Information Processing Systems, NIPS ’22, Red Hook, NY, USA. External Links: ISBN 9781713871088 Cited by: §1, §2, §4.1, §4.1. Landajuela et al. (2022) M. Landajuela, C. S. Lee, J. Yang, R. Glatt, C. P. Santiago, I. Aravena, T. Mundhenk, G. Mulcahy, and B. K. Petersen A unified framework for deep symbolic regression. In Advances in Neural Information Processing Systems, S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Eds.), Vol. 35, p. 33985–33998. External Links: Document, Link Cited by: §1, §4.1. Lee et al. (2019) J. Lee, Y. Lee, J. Kim, A. Kosiorek, S. Choi, and Y. W. Teh Set transformer: a framework for attention-based permutation-invariant neural networks. In Proceedings of the 36th International Conference on Machine Learning, K. Chaudhuri and R. Salakhutdinov (Eds.), Proceedings of Machine Learning Research, Vol. 97, p. 3744–3753. External Links: Link Cited by: §1, §3.2.1. Li and Liang (2021) X. L. Li and P. Liang Prefix-tuning: optimizing continuous prompts for generation. In Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers), C. Zong, F. Xia, W. Li, and R. Navigli (Eds.), Online, p. 4582–4597. External Links: Link, Document Cited by: §1. Loshchilov and Hutter (2019) I. Loshchilov and F. Hutter Decoupled weight decay regularization. In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019, External Links: Link Cited by: §3.3.3. O’Leary and Rust (2013) D. P. O’Leary and B. W. Rust Variable projection for nonlinear least squares problems. Comput. Optim. Appl. 54 (3), p. 579–593. External Links: ISSN 0926-6003, Link, Document Cited by: §3.2.3. Petersen et al. (2021) B. K. Petersen, M. Landajuela, T. N. Mundhenk, C. P. Santiago, S. Kim, and J. T. Kim Deep symbolic regression: recovering mathematical expressions from data via risk-seeking policy gradients. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021, External Links: Link Cited by: §1, §4.1. Romera-Paredes et al. (2024) B. Romera-Paredes, M. Barekatain, A. Novikov, M. Balog, M. P. Kumar, E. Dupont, F. J. Ruiz, J. S. Ellenberg, P. Wang, O. Fawzi, et al. Mathematical discoveries from program search with large language models. Nature 625 (7995), p. 468–475. External Links: Document Cited by: §1, §2, §3.1. Schmidt and Lipson (2009) M. Schmidt and H. Lipson Distilling free-form natural laws from experimental data. Science 324 (5923), p. 81–85. External Links: Document, Link, https://w.science.org/doi/pdf/10.1126/science.1165893 Cited by: §1, §2. Shojaee et al. (2023) P. Shojaee, K. Meidani, A. Barati Farimani, and C. Reddy Transformer-based planning for symbolic regression. In Advances in Neural Information Processing Systems, A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (Eds.), Vol. 36, p. 45907–45919. External Links: Document, Link Cited by: §1, §2. Shojaee et al. (2025a) P. Shojaee, K. Meidani, S. Gupta, A. B. Farimani, and C. K. Reddy LLM-SR: scientific equation discovery via programming with large language models. In The Thirteenth International Conference on Learning Representations, ICLR 2025, Singapore, April 24-28, 2025, External Links: Link Cited by: §A.1, §1, §2, §3.1, §4.1. Shojaee et al. (2025b) P. Shojaee, N. Nguyen, K. Meidani, A. Barati Farimani, K. D. Doan, and C. K. Reddy LLM-SRBench: a new benchmark for scientific equation discovery with large language models. In Proceedings of the 42nd International Conference on Machine Learning, A. Singh, M. Fazel, D. Hsu, S. Lacoste-Julien, F. Berkenkamp, T. Maharaj, K. Wagstaff, and J. Zhu (Eds.), Proceedings of Machine Learning Research, Vol. 267, p. 55325–55359. External Links: Link Cited by: §1, §4.1. Udrescu and Tegmark (2020) S. Udrescu and M. Tegmark AI feynman: a physics-inspired method for symbolic regression. Science Advances 6 (16), p. eaay2631. External Links: Document, Link, https://w.science.org/doi/pdf/10.1126/sciadv.aay2631 Cited by: §1, §2. Xia et al. (2026) S. Xia, Y. Sun, and P. Liu Sr-scientist: scientific equation discovery with agentic ai. In International Conference on Learning Representations, Vol. 2026, p. 75787–75811. Cited by: Appendix E, §1, §2, §4.1, §4.1, §4.1, §4.4. Yang et al. (2025) A. Yang, A. Li, B. Yang, B. Zhang, B. Hui, B. Zheng, B. Yu, C. Gao, C. Huang, C. Lv, et al. Qwen3 technical report. arXiv preprint arXiv:2505.09388. Cited by: Table 1, Table 1. Zaheer et al. (2017) M. Zaheer, S. Kottur, S. Ravanbhakhsh, B. Póczos, R. Salakhutdinov, and A. J. Smola Deep sets. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, Red Hook, NY, USA, p. 3394–3404. External Links: ISBN 9781510860964 Cited by: §3.2.1. Zeng et al. (2025) A. Zeng, X. Lv, Q. Zheng, Z. Hou, B. Chen, C. Xie, C. Wang, D. Yin, H. Zeng, J. Zhang, et al. Glm-4.5: agentic, reasoning, and coding (arc) foundation models. arXiv preprint arXiv:2508.06471. Cited by: Table 1. Appendix A Details of Residual-Informed Formula Search This section describes the configuration used for RISR. Stage A combines residual-conditioned LLM proposals and learned structural edits. A.1 Formula search and numerical inputs The search extends the evaluator–buffer–sampler design of LLM-SR (Shojaee et al., 2025a). Executable candidates are fitted to the training observations and stored with their training MSE, predictions, and residuals. The buffer maintains ten islands and samples two equation programs from one island to form each prompt. The highest-scoring equation included in that prompt supplies the parent predictions for numerical conditioning. The textual input contains an equation-completion instruction, variable descriptions, and the sampled code context. Pointwise observations are supplied through continuous embeddings rather than a natural-language account of the residuals. The residual encoder receives aligned coordinates xix_i, targets yiy_i, parent predictions pip_i, and residuals ri=yi−pir_i=y_i-p_i. It has six coordinate slots. For a problem with more than six variables, each coordinate is scored by its largest absolute correlation with the residual, considering raw, absolute-value, and squared transformations on both sides. The six highest-scoring coordinates are retained in their original order, with coordinate-index tie breaking. Routing affects only this encoder: the equation signature and numerical fitting retain all original variables. The search configuration allows up to 4,000 encoding rows. If subsampling is required, rows are sorted by residual and one aligned row is sampled from each equal-count stratum. Normalization statistics are computed on the retained rows. Coordinate j is centered by its median and divided by the larger of its interquartile range and standard deviation. A scale at most 10−810^-8 is replaced by one. Coordinates are clipped to [−12,12][-12,12] and padded to six slots. Let μy _y be the target mean and σy=maxstd(y),10−8 _y= \std(y),10^-8\. Targets and predictions are centered by μy _y and divided by σy _y. Residuals are divided by σy _y without centering. These three channels are clipped to [−20,20][-20,20]. The normalized residual’s magnitude and signed log(1+|r|) (1+|r|) complete five observation channels. Each row therefore has eleven entries, with a separate mask for active coordinate slots. A.2 Point-attention and spectral-field paths Figure 4 shows how the residual encoder maps aligned numerical observations to continuous LLM tokens through point-attention and spectral-field paths. Point-attention path. A two-layer →25611\!→\!256\!→\!256 MLP embeds each row. Two eight-head self-attention blocks mix point features, followed by layer normalization. For large observation sets, attention operates within chunks after sorting by normalized residual plus a small deterministic projection of the full row. Search-time chunks contain at most 256 points. Sixteen learned queries aggregate the valid points through two cross-attention blocks. Their shared context combines an MLP of feature means, standard deviations, and the variable mask with an MLP of the mask alone. A final layer normalization yields sixteen vectors of width 256. Spectral-field path. The second path adds an explicit feature lift to the point representations. It contains 24 coordinate channels (each coordinate alone and multiplied by the normalized residual, its magnitude, or its square), 21 coordinate-pair products multiplied by the residual, 72 sine/cosine channels, and the five observation channels. Thus, each point supplies 24+21+72+5=12224+21+72+5=122 features. The coordinate-pair channels include squares. Six positive frequencies, shared across coordinates, are learned in log space and initialized on a logarithmic grid from one to eight. Their sinusoidal arguments include a factor of π. Features associated with padded coordinates are masked. An MLP maps the lift to width 256 and adds it to the point features. These are learned frequency features, not an FFT of the observations. Seven learned residual centers, initialized uniformly from −3-3 to 33, define additional summaries. For each center, Gaussian-shaped weights depend on normalized residual distance and a shared learned positive precision. Weights sum to one over valid points. A weighted field-feature average and a learned region embedding form each summary. Sixteen field queries attend to both point features and the seven summaries. Two cross-attention/latent-self-attention pairs, followed by layer normalization, produce sixteen width-256 vectors. Fusion and projection. One learned sigmoid gate per token interpolates between the two paths. The gates are shared across examples. A two-layer GELU projector with widths →2048256\!→\!2048\!→\!2048, layer normalization, and a learned positive output scale maps the fused vectors to the Qwen embedding space. The scale starts at 0.020.02. One active marker embedding is replaced by the sixteen continuous vectors, and the sequence mask is updated. Inserted positions do not receive target-formula token labels. The LLM is frozen during encoder training, and both networks operate in evaluation mode during equation search. Figure 4: Residual encoder architecture. A.3 Auxiliary learned structural-edit generator Figure 5 shows the structural editor’s numerical and program encoders. The editor takes the training observations and the current formula. Its own residual encoder summarizes coordinates, targets, predictions, and residuals into sixteen latent vectors, projected to the editor’s hidden width. A program encoder represents the formula’s abstract syntax tree, including node types, variables, constants, parent relations, and depths. Cross-attention combines the numerical and program representations to predict an edit operation, target node, function family, variable support, and continuous parameters. The available operations are addition, replacement, deletion, and stopping. The input program is limited to 96 nodes and numerical encoding uses at most 512 observations. Predicted edits are converted into executable formulas and evaluated with lightweight training-data fitting. Their resulting residuals inform learned outcome scores used for reranking. A selected edit is then submitted to the common outer evaluator, which fits its coefficients and scores its training MSE. The editor’s predicted extrapolation-related score is learned from synthetic data. It does not evaluate the benchmark OOD split during search. All editor weights are frozen at deployment. Its pretraining and edit supervision are described in Appendix C.2. Figure 5: Numerical and program encoders within the structural edit generator. A.4 Stage-A LLM Calling LLM decoding uses temperature 0.80.8, top-k=30k=30, top-p=0.9p=0.9, repetition penalty 1.01.0, and at most 256 new tokens per completion. The buffer cluster temperature is 0.10.1. Outer coefficient fitting uses BFGS with one configured restart and a 120-second evaluation timeout. Among executable logged candidates that reproduce their recorded training fit, the equation with the lowest training MSE is selected, retaining the earlier candidate on ties. This selected equation is the Stage-A output p0p_0. Appendix B Details of Residual-Guided Refinement B.1 Dual-view construction and relational statistics Figure 6 details the dual-view refinement encoder and its relation-attention blocks. The encoder maps additive and multiplicative residual statistics to variable representations used by the correction utility scorer. For predictions pip_i and residuals ri=yi−pir_i=y_i-p_i, the regularized multiplicative coordinate solves ui=argminu∈ℝ(ri−piu)2+λ2u2=piripi2+λ2,λ=0.05maxRMS(p),10−8.u_i= _u \(r_i-p_iu)^2+λ^2u^2\= p_ir_ip_i^2+λ^2, λ=0.05 \RMS(p),10^-8\. (8) Prediction RMS is computed on the encoding rows. The positive floor makes uiu_i defined even when all predictions vanish. The additive view uses ri/σyr_i/ _y; the multiplicative view uses ui/maxstd(u),10−8u_i/ \std(u),10^-8\, without centering the signal. Both signals are clipped to [−20,20][-20,20]. Coordinates use median/IQR-or-standard-deviation normalization and clipping as in the residual encoder, but are padded to eight slots. Each view also contains target and prediction channels normalized by the target mean and standard deviation, its correction signal, the signal magnitude, and its signed logarithm. Squaring and cubing the signal yields seven observation features. For a normalized coordinate z, the fourteen coordinate transforms are Ψ(z)=[z,z2,z3,|z|,log(1+z2),e−z2,sin(πfz),cos(πfz)f∈0.5,1,2,4]. split (z)= [&z,z^2,z^3,|z|, (1+z^2),e^-z^2,\\ &\ (π fz), (π fz)\_f∈\0.5,1,2,4\ ]. split (9) Cross-correlations of the fourteen transforms with seven signals, plus their means and log standard deviations, provide 14×7+2×14+2×7=14014× 7+2× 14+2× 7=140 features per variable. Moment standardization floors the variance at 10−810^-8. For each directed variable pair (j,k)(j,k), correlations of zjzkz_jz_k, zj2zkz_j^2z_k, and zjzk2z_jz_k^2 with the correction signal, its magnitude, and its square give nine edge features. Separately, each unordered pair has 64 operator-product correlations in each view. These use the eight correction operators evaluated on normalized coordinates, standardize each operator response, and correlate their products with the view signal. Operator-response and signal standard deviations are floored at 10−410^-4; product variances are floored at 10−810^-8, and correlations are clipped to [−1,1][-1,1]. Exponentials in this statistical probe use arguments clipped to [−5,5][-5,5]. This probe clipping does not redefine the exponential operator in fitted correction formulas. A shared →128140\!→\!128\!→\!128 embedding and two four-head relation-attention blocks encode each view. Linear maps of the nine edge features contribute attention biases and values. Each block has residual connections and a width-256256 feed-forward hidden layer. Masked mean pooling of the two views supplies a global width-128128 context. Concatenating this context with each variable’s two representations gives a width-384384 vector, which a →128384\!→\!128 fusion layer maps to jh_j. Figure 6: Dual-view refinement encoder. The two views share the embedding network and relation-attention weights. The lower panel shows how pair features modify attention logits and values. B.2 Action vocabulary and shared utility head Each action specifies an additive or multiplicative update, one variable or two distinct variables, and one operator per variable. The operator order is Φ=z,z2,z3,sinz,tanhz,expz,log(1+z2),(1+z2)−1. =\z,z^2,z^3, z, z, z, (1+z^2),(1+z^2)^-1\. (10) There are 16d+128(d2)16d+128 d2 valid actions at dimension d, or 3,712 for eight variables. Product factors are stored in canonical variable order, and actions involving padded variables are masked. A shared network maps the concatenation of jh_j and a 16-dimensional operator embedding to a 64-dimensional query j,kq_j,k. The action representation is a=[L+R,L⊙R,|L−R|,m,Iproduct,ca(+),ca(×)], splitv_a= [&q_L+q_R,\;q_L _R,\;|q_L-q_R|,\\ &e_m,\;I_product,\;c_a^(+),\;c_a^(×) ], split (11) where me_m is a 16-dimensional mode embedding and the last two scalars are the action’s product correlations. Unary actions set R=0q_R=0, Iproduct=0I_product=0, and both correlations to zero. A shared →1211\!→\!128\!→\!1 head predicts fitted utility. The refinement encoder, queries, and utility head together contain 389,297 trainable parameters. B.3 Correction coefficient fitting For an observation subset S, the fitting loss is LS(f)=1|S|∑i∈S(f(xi)−yisi,S)2,si,S=max|yi|,0.01RMS(yS),ϵ64,L_S(f)= 1|S| _i∈ S ( f(x_i)-y_is_i,S )^2, s_i,S= \|y_i|,0.01RMS(y_S), _64\, (12) where ϵ64 _64 is double-precision machine epsilon. This weighted squared-error objective is different from the test-time accuracy-to-tolerance criterion. At the start of Stage B, all adjustable coefficients of p0p_0 are fitted on the full training set with its structure fixed. The trust-region reflective least-squares solver starts from the saved coefficients, uses a two-point numerical Jacobian, at most 500 function evaluations, and 10−1010^-10 for each termination tolerance. Parameter scales are max|cj|,1 \|c_j|,1\. A finite strict improvement in LSL_S is retained; otherwise the original coefficients are kept. This initial calibration is separate from the four-attempt correction budget and does not use the 0.1%0.1\% minimum-improvement threshold required for accepting a new correction. An action defines ha=αga+βh_a=α g_a+β, where ga(x)=∏j∈Jϕkj(bjzj+cj)g_a(x)= _j∈ J _k_j(b_jz_j+c_j) and zj=(xj−μj)/maxσj,10−8z_j=(x_j- _j)/ \ _j,10^-8\ uses fitting-row means and standard deviations. For fixed nonlinear parameters, the two linear design columns are (ga,1)(g_a,1) for an additive update and (pga,p)(pg_a,p) for a multiplicative update, with response y−py-p. Dividing the design and response by si,Ss_i,S gives the weighted least-squares problem for α,βα,β. The current equation’s existing coefficients remain fixed during every new correction fit. The nonlinear parameters use three deterministic starts. Unary starts for (b,c)(b,c) are (1,0)(1,0), (−1,0)(-1,0), and (2,0.5)(2,0.5); product starts for (b1,c1,b2,c2)(b_1,c_1,b_2,c_2) are (1,0,1,0)(1,0,1,0), (−1,0,1,0)(-1,0,1,0), and (2,0.5,0.5,−0.5)(2,0.5,0.5,-0.5). Slopes are bounded by [−20,20][-20,20] and shifts by [−10,10][-10,10]. Each start uses trust-region reflective least squares, a two-point Jacobian, and tolerances of 10−910^-9. Each action has a 15-second overall deadline across its starts. The best finite fit completed within the budget is retained. An attempt with no finite candidate leaves the equation unchanged; solver convergence alone is not the acceptance criterion. B.4 Encoding rows, action selection, and stopping The current refinement configuration encodes 256 training rows. When more rows are available, a deterministic sample without replacement is drawn using seed 20260831 and its indices are sorted. The same indices are used after accepted updates; their predictions, residuals, and normalization statistics are recomputed. This fixed encoding-row seed is distinct from the Stage-A search seed. Every coefficient fit and acceptance decision uses the full training set. During synthetic utility training, encoding and fitting instead use the designated 192 fitting rows. At each attempt, the frozen model ranks dimension-valid actions and selects the highest-scoring untried action, using a stable order for ties. For the fitted equation p~a p_a, acceptance requires finite training predictions and LS(p~a)<0.999LS(p)L_S( p_a)<0.999L_S(p). An accepted action updates p and its residual state; a rejection retains p. Every attempted action is marked as tried, and all attempts count toward the cap of four. Refinement stops before a new attempt when training NMSE is at most 10−1210^-12. Thus, a task may require zero new correction fits even though its initial coefficients were refitted. No Stage-B operation calls the LLM or uses test observations for acceptance. Appendix C Synthetic Data and Offline Training The two main training corpora supervise residual-to-formula generation and fitted correction utility. The auxiliary structural editor has additional residual-pretraining and edit-supervision datasets. All are generated from elementary functions, sampled transformations, coefficients, and coordinates. They do not use benchmark equations or observation arrays as generation inputs. Table 4 separates their sizes and learning roles. Table 4: Synthetic training collections. Counts are examples except for the paired utility corpus, whose counts are pairs. The three utility test sets are additional to its training and validation splits. Collection Train Validation Supervision Search residual-to-formula 6,000 768 Target formula tokens Editor residual pretraining 50,000 5,000 Families, variables, gains Structural edits 20,000 2,000 Edits and fitted outcomes Paired correction utility 8,192 1,024 Eight actions per state C.1 Search encoder data and training Three families contribute 2,000 training and 256 validation examples each. Compositional examples multiply two to six transformed factors with a global amplitude and at most ten target-program parameters; the current equation omits a factor or uses an affine approximation. Smooth-regime examples combine coordinate and condition responses through sigmoid or tanh gates; the current equation omits the late regime or uses an affine approximation. Factorized responses combine onset, hardening or softening, condition-dependent factors, and a smooth damage gate; current equations omit a condition factor or damage gate, or use an affine approximation. Each example has 256 observations and one to six variables. Compositional coordinates are sampled uniformly from randomly chosen intervals among [−3,3][-3,3], [0.05,5][0.05,5], [0,1][0,1], and [0,40][0,40]. The other families use two randomly assigned active coordinates in [0,1][0,1], with additional coordinates in [−1,1][-1,1]. Elementary transformations include shifts, rescaling, powers, saturation, exponential decay or onset, squared sine, and hinges. Generation rejects nonfinite targets, target standard deviation at most 10−410^-4, and absolute outputs exceeding 10510^5 for compositional examples or 10410^4 for the other families. Validation checks row alignment, finite inputs, and r=y−pr=y-p. All 6,768 target programs pass parsing and compilation; this is a syntactic validation and is not a proof of identifiability from observations. The residual encoder and projector are randomly initialized and trained through target-formula token cross-entropy, with all LLM weights frozen and no LoRA adapter. Training uses AdamW, one epoch over 6,000 examples, an effective batch size of eight, and 750 updates. Learning rates are 2×10−52× 10^-5 for the encoder trunk and 10−410^-4 for the projector, with weight decay 0.010.01. C.2 Structural-editor data and training The editor’s residual encoder is separately pretrained on 50,000 training and 5,000 validation examples. A sampled sum of elementary terms forms a current equation; additional sampled terms form the clean target. The library includes constants, powers, trigonometric and saturating functions, bounded exponentials, protected inverses and ratios, and variable products. Labels identify function families, active variables, and numerical fitting gains. This encoder is frozen during subsequent edit training. Edit training uses 20,000 training and 2,000 validation examples, each with 256 observations and one to six variables. Current programs omit a target component, contain an extra or incorrect component, or already match the target. These constructions provide addition, deletion, replacement, and stop labels. Up to four valid candidate slots provide measured fitting gains, gains on separately sampled wider-range synthetic inputs, and changes in expression complexity. Gaussian observation noise has standard deviation equal to 0.5%0.5\% of the clean-target standard deviation, with a positive floor. Structural-record hashes separate training and validation. The program encoder and editing modules learn supervised operation, location, family, and variable predictions, continuous parameters, outcome regressions, and candidate rankings. AdamW training uses ten epochs, batch size sixteen, learning rate 2×10−42× 10^-4, and weight decay 0.010.01. Selection uses validation loss. The synthetic wider-range observations provide offline supervision; they are not the benchmark OOD arrays. C.3 Paired-state refinement data Targets contain three to six unary or two-variable product components using the eight correction operators. Affine input slopes are sampled in [0.6,2.2][0.6,2.2] and shifts in [−1,1][-1,1]. Components are combined by addition, multiplication by one plus a component, saturation through tanh , or division by one plus a component’s square. An affine output transformation varies scale from 10−310^-3 to 10310^3. Each problem has 256 points drawn from uniform, clipped-normal, or standardized-lognormal distributions, followed by coordinate rescaling and translation. Training and validation use two to six variables; the dimensional test uses eight. Two current equations share each observation cloud and noisy target. Coarse states are formed by removing an additive component, dividing out a multiplicative component, fitting a subset of components, or combining these operations. One quarter of pairs remain coarse. For the others, each prediction is rescaled toward the clean target, p(s)=fclean+αs(pcoarse(s)−fclean),s∈0,1,p^(s)=f_clean+ _s(p^(s)_coarse-f_clean), s∈\0,1\, (13) placing its clean-error contribution, measured on the first 192 rows, into one of eight logarithmic bins spanning 10−810^-8 to 11. Observation noise is unchanged. The encoder receives the numerical state, not the target expression. Each pair shares eight dimension-valid actions. Frozen earlier residual-based policies propose top actions for both current states, and a frozen target-only policy supplies additional candidates. After deduplication, random valid actions fill the set to eight, retaining at least two random candidates. Table 5: Paired-state utility splits. Each pair has two states and eight actions per state. Only the training split contributes labels to gradient updates. Split Pairs States Utility labels Training 8,192 16,384 131,072 Validation 1,024 2,048 16,384 New observation clouds 1,024 2,048 16,384 Held-out product actions 1,024 2,048 16,384 Eight-variable inputs 1,024 2,048 16,384 Total 12,288 24,576 196,608 Splits use separate random streams and keep both states of a pair together. In the held-out-product test, operator-index pairs with (o1+o2)mod4=0(o_1+o_2) 4=0 are reserved for candidate actions and excluded from candidates in the other splits. Target-expression construction can still use all operator combinations. This tests held-out correction actions rather than disjoint target-function families. Validation checks shared-action alignment, variable eligibility, acceptance decisions, and utility arithmetic. Positive, zero, and negative labels are retained. C.4 Utility targets, losses, and optimization The first 192 rows form F and supply encoder inputs, coefficient fits, and acceptance decisions. The remaining 64 rows form H and measure fitted utility. Some target-generation normalizations use the full synthetic cloud before this assignment; H is therefore held out from correction fitting, not from every operation used to construct the synthetic problem. Neither subset contains benchmark test observations. After fitting on F, define pa′=p~a,LF(p~a)<0.999LF(p) and the fit is finite,p,otherwise.p _a= cases p_a,&L_F( p_a)<0.999L_F(p) and the fit is finite,\\ p,&otherwise. cases (14) The supervision target is U(p,a)=clip[−3,3][log10LH(p)+10−8LH(pa′)+10−8].U(p,a)=clip_[-3,3] [ _10 L_H(p)+10^-8L_H(p _a)+10^-8 ]. (15) The same loss in Eq. 12 is evaluated separately on F and H. Rejection gives zero utility; nonfinite utility-row predictions receive −3-3. An accepted action may have negative utility if it worsens prediction on H. Let U^s,a U_s,a be the predicted utility and ℓ the smooth-L1L_1 loss with transition at absolute error one. The value term averages ℓ(U^s,a−Us,a) ( U_s,a-U_s,a). For actions a,ba,b in the same state, define ds,ab=Us,a−Us,bd_s,ab=U_s,a-U_s,b, d^s,ab=U^s,a−U^s,b d_s,ab= U_s,a- U_s,b, and ws,ab=min|ds,ab|,3[|ds,ab|>0.05]w_s,ab= \|d_s,ab|,3\1[|d_s,ab|>0.05]. The ranking loss is ℒrank=∑s,a,bws,ablog(1+exp[−sign(ds,ab)d^s,ab])max1,∑s,a,bws,ab.L_rank= _s,a,bw_s,ab (1+ [-sign(d_s,ab) d_s,ab]) \1, _s,a,bw_s,ab\. (16) For each paired current state and shared action, the paired term averages ℓ((U^0,a−U^1,a)−(U0,a−U1,a)). (( U_0,a- U_1,a)-(U_0,a-U_1,a) ). (17) The total is ℒvalue+0.5ℒrank+0.5ℒΔL_value+0.5L_rank+0.5L_ . The refinement encoder, action-query network, and utility head are trained jointly from random initialization. AdamW uses learning rate 3×10−43× 10^-4, weight decay 0.010.01, and gradient-norm clipping at one. A batch contains 64 pairs, or 128 states. Sixteen epochs over 8,192 pairs give 2,048 updates. The selected model remains frozen during refinement. Appendix D Result statistics The main table and ID/OOD plot average three complete results with seeds 0, 42, and 20260831 (Table 6). We first compute the success rate for each seed, then take the mean of the three rates. Standard deviations are sample standard deviations with two degrees of freedom. Domain results use the corresponding domain’s task count. Table 6: Detailed results for the three seeds. Result set ID, 1% ID, 0.1% OOD, 1% OOD, 0.1% Seed 0 65.89 (85) 42.64 (55) 57.36 (74) 41.09 (53) Seed 42 61.24 (79) 35.66 (46) 55.81 (72) 34.88 (45) 20260831 63.57 (82) 37.21 (48) 55.04 (71) 38.76 (50) Mean ± SD 63.57±2.3363.57± 2.33 38.50±3.6638.50± 3.66 56.07±1.1856.07± 1.18 38.24±3.1338.24± 3.13 Appendix E Symbolic Accuracy: Protocol and Detailed Results SA uses the constant-assignment criterion of SR-Scientist (Xia et al., 2026): whether some assignment of the candidate’s adjustable constants makes it represent the target relationship. The original program parameters and the parameters of accepted corrections remain adjustable for this assessment. Coordinate normalization, operator choices, powers, clipping guards, and fixed literals remain unchanged. SA uses all 129 final programs and includes every accepted correction. With their saved coefficients, the exported programs reproduce stored predictions in all 387 training/ID/OOD replay checks. GPT-6-Astra evaluates each candidate–target pair in ten separate requests. The recorded requests specify high reasoning effort and a limit of 8,192 completion tokens, without an explicit temperature or sampling seed. The prompt asks about the existence of suitable constant assignments and requests a structured Yes/No answer. Unanimous Yes and No outcomes are accepted; non-unanimous outcomes are reserved for human review. The completed assessment has 1,290 valid judgments: 17 tasks receive ten Yes votes and 112 receive ten No votes, with no unresolved cases (Table 7). Table 7: Symbolic accuracy for the seed 20260831. All positive tasks receive ten Yes votes; all remaining tasks receive ten No votes. Domain Positive Tasks SA (%) Materials 1 25 4.00 Chemistry 5 36 13.89 Biology 3 24 12.50 Physics 8 44 18.18 Overall 17 129 13.18 Appendix F More Case Studies The examples below report saved formulas, selected actions, acceptance decisions, and replayed errors. The residual encoder supplies continuous numerical tokens, and the Stage-B model outputs utility scores rather than textual reasoning. F.1 CRK22: a four-step numerical refinement The chemistry reference has the structure A˙=−kA2−kAexp(−kst)+kzA21+qA4, A=-kA^2-kA (-k_st)+ k_zA^21+qA^4, (18) where t is time and A is concentration. Table 8: CRK22 action sequence. All four actions are additive and accepted by the training objective. Define rq(z)=(1+z2)−1rq(z)=(1+z^2)^-1 and lsq(z)=log(1+z2)lsq(z)= (1+z^2). Attempt Mode Correction family OOD NMSE 0 — Fitted parent 2.09172.0917 1 Additive rq(t)tanh(A)rq(t)\, (A) 4.0358×10−44.0358× 10^-4 2 Additive rq(t)lsq(A)rq(t)\,lsq(A) 2.9398×10−42.9398× 10^-4 3 Additive rq(t)sin(A)rq(t)\, (A) 1.7143×10−41.7143× 10^-4 4 Additive sin(t)rq(A) (t)\,rq(A) 6.4102×10−66.4102× 10^-6 ID NMSE changes from 1.6757×10−71.6757× 10^-7 to 4.7072×10−84.7072× 10^-8. The OOD point fraction below each of the two main tolerances changes from zero to one. Each accepted correction changes the residual state used for the next ranking. The final formula is nevertheless not counted as symbolically equivalent by the ten-vote assessment. F.2 MatSci9: additive and multiplicative corrections MatSci9 (Table 9) alternates between multiplicative and additive actions on strain ϵε and temperature T. A multiplicative action means p←p(1+ha)p← p(1+h_a), not multiplication by the displayed operator product alone. Every listed factor again has fitted affine parameters, and all four attempts are accepted. Table 9: MatSci9 refinement trajectory. Action names omit fitted affine parameters, amplitudes, and offsets. Attempt Mode Correction family OOD NMSE 0 — Fitted parent 5.478×10−75.478× 10^-7 1 Multiplicative cube(ϵ)square(T)cube(ε)square(T) 1.471×10−81.471× 10^-8 2 Multiplicative tanh(ϵ)linear(T) (ε)linear(T) 3.322×10−93.322× 10^-9 3 Additive lsq(ϵ)lsq(T)lsq(ε)lsq(T) 2.878×10−92.878× 10^-9 4 Multiplicative cube(T)cube(T) 7.566×10−97.566× 10^-9 ID NMSE decreases from 1.490×10−91.490× 10^-9 to 2.580×10−122.580× 10^-12. Final OOD NMSE is lower than the parent’s, but higher than at the third prefix. The final OOD point fraction is one at both main tolerances. This nonmonotonic test error is consistent with training-only acceptance and is retained rather than replaced by the best test-scoring prefix. Appendix G Utility-Ranking Experiments Each synthetic example contains two current equations sharing the same observation cloud and the same eight candidate actions. Let s denote one state and s¯ s its paired state. For the paired-state swap diagnostic, we use the scorer outputs produced for s¯ s to select an action for s: asswap=argmaxaU^s¯,a,Us(asswap) is then evaluated using the utility labels of s.a_s^swap= _a U_ s,a, U_s(a_s^swap) is then evaluated using the utility labels of s. The same operation is performed in the opposite direction for s¯ s. Thus, the action pool and fitting protocol are unchanged; only the predicted scores are taken from the other current equation. Table 10: Held-out utility of the selected action in fixed eight-action candidate pools. Each column contains 2,048 states. Utilities are clipped base-10 log error reductions measured after coefficient fitting. Higher values are better. Selector New observation clouds Unseen product actions Eight-variable inputs Learned utility scorer 0.46273 0.40994 0.45324 Residual-correlation ranking 0.39336 0.37877 0.39710 Uniform random selection 0.24101 0.22907 0.19557 Paired-state swap 0.27099 0.25677 0.22226 To directly test whether the learned scorer identifies promising corrections, we evaluate it on three held-out synthetic splits with a fixed pool of eight dimension-valid actions per state. Each split contains 1,024 paired problems, or 2,048 current-equation states. Every state is associated with a fixed pool of eight dimension-valid correction actions. The refinement model is frozen during this evaluation. For every state and candidate action, the correction is fitted on the designated fitting rows using the same variable-projection procedure as in the refinement stage. The fitted equation is accepted only when it reduces the fitting loss by more than 0.1%0.1\%. Its utility is then measured on the disjoint utility rows using the clipped logarithmic error reduction defined in Eq. 5. The reported score is the mean utility of the action selected by each selector. The results are shown in Table 10. The learned scorer achieves mean utilities of 0.46270.4627, 0.40990.4099, and 0.45320.4532 on new observation clouds, unseen product actions, and eight-variable inputs, respectively. The learned scorer outperforms residual-correlation ranking and uniform random selection on all three held-out splits. The paired-state swap is substantially worse, which is consistent with the scorer using information about the current equation rather than only the shared action pool. This result provides direct evidence that the learned utility scorer can identify more promising corrections than simple selection rules in controlled candidate pools. Appendix H Computational Accounting Search budget. Stage A permits at most 64 outer candidate submissions, including the initial seed equation. Invalid submissions also consume an outer slot. After the outer search terminates, executable candidates that reproduce their recorded training fit are compared using the training objective, and the candidate with the lowest training MSE is selected as p0p_0. Stage B first fits the coefficients of p0p_0 under the relative-error objective. This calibration is separate from the subsequent budget of at most four new correction-fitting attempts. Measured cost. The Stage-A configuration generates four continuations per LLM request. A complete 64-slot task search therefore uses at most 64 LLM generation requests. For the seed 20260831, in Stage B, the refinement logs contain 424 new correction attempts, of which 400 are accepted. The projected fitting objective is evaluated 461,860 times, and the recorded correction-fit durations sum to 433.98 seconds.