Symbolic Regression for Scientific Discovery

• SRSD is a method for deriving explicit, human-readable mathematical expressions from data that serve as interpretable scientific laws.
• It employs diverse approaches—such as genetic programming, sparse regression, and Bayesian methods—to balance data fit with complexity for robust model discovery.
• SRSD integrates domain-specific priors and constrained search strategies to ensure the resulting models are both accurate and physically plausible, aiding iterative scientific hypothesis formation.

Symbolic Regression for Scientific Discovery (SRSD) is the use of symbolic regression to recover explicit, human-readable mathematical expressions from data with the explicit aim of discovering scientific laws rather than merely achieving predictive accuracy. In this setting, the target is not simply a predictor but an interpretable equation that can function as a descriptor, a governing law, or a mechanistic hypothesis. Across the recent literature, SRSD is characterized by an optimization that balances data fit against expression complexity, by a strong emphasis on extrapolation and physical plausibility, and by a workflow in which discovered equations are inspected, validated, and iteratively refined as scientific objects rather than treated as opaque regressors (Makke et al., 2022).

1. Definition, aims, and problem formulation

At its most general, symbolic regression seeks an explicit function $f$ from data $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$ by optimizing over a hypothesis space $\mathcal{F}$ of symbolic expressions. A standard formulation minimizes a loss plus a complexity penalty,

$$\min_{f \in \mathcal{F}} \mathcal{L}(f;\mathcal{D}) + \lambda \Omega(f),$$

or, in multi-objective form, searches for a Pareto set that trades off error and complexity (Makke et al., 2022). In SRSD, the relevant solution is not simply the most accurate expression on observed data, but one whose structure is concise enough to support interpretation, critique, reuse, and extrapolation.

This scientific orientation is reflected in several specialized formulations. In sparse-regression variants, the target can be written as a linear combination of symbolic basis functions, as in

$$f(\mathbf{x}) = \sum_{j=1}^{M}\beta_j \Phi_j(\mathbf{x}),$$

with sparsity constraints such as $\|\boldsymbol{\beta}\|_0 \le D$ to induce low-dimensional descriptors (Gajera et al., 27 Jan 2026). In dynamical-systems settings, the objective may be to identify a sparse right-hand side $\dot{x}=\Theta(x)\xi$ from state observations (Makke et al., 2022). In unsupervised implicit equation discovery, the target becomes an implicit relation $f(\mathbf{x})=0$ over observed points $\mathbf{x}_n \in \mathbb{R}^D$, with no designated output variable (Yufei et al., 6 May 2025). In hierarchical Bayesian formulations, the law itself is modeled as a symbolic forest,

$$y = \beta_0 + \sum_{j=1}^{K} t_j(\mathbf{x}) \beta_j + \varepsilon,$$

where each $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$0 is a symbolic tree and posterior inference is performed over both structure and coefficients (Roy et al., 24 Sep 2025).

The field repeatedly emphasizes why this problem is scientifically distinctive. SR outputs equations in human-readable form, aligns naturally with the fact that many scientific laws are concise mathematical expressions, can extrapolate better than black-box predictors when structure matches mechanism, and supports hypothesis generation in the style of a “machine scientist” (Makke et al., 2022). At the same time, SR is computationally hard; both the survey literature and recent framework papers note its NP-hardness, making search strategy, regularization, and prior knowledge central design choices (Makke et al., 2022).

2. Expression classes and representational regimes

SRSD now spans several representational regimes, each suited to a different notion of scientific structure. Classical symbolic regression represents expressions as syntax trees whose internal nodes are operators and whose leaves are variables or constants (Makke et al., 2022). More recent work uses graphs, token sequences, or ensembles of trees; the common objective is still to recover an explicit symbolic law rather than a latent black-box map.

Regime	Core representation	Typical use
Tree-based	Syntax trees of operators and operands	Classical GP and Bayesian SR
Sparse library-based	Linear combinations of basis functions	SINDy, FFX, descriptor discovery
Sequence-based	Prefix or token sequences	Transformer and RNN symbolic regression
Forest-based	Linear combinations of symbolic trees	Hierarchical Bayesian symbolic forests
Implicit-form	Equation $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$1 on unlabeled points	Unsupervised scientific data

The tree representation remains fundamental because it matches the compositional structure of equations. In Bayesian work based on probabilistic Regular Tree Expressions, symbolic expressions are explicitly treated as trees over a ranked alphabet, with structure $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$2 separated from continuous and discrete parameters $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$3 (Schneider et al., 2023). In transformer-based models for SRSD, equations are converted from SymPy trees to pre-order token sequences, preserving a one-to-one mapping between tree and sequence (Lalande et al., 2023). In “Neural Symbolic Regression that Scales,” equations are represented in prefix notation and decomposed into a skeleton plus constant placeholders, so the model first predicts structure and then fits constants numerically (Biggio et al., 2021). In Deep Generative Symbolic Regression, the equation is likewise treated as a latent symbolic object $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$4 and SR is cast as probabilistic inference over closed-form expressions (Holt et al., 2023).

A distinct regime is descriptor-centric SR, especially in materials science. There, the learned expression is often a low-dimensional descriptor: a “fundamental materials property that is directly correlated to a desired or undesired functional property of the material” (Gajera et al., 27 Jan 2026). DISCOVER formalizes this setting as sparse linear regression over a symbolic feature library generated from primary variables and operators, making descriptor discovery a constrained subset-selection problem over candidate symbolic features (Gajera et al., 27 Jan 2026).

Another extension is parametric SR. In SEQL and HEQL, the structural form of an equation is kept fixed while coefficients vary with a parameter such as time or space, allowing discovery of parametric families such as $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$5 or PDEs with variable coefficients (Zhang et al., 2022). This is materially different from ordinary SR on a single static formula and is particularly relevant when the “law” is a structured family rather than one equation.

3. Methodological landscape

The modern SRSD literature contains several major methodological families, differing in how they represent expressions, search the hypothesis space, and incorporate scientific structure (Makke et al., 2022).

Family	Search mechanism	Representative examples
Genetic programming	Evolutionary mutation and crossover over trees	Eureqa, AFP, gplearn
Sparse regression / compressed sensing	Selection over a predefined library	SINDy, FFX, SISSO-style descriptor discovery
Bayesian symbolic regression	Posterior inference over expressions	Bayesian SR, pRTE/PTA priors, HierBOSSS
Neural / deep SR	Learned equation generation from data	EQL, DSR, NeSymReS, transformer SR
LLM- or prior-guided SR	Language-guided or prior-constrained search	DrSR, PG-SR, LLM-SR-style systems

Genetic-programming approaches remain the canonical baseline. They optimize populations of expression trees using crossover, mutation, and selection, often with explicit multi-objective pressure toward accuracy and simplicity (Makke et al., 2022). Their strengths are flexibility and direct symbolic search; their recurring limitations are computational cost, bloat, and sensitivity to hyperparameters.

Sparse-regression approaches reduce SR to regression over a large library of candidate terms. SINDy is the canonical example for dynamical systems, while descriptor-discovery systems such as SISSO and DISCOVER adopt related ideas for materials-science workflows (Makke et al., 2022). DISCOVER exemplifies a constrained, sparse linear view in which candidate symbolic features are generated by composition of operators and then selected via OMP, MIQP, or simulated annealing under dimensional and complexity constraints (Gajera et al., 27 Jan 2026).

Bayesian SR has become more structurally sophisticated. Probabilistic Regular Tree Expressions and probabilistic tree automata provide composable priors over equation trees, while HierBOSSS places a hierarchical Bayesian prior over ensembles of symbolic regression trees, enabling full posterior inference over structure and coefficients (Schneider et al., 2023, Roy et al., 24 Sep 2025). A notable distinction of these approaches is that they provide uncertainty quantification over both predictive outputs and symbolic structure.

Neural and deep SR methods use neural networks either as symbolic models or as search engines. Equation Learner architectures embed symbolic primitives directly into network layers and induce equations through sparsity (Kim et al., 2019, Zhang et al., 2022). Neural Symbolic Regression that Scales pre-trains a Transformer to predict symbolic skeletons from input-output pairs, then uses beam search and constant fitting at test time (Biggio et al., 2021). The SRSD-specific Transformer of Biggio et al. employs an encoder–decoder architecture over tabular datasets and evaluates on normalized tree-based edit distance, achieving state-of-the-art structural performance on medium and hard SRSD problems with almost instantaneous inference once trained (Lalande et al., 2023). Deep Generative Symbolic Regression further recasts SR as Bayesian inference over equations with a learned conditional generative model $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$6, and reports stronger recovery rates in higher-dimensional settings than RL-based baselines (Holt et al., 2023).

LLM-based and prior-guided systems add a different mechanism: structured reasoning about data and constraints. DrSR combines a main equation generator with a data-analysis module and an idea-extraction module, creating a closed loop in which residuals, structural patterns, and prior failures inform subsequent generations (Wang et al., 4 Jun 2025). PG-SR introduces a three-stage pipeline—warm-up, evolution, and refinement—together with executable prior constraints and a Prior Annealing Constrained Evaluation mechanism that gradually steers search toward scientifically consistent regions (Xiao et al., 13 Feb 2026).

4. Physical priors, scientific consistency, and constrained search

A defining feature of SRSD is that search is rarely left unconstrained. Domain knowledge enters through units, dimensional consistency, monotonicity, symmetry, operator restrictions, conservation laws, and executable priors (Makke et al., 2022). In scientific settings, these are not optional embellishments; they reshape the hypothesis space itself.

DISCOVER is explicit on this point. Base features are assigned physical units through pint, units are tracked through all symbolic operations, and candidate expressions are discarded immediately when dimensional rules are violated. Addition and subtraction are restricted to compatible units, multiplication and division compose units, and functions such as $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$7 and $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$8 require dimensionless arguments (Gajera et al., 27 Jan 2026). The same framework allows users to restrict variable combinations, cap expression depth, choose operator libraries, and encode shape constraints or monotonic dependence through configuration files rather than source-code edits. This yields a scientist-in-the-loop workflow in which symbolic search is guided toward “physically meaningful models” (Gajera et al., 27 Jan 2026).

The recent prior-guided literature sharpens this into a critique of empirical-risk-only SR. PG-SR formalizes the “Pseudo-Equation Trap”: expressions can achieve very low empirical risk while remaining scientifically inconsistent, especially in large symbolic hypothesis spaces and under noise or data scarcity (Xiao et al., 13 Feb 2026). Its remedy is a prior constraint checker $\mathcal{D}=\{(x_i,y_i)\}_{i=1}^N$9, where $\mathcal{F}$0 may encode causality, boundedness, monotonicity, equilibrium behavior, non-negativity, or thermo-mechanical structure, together with PACE, which anneals the influence of these priors over the search. The theoretical claim is that restricting search to the prior-consistent subspace $\mathcal{F}$1 reduces Rademacher complexity and yields tighter generalization bounds, thereby reducing the incidence of pseudo-equations (Xiao et al., 13 Feb 2026).

A complementary formalism appears in probabilistic regular tree priors. Probabilistic Regular Tree Expressions allow experts to encode prior beliefs directly over tree-structured equations, and—because regular tree languages are closed under union, intersection, difference, and complement—independent scientific constraints can be combined compositionally (Schneider et al., 2023). The paper’s Newtonian example uses separate tree priors for translational invariance and dimensional consistency, then intersects them to produce a refined prior over admissible laws (Schneider et al., 2023). This is a more structural notion of scientific prior than simple complexity penalties.

For PDE-oriented SRSD, StruSR introduces another form of prior: locally structured physical information extracted from a trained Physics-Informed Neural Network. It computes local Taylor expansions of the PINN solution, compares Taylor coefficients of candidate symbolic expressions to those of the PINN, and combines Taylor mismatch with physics residuals in a hybrid fitness function $\mathcal{F}$2 (Gong et al., 8 Oct 2025). A masking-based attribution mechanism then identifies which subtrees contribute most to structural alignment and residual reduction, so mutation and crossover preserve physically significant substructures (Gong et al., 8 Oct 2025).

Taken together, these methods reflect a shift from generic equation fitting toward scientific consistency as a first-class objective. A plausible implication is that “symbolic regression” in SRSD increasingly denotes constrained scientific model discovery rather than unrestricted search over formulas.

5. Benchmarks, metrics, and evaluation controversies

SRSD evaluation is unusually sensitive to how “success” is defined. Predictive accuracy alone is insufficient when the goal is structural recovery of a scientific law. The survey literature accordingly distinguishes predictive error, expression complexity, robustness, and recovery of ground-truth structure as separate metrics (Makke et al., 2022).

A major benchmark in this area is the SRSD dataset introduced by Matsubara et al. and used by subsequent work. It consists of 120 equations derived from Feynman-style laws, partitioned into 30 easy, 40 medium, and 50 hard problems, with realistic sampling and standardized train/validation/test splits (Lalande et al., 2023). For this benchmark, the central metric is the normalized tree-based edit distance

$\mathcal{F}$3

where $\mathcal{F}$4 is the Zhang–Shasha tree edit distance and $\mathcal{F}$5 is the number of nodes in the ground-truth tree (Lalande et al., 2023). This metric explicitly evaluates structural similarity rather than merely output fit. On this benchmark, the transformer model of Biggio et al. reported the lowest normalized tree edit distance on the medium and hard subsets, while maintaining almost instantaneous inference once trained (Lalande et al., 2023).

Another evaluation issue is algebraic equivalence. Deep Generative Symbolic Regression argues that token-level cross-entropy is fundamentally misaligned with SR because mathematically equivalent expressions can have different token sequences; it instead trains against an NMSE-derived reward so that equivalent formulas are not penalized for syntactic differences (Holt et al., 2023). This addresses a longstanding tension in SRSD between canonical expression recovery and functional equivalence.

That tension is explicit in “Fast Symbolic Regression Benchmarking.” Martinek argues that existing rediscovery benchmarks overemphasize recovery of “the one” expression form and rely too heavily on SymPy-based equivalence checks, while continuing search even after a correct expression has been found (Martinek, 20 Aug 2025). Using curated lists of acceptable expressions and callback-based early termination on the SRSD benchmark, the paper reports that the rediscovery rate of SymbolicRegression.jl increases from 26.7% to 44.7%, with 41.2% less computational expense; TiSR reaches 69.4% rediscovery with 63% time savings (Martinek, 20 Aug 2025). This establishes that benchmark design can materially change reported SRSD performance.

The evaluation literature therefore contains an active controversy: whether success should mean exact symbolic recovery, structural proximity, functional equivalence, or scientifically acceptable equivalence classes. The recent trend is toward explicitly structural metrics, equivalence-aware evaluation, and early stopping once scientific adequacy is reached.

6. Scientific applications and software ecosystems

SRSD has been applied across physics, chemistry, materials science, biology, engineering, and environmental modeling (Makke et al., 2022). In each case, the scientific value lies in discovering a law or descriptor that is both predictive and interpretable.

In physics and dynamical systems, landmark examples include the recovery of conserved quantities and governing equations from motion data, as well as sparse identification of nonlinear systems such as the Lorenz equations (Makke et al., 2022). StruSR targets PDE solution discovery for advection, diffusion, Poisson, heat, and wave equations, using PINN-derived Taylor structure and PDE residuals to recover symbolic solutions with strong structural fidelity (Gong et al., 8 Oct 2025). PIE extends this frontier to unsupervised implicit discovery, addressing settings where data are sampled from an implicit manifold $\mathcal{F}$6, as in particle trajectories or orbital constraints, and reports clear gains over GP- and RL-based baselines on synthetic and AI-Feynman-derived implicit datasets (Yufei et al., 6 May 2025).

In materials science, descriptor discovery is a central application. DISCOVER is explicitly designed for descriptor and equation discovery in materials science, physics, and chemistry, with examples situated around crystal structure stability, ion mobility, and magnetic materials (Gajera et al., 27 Jan 2026). The framework targets low-dimensional, physically meaningful descriptors generated from atomic or thermodynamic variables and selected by sparse regression under units and operator constraints (Gajera et al., 27 Jan 2026). HierBOSSS likewise demonstrates SRSD on a single-atom catalysis dataset, identifying compact descriptors involving heat of sublimation, oxidation energy, and vacancy energy that align with physically motivated formulations (Roy et al., 24 Sep 2025).

Biology, chemistry, and engineering provide cases where priors and extrapolation are particularly important. PG-SR evaluates on E. coli growth, material stress–strain data for Al 6061-T651, chemical reaction kinetics, and nonlinear oscillators, reporting strong in-distribution and out-of-distribution NMSE while remaining robust to noisy data, data scarcity, and varying prior quality (Xiao et al., 13 Feb 2026). DrSR is likewise evaluated across oscillators, E. coli growth, stress–strain, transformed physics equations, and synthetic chemistry, and reports improved valid equation rate, accuracy, OOD generalization, and search efficiency relative to both classical and earlier LLM-based baselines (Wang et al., 4 Jun 2025). In environmental science, EQL-based symbolic regression has been used for short-term wind speed forecasting, yielding closed-form equations that use only a few selected features and support interpretation of spatial-temporal dependencies across Danish cities (Abdellaoui et al., 2021).

The software ecosystem increasingly reflects practical scientific workflows. DISCOVER emphasizes Python-native integration, configuration-based search constraints, and optional GPU acceleration with CUDA and Metal Performance Shaders, positioning SRSD as a reproducible component in data-processing and screening pipelines rather than an isolated executable (Gajera et al., 27 Jan 2026). This suggests a shift from stand-alone symbolic discovery programs toward modular scientific ML infrastructure.

7. Open problems and research directions

The current literature presents SRSD as a mature but rapidly evolving area rather than a solved problem (Makke et al., 2022). Several limitations recur across methods.

Scalability remains fundamental. Search spaces grow combinatorially with variable count, operator breadth, and expression depth; even methods designed for scale typically rely on strong priors, synthetic pretraining, or constrained grammars (Makke et al., 2022). DISCOVER explicitly notes that it is not optimized for fully unconstrained SISSO-level massive feature spaces and is instead intended for guided, constrained exploration (Gajera et al., 27 Jan 2026). PIE currently assumes implicit laws in at most three variables and a fixed operator set, illustrating that even strong pretraining-based methods remain domain-bounded (Yufei et al., 6 May 2025).

Noise, false discovery, and validation are persistent concerns. The survey literature highlights spurious correlations, overfitting, and reproducibility as open challenges, with different algorithms and hyperparameters often yielding different equations on the same data (Makke et al., 2022). PG-SR reframes this as a scientific-consistency problem, arguing that empirical risk minimization alone is insufficient and that explicit priors are needed to avoid pseudo-equations (Xiao et al., 13 Feb 2026). A plausible implication is that future SRSD systems will increasingly resemble constrained scientific modeling platforms rather than generic regressors.

Uncertainty quantification is still unevenly developed. Bayesian systems such as HierBOSSS offer full posterior inference and structural uncertainty, but many high-performing neural and LLM-based systems still output point estimates or ranked candidates without calibrated posterior mass (Roy et al., 24 Sep 2025). The survey literature identifies principled uncertainty over both equations and constants as a major future need (Makke et al., 2022).

Another frontier is expanding the scope of discovery. Parametric equation learning extends SR from single formulas to structured families (Zhang et al., 2022). StruSR suggests a pathway from neural PDE surrogates to symbolic PDE solution discovery (Gong et al., 8 Oct 2025). PIE opens unsupervised implicit SR for scientific data lacking input-output labels (Yufei et al., 6 May 2025). DISCOVER points toward richer operator libraries, including integrals, differential operators, and special functions, while PG-SR and DrSR suggest deeper integration of executable priors, residual analysis, and reflective reasoning (Gajera et al., 27 Jan 2026, Xiao et al., 13 Feb 2026, Wang et al., 4 Jun 2025).

Finally, benchmarking itself remains unsettled. Structural metrics, acceptable-form equivalence classes, out-of-distribution evaluation, and domain-specific testbeds are all still being refined (Lalande et al., 2023, Martinek, 20 Aug 2025). As a result, SRSD is increasingly evaluated not only by rediscovery of canonical laws, but by whether it yields parsimonious, physically plausible, and experimentally actionable equations under realistic data conditions.