Physics-Informed Symbolic Regression (PISR)

• PISR is a methodology that combines symbolic regression with physical constraints (e.g., conservation laws, dimensional consistency) to yield models that are both accurate and physically plausible.
• It employs diverse approaches such as hard and soft constraints, teacher-student pipelines, and sparse library selection to control search space and improve efficiency.
• Evaluation in PISR balances data fidelity, model simplicity, and physics compliance, achieving lower error metrics and robust recovery even under noise.

Physics-Informed Symbolic Regression (PISR) denotes symbolic regression in which the search for analytical expressions is constrained or steered by physical knowledge rather than by data fit alone. In the broad literature represented here, that physical knowledge includes conservation laws, symmetries, dimensional consistency, PDE/ODE structure, boundary and asymptotic conditions, partial theoretical models, and other domain constraints, so that the recovered expression is intended to be not only numerically accurate but also physically plausible (Taskin et al., 3 Sep 2025, Bartlett et al., 17 Dec 2025).

1. Scope, definitions, and conceptual variants

PISR is not a single algorithmic recipe. It is a family of formulations in which symbolic regression is modified by physically motivated constraints, priors, or auxiliary models. In one explicit formulation, symbolic regression is the problem of finding an analytical function $f$ from data $D=\{(x_i,y_i)\}_{i=1}^n$, often with a data-fitting loss such as mean squared error; PISR augments this by imposing physical knowledge on top of the data term (Taskin et al., 3 Sep 2025). In another formulation, PISR is cast as sparse regression over a library of symbolic features, with dimensional consistency and user rules constraining which features may be generated or selected (Gajera et al., 27 Jan 2026). In yet another, PISR seeks closed-form symbolic approximations of solution functions to known differential equations, using physics-based residuals, boundary conditions, and simulation data as the search criterion (Finkelstein, 2 Oct 2025).

A useful distinction across the literature is between hard constraints, soft constraints, and teacher-student pipelines. Hard-constraint approaches reject dimensionally invalid or rule-violating expressions outright. Soft-constraint approaches add physics-related terms to the loss or objective. Teacher-student approaches first train a physics-informed neural model and then perform symbolic regression on the resulting physics-consistent surrogate.

PISR mode	Physical information	Representative papers
Hard-constrained search	Dimensional consistency, operator/feature restrictions, allowed variable combinations	(Tenachi et al., 2023, Gajera et al., 27 Jan 2026, Austel et al., 2020)
Soft-constrained optimization	Physics scores, PDE residuals, Taylor mismatch, constraint-violation objectives	(Taskin et al., 3 Sep 2025, Finkelstein, 2 Oct 2025, Gong et al., 8 Oct 2025, Ünal et al., 19 Feb 2026)
Teacher-student or gray-box pipelines	PINN, X-PINN, or UPINN produces physics-consistent targets for SR	(Kiyani et al., 2023, Majumdar et al., 2023, Podina et al., 2024)

This diversity is central to the topic. “Physics-informed” may mean dimensional analysis by construction, explicit PDE residuals, natural-language domain guidance, sparse descriptor discovery, or symbolic extraction from a physics-informed neural surrogate. A recurring misconception is that PISR must enforce the governing equations directly inside the symbolic model. Several important formulations instead encode physics through the search space, a prior over expressions, or an intermediate model (Majumdar et al., 2023).

2. Mechanisms for injecting physical knowledge

Dimensional analysis is the most common and the most explicit physical prior. In Physical Symbolic Optimization, $\Phi$-SO, each variable, constant, and target is assigned a unit vector, and dimensional rules for addition, multiplication, division, powers, and dimensionless functions are propagated through a partial expression tree. At each generation step, a mask forbids tokens whose units would violate the local requirement, so dimensional consistency is enforced by construction through an effective sampling rule of the form $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$ (Tenachi et al., 2023). A related idea appears in mixed-integer nonlinear symbolic regression, where dimensional consistency becomes a system of linear constraints on integer exponents inside a MINLP formulation (Austel et al., 2020).

A second mechanism is to constrain the search space itself. DISCOVER builds symbolic feature libraries under user-defined restrictions on operators, variable combinations, expression complexity, and dimensional consistency, and then solves a sparse regression problem

$$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$

so the physics is encoded primarily through the construction and pruning of $\boldsymbol{\Phi}$, not through an added penalty term (Gajera et al., 27 Jan 2026). A closely related sparse-library design appears in the physics-informed symbolic learner for ground-motion models, where a library $\Theta(X)$ contains physically motivated basis functions such as magnitude terms, distance saturation terms, and site terms, and Sequential Threshold Ridge Regression distills a sparse model from that library (Chen et al., 2023).

A third mechanism is to add explicit physics terms to the objective. In the LLM-based PiSR formulation, the loss is

$L = w_1 e + w_2 s + w_3 c,$

where $e$ is mean squared error, $s$ is expression-tree size, and $D=\{(x_i,y_i)\}_{i=1}^n$0 is an LLM-derived penalty computed from dimensional correctness, simplicity, and physical realism scores assigned to each candidate equation (Taskin et al., 3 Sep 2025). In the bright-soliton plasma study, the loss is decomposed as

$D=\{(x_i,y_i)\}_{i=1}^n$1

with symbolic differentiation used to evaluate PDE residuals and boundary conditions exactly on candidate symbolic expressions (Finkelstein, 2 Oct 2025). In StruSR, the objective couples a PDE residual term with a Taylor-structure mismatch term extracted from a trained PINN,

$D=\{(x_i,y_i)\}_{i=1}^n$2

so expressions are judged by both governing-equation consistency and local derivative structure (Gong et al., 8 Oct 2025).

A fourth mechanism is domain-aware priors over symbolic syntax. AI Feynman reduces the search space through dimensional analysis, symmetry tests, separability tests, and compositional decompositions, all guided by a neural surrogate that detects physical structure such as translational symmetry or multiplicative separability (Udrescu et al., 2019). More recently, domain-aware symbolic priors have been learned from 10,000 arXiv papers per domain, producing probability distributions over unary operators, binary operators, variables, depth, width, and sibling/parent relationships in expression trees; these priors are injected through KL-divergence regularization and hard masks in a tree-structured RNN (Huang et al., 12 Mar 2025). This suggests a form of PISR in which “physics” is represented not by a specific conservation law but by the statistical language of physics expressions.

Finally, dimensional analysis can act as preprocessing rather than as a runtime constraint. Buckingham-$D=\{(x_i,y_i)\}_{i=1}^n$3 groups and Ipsen’s method are used to transform algebraic equations and hidden terms in differential equations into dimensionless form before symbolic regression or UPINN-based discovery, thereby reducing input dimension and simplifying the hypothesis class (Podina et al., 2024).

3. Algorithmic families and search strategies

Genetic programming remains the dominant search mechanism in several PISR systems. DEAP, gplearn, and PySR are used as black-box symbolic optimizers whose fitness functions are modified to include LLM-derived physics scores (Taskin et al., 3 Sep 2025). SciMED layers a GA-based feature-selection stage, GA-based symbolic regression, and a Las Vegas search over full binary trees, with five explicit “knowledge injection junctions” for feature grouping, loss design, feature weights, sampling strategy, and structural constraints (Keren et al., 2022). Turbulent pipe-flow PISR modifies GPTIPS2 so that model selection is driven by a three-objective tradeoff among fit error, complexity, and physics-constraint violation derived from order-of-magnitude analysis (Ünal et al., 19 Feb 2026).

Sparse regression over symbolic libraries is the principal alternative. DISCOVER uses Orthogonal Matching Pursuit, Mixed-Integer Quadratic Programming, or Simulated Annealing to solve variants of an $D=\{(x_i,y_i)\}_{i=1}^n$4-regularized least-squares problem on a constrained symbolic feature matrix (Gajera et al., 27 Jan 2026). The ground-motion PISL framework uses Sequential Threshold Ridge Regression to select a compact subset of magnitude, distance, site, and saturation terms from a predetermined library, yielding explicit PGA and PGV equations (Chen et al., 2023).

Reinforcement learning provides another major line of development. $D=\{(x_i,y_i)\}_{i=1}^n$5-SO generates prefix expressions token by token with an LSTM, but dimensional analysis acts as an action mask and as a supervised unit-label signal, so the policy learns dimensional consistency while maximizing reward based on normalized root mean squared error (Tenachi et al., 2023). Domain-aware symbolic priors extend this idea to a tree-structured RNN, where symbol-generation probabilities are regularized toward empirical physics-domain distributions extracted from the literature (Huang et al., 12 Mar 2025).

Constraint-based global optimization supplies a more exact but less scalable alternative. The MINLP-based approach enumerates generalized expression trees whose leaves are Laurent monomials and solves, for each structure, a mixed-integer nonlinear program with dimensional constraints encoded as linear equalities on exponent vectors (Austel et al., 2020). This formulation emphasizes global optimality within a fixed structure, at the cost of rapid growth in the number of candidate trees.

Neural-symbolic hybrids form a distinct family. One variant trains a PINN to generate dense, physics-consistent data and then fits a Differentiable Program Architecture, followed by depth-first pruning; on average this pruning gives a 95.3% reduction in parameters while maintaining accuracy at par with PINNs, and on average pruning improves accuracy by 7.81% (Majumdar et al., 2023). Another gray-box variant uses X-PINNs with flux continuity to learn an unknown reaction term $D=\{(x_i,y_i)\}_{i=1}^n$6 in the Allen–Cahn equation and then applies gplearn to infer a closed-form expression from $D=\{(x_i,y_i)\}_{i=1}^n$7 pairs (Kiyani et al., 2023). A related dimensionless pipeline combines UPINNs with Ipsen’s method so that the hidden term is learned in a reduced dimensionless coordinate system before symbolic regression is applied (Podina et al., 2024). StruSR uses a trained PINN as a source of local Taylor expansions, which then guide subtree-level crossover and mutation in GP (Gong et al., 8 Oct 2025).

4. Objectives, constraints, and evaluation regimes

PISR objectives typically combine three ingredients: data fidelity, model parsimony, and physical admissibility. The exact balance varies. The LLM-integrated loss $D=\{(x_i,y_i)\}_{i=1}^n$8 is a direct weighted sum of those components (Taskin et al., 3 Sep 2025). The bright-soliton framework instead treats governing-equation residuals, boundary conditions, symmetry, and data fidelity as separate squared-residual blocks, reflecting a forward problem in which the physics constraints dominate the search (Finkelstein, 2 Oct 2025). StruSR’s hybrid objective couples PDE residuals to Taylor-coefficient alignment, and then uses masking-based attribution to preserve subtrees that are important for either physical residual reduction or structural fidelity (Gong et al., 8 Oct 2025). In turbulent pipe-flow PISR, the physics objective is not an equation residual but a violation score for envelopes on local exponents $D=\{(x_i,y_i)\}_{i=1}^n$9, derived from order-of-magnitude analysis of the Reynolds-averaged equations (Ünal et al., 19 Feb 2026).

Evaluation likewise depends on the formulation. The LLM-based PiSR study reports MAE, MSE, $\Phi$0, and an expression-tree score defined from a structural tree distance, so both numerical fit and symbolic fidelity are measured (Taskin et al., 3 Sep 2025). $\Phi$1-SO uses exact symbolic recovery, the rate of expressions with $\Phi$2, and the same rate normalized by the number of free parameters on the SRBench implementation of the Feynman benchmark (Tenachi et al., 2023). The X-PINN gray-box pipeline evaluates the learned field and hidden term by Frobenius norms over the 2D spatial grids, then uses symbolic regression to verify whether the correct functional structure $\Phi$3 is recovered (Kiyani et al., 2023). PINN-to-DPA teacher-student systems report relative $\Phi$4 errors against analytical or numerical solutions (Majumdar et al., 2023). Ground-motion PISL compares intra-event and inter-event residual statistics and their magnitude-dependent standard deviations against empirical regression models and an ANN (Chen et al., 2023).

Several representative results illustrate what “physics-informed” contributes empirically. In the LLM-based formulation, PySR on the electromagnetic-wave benchmark improves from baseline MAE $\Phi$5, MSE $\Phi$6, $\Phi$7, tree score $\Phi$8 to MAE $\Phi$9, MSE $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$0, $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$1, tree score $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$2 when combined with Mistral (Taskin et al., 3 Sep 2025). In $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$3-SO, noiseless performance ranks second in exact symbolic recovery on the Feynman benchmark, but under noise the advantage is larger: at 10% noise, when scores for all other methods drop below 20%, $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$4-SO consistently maintains a score above 50% exact recovery (Tenachi et al., 2023). In teacher-student PDE regression, pruning the symbolic architecture after PINN fitting dramatically reduces complexity while preserving accuracy (Majumdar et al., 2023). In StruSR, benchmark PDE experiments report lowest or second-lowest MAE on almost all test systems and recover expressions structurally identical to the ground truth up to minor coefficient perturbations (Gong et al., 8 Oct 2025).

These evaluation regimes underscore a central point: in PISR, accuracy alone is insufficient. Structural similarity, exact recovery, constraint satisfaction, parameter count, extrapolation behaviour, and regime-specific residuals all function as scientific criteria.

5. Application domains and representative case studies

One major application class is governing-equation discovery from algebraic or low-dimensional physical data. AI Feynman discovers all 100 equations in the Feynman symbolic regression database by combining dimensional analysis, symmetry, separability, and brute-force symbolic search after neural simplification (Udrescu et al., 2019). $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$5-SO uses dimensional analysis by construction for the same benchmark and shows particular robustness to noise (Tenachi et al., 2023). Dimensional-analysis preprocessing with Buckingham-$$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$6 groups can reduce runtime and sample requirements further; for example, on the Darcy–Weisbach equation the dimensionless AI Feynman formulation uses 10 points and 603 s versus 100 points and 3600 s in the traditional form, while maintaining a comparably small MAE (Podina et al., 2024).

A second class is descriptor discovery in materials, chemistry, and condensed-matter settings. DISCOVER is explicitly oriented toward low-dimensional descriptor formulas in computational physics, computational chemistry, and materials science, with user-configurable operator libraries, unit handling, and sparse model selection (Gajera et al., 27 Jan 2026). The review of symbolic regression in the physical sciences highlights descriptor discovery, effective theories, and compact emulators as major scientific use cases, and emphasizes that search-space design, operator selection, complexity control, and integration of theoretical information are central methodological considerations (Bartlett et al., 17 Dec 2025).

A third class is gray-box learning of missing terms in known differential equations. In the Allen–Cahn study, the PDE is written as

$$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$7

with the diffusion operator treated as known and the reaction term treated as unknown. X-PINNs with flux continuity infer $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$8 numerically from noisy or downsampled data, and symbolic regression subsequently identifies expressions close to $$\tilde{\pi}_\theta(a_t\mid s_t)\propto \pi_\theta(a_t\mid s_t)\, m_t$$9; the abstract states that the framework is able to predict the correct form and coefficients of the underlying dynamical equation when at least 50% data is used for training, while the detailed discussion further reports that about 60% of the original data should be used to recover both the correct form and reasonably accurate coefficients (Kiyani et al., 2023). In the UPINN framework, dimensional analysis and hidden-term networks recover dimensionless hidden terms such as $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$0 in an extended logistic-growth model and the correct $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$1-type structure in a rotating-bead model (Podina et al., 2024).

A fourth class is symbolic approximation of PDE solutions rather than PDE operators. PINN-generated datasets have been converted into symbolic solutions for diffusion, Kovasznay flow, Taylor–Green vortex, diffusion–reaction systems, and an industrial air-preheater, with the pruned symbolic model often achieving PINN-level accuracy while remaining interpretable (Majumdar et al., 2023). StruSR advances this line by extracting Taylor expansions from the PINN and using them to guide GP toward structurally faithful symbolic solutions (Gong et al., 8 Oct 2025). The bright-soliton plasma study goes further by directly searching the space of symbolic solution functions $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$2 and $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$3 under Maxwell–fluid residuals, boundary conditions, symmetry, and simulation data, yielding explicit closed-form approximations and a total squared-norm error of $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$4 across the enforced constraints (Finkelstein, 2 Oct 2025).

A fifth class concerns empirical correlations that remain physically constrained. In earthquake engineering, PISL constructs sparse symbolic ground-motion models whose libraries contain magnitude scaling, distance attenuation, distance saturation, site amplification, and other physically motivated terms; the resulting PGA and PGV formulas are more concise than the empirical regression and ANN baselines, while the extrapolation capability is more accurate than that of the ANN (Chen et al., 2023). In turbulent pipe flow, order-of-magnitude analysis of the Reynolds-averaged equations yields exponent envelopes that become the physics prior for symbolic regression, leading to compact friction-factor correlations validated up to $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$5 (Ünal et al., 19 Feb 2026).

These case studies show that PISR is not limited to rediscovering closed-form textbook equations. It is equally applicable to sparse descriptors, hidden constitutive terms, symbolic surrogates of PDE solutions, and semi-empirical correlations constrained by known physics.

6. Limitations, controversies, and active directions

A first limitation is that physical knowledge can be encoded too weakly or too strongly. DISCOVER states explicitly that overly restrictive constraints may exclude valid expressions, while insufficient constraints can lead to large search spaces with increased computational cost (Gajera et al., 27 Jan 2026). The same tension appears in SciMED’s scientist-in-the-loop design, where domain knowledge improves success rate and runtime, but incorrect hard constraints can misdirect the search (Keren et al., 2022). This is not a marginal issue: PISR is fundamentally a search-space engineering problem.

A second limitation is that many methods require metadata or structural assumptions that may be unavailable. $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$6-SO requires all variables and constants involved to have defined physical units (Tenachi et al., 2023). Ground-motion PISL depends on a carefully curated basis library containing the appropriate saturation and site terms (Chen et al., 2023). UPINN-based hidden-term discovery presumes a clean decomposition into known and unknown operators (Podina et al., 2024). PINN-driven symbolic extraction presumes that a sufficiently accurate physics-informed neural surrogate can be trained first (Majumdar et al., 2023).

A third limitation is computational. LLM-guided PiSR adds an LLM call to candidate evaluation and does not deeply quantify that overhead (Taskin et al., 3 Sep 2025). MINLP-based symbolic regression becomes difficult beyond modest tree depth and data size (Austel et al., 2020). PINN-first pipelines inherit the training cost and failure modes of PINNs before symbolic search even begins (Gong et al., 8 Oct 2025). Genetic-programming approaches remain sensitive to expression bloat and can require substantial evaluation budgets (Keren et al., 2022).

A fourth issue concerns what counts as “physics-informed.” Hard dimensional constraints, soft residual penalties, teacher-student pipelines, symbolic priors from literature, and prompt-based LLM critics all appear under the same umbrella, but they are not equivalent. The review literature frames this plurality directly, emphasizing that symmetries, conservation laws, asymptotic behaviour, dimensional homogeneity, and partial theoretical models can all serve as constraint sources, and that evaluation must consider interpretability, extrapolation, and physical plausibility rather than mean-squared error alone (Bartlett et al., 17 Dec 2025). A common misconception is therefore to equate PISR with only one subtype, such as PINN-coupled regression or unit-aware grammars.

A fifth issue is benchmark realism. Several prominent results are obtained on synthetic Newtonian systems, benchmark PDEs, or curated equation corpora. The LLM-based PiSR paper itself notes synthetic benchmarks, zero-shot LLM scoring, and the artificiality of prompt settings that include the ground-truth equation, which unsurprisingly yield exact recovery and should be read as an upper bound rather than as realistic discovery (Taskin et al., 3 Sep 2025). Similar caution applies whenever a physics-informed teacher supplies a near-noiseless target to the symbolic stage (Majumdar et al., 2023).

The most active directions in the cited literature are convergence of hard and soft priors, richer operator libraries, and more explicit integration with modern AI. DISCOVER identifies expanded operator libraries, improved benchmarking, and scalability enhancements as ongoing development targets (Gajera et al., 27 Jan 2026). $$\min_{\boldsymbol{\beta}} \left\| \mathbf{y} - \boldsymbol{\Phi}\boldsymbol{\beta} \right\|_2^2 \quad \text{subject to} \quad \left\|\boldsymbol{\beta}\right\|_0 \le D,$$7-SO explicitly points to adding conservation laws, symmetries, and PDE constraints to its dimensional-analysis module (Tenachi et al., 2023). The review literature emphasizes deeper integration of domain-specific knowledge, hierarchical models that share global structure across datasets, and coupling SR with LLMs or other AI systems that can formalize constraints and search drivers from qualitative theory (Bartlett et al., 17 Dec 2025). The LLM-based PiSR work is one concrete realization of that direction, replacing hand-coded physical critics with a learned evaluator trained on scientific text (Taskin et al., 3 Sep 2025).

Taken together, the literature defines PISR less as a single method than as a research program: symbolic regression shaped by physical admissibility, whether through units, constraints, sparse libraries, governing-equation residuals, Taylor structure, learned symbolic priors, or hybrid neural-symbolic workflows.