Moving from table to graph in physics-informed spatio-temporal symbolic regression

The pursuit of mathematical models that capture the governing principles of natural and engineered systems is a cornerstone of scientific inquiry^1–5. Across disciplines such as physics⁶, biology⁷, and economics⁸, researchers aim to identify concise and interpretable formulations that describe observed dynamics. To be exact, scientific inquiry, from as early as Isaac Newton, proceeds through a three-step cycle of observation, hypothesis generation, and hypothesis validation^9,10. Traditionally, researchers begin by collecting data about the world (observation), which they then use to develop explanatory hypotheses (hypothesis generation). A good hypothesis is required to allow for extrapolation and the prediction of new data within the observed system during the hypothesis validation stage^11–13. Modern science, in general, and modern exact science, in particular, is based on an equation-based formalization of hypothesis^14,15. This process, more often than not, requires a comprehensive understanding of a field, deep familiarity with the collected observations, and a decent amount of creativity and mathematical expertise^16,17.

Recently, with the rapid increase in both the volume and complexity of available data, there is a growing need for automated methods that can discover such models directly from observations^18–23. To this end, scholars developed a wide range of symbolic regression (SR) models, which allow the discovery of analytical expressions that fit given data^24–28. Unlike black-box machine learning (ML) approaches on the one hand, and simplistic statistical models (such as linear regression), which are explainable but more often than not lacking expressiveness, on the other hand, SR produces interpretable symbolic forms that can be readily analyzed and understood by scientists²⁹.

One field in which SR has gained much attention is physics^30–32. The inherent goal of physics is to uncover governing equations that describe the laws of nature, making SR a natural candidate for accelerating this process^33–35. Multiple studies have demonstrated the ability of SR to rediscover known physical laws such as Newton’s equations of motion, conservation laws, and simple partial differential equations directly from data^36,37.

Nonetheless, existing SR models suffer from two limitations that hinder their practical applicability in physics-driven domains. First, SR models assume that data are presented in a tabular form, where each sample is treated independently, without accounting for the spatio-temporal dependencies intrinsic to physical systems. Recent SR models, especially those using Recurrent Neural Network (RNN) based models, partially tackle this limitation as they do consider the temporal context of samples^38–42. Other works use large language models as generalized and implicit models to support the SR search process^43–45. In reality, physical processes evolve dynamically in time and space, with local interactions shaping the behavior of the global system. Ignoring these structures leads to models that may fit data but fail to capture the underlying dynamics. Second, SR searches for algebraic expressions mapping variables but does not explicitly account for the governing equations being ordinary or partial differential equations (ODEs/PDEs). Since most physical systems are defined through such equations, neglecting their structure often results in models that lack physical meaning, generalizability, and predictive power.

To this end, in this study, we propose a simple yet powerful data representation approach with SR fitting procedure that aims to tackle these two limitations regardless of the SR method used under the hood. Namely, the proposed method leverages a dual representation of the raw data - one using the traditional tabular representation and the additional graph-based data representation, where nodes encode spatio-temporal states and edges capture spatial or temporal dependencies. This representation allows us to move beyond the traditional tabular paradigm and to directly model the interconnected nature of physical systems. Through a series of benchmarks on synthetic and real-world datasets, we demonstrate that the proposed method is capable of improving the performance of out-of-the-box SR models in recovering governing equations with high fidelity, even under noisy conditions. Hence, the main contributions of this work are:

• We propose a dual tabular–graph representation for spatio-temporal physical data, where nodes represent space–time observations and edges encode spatial/temporal dependencies, enabling SR to leverage physical locality.
• We introduce a solver-agnostic SR fitting objective that combines consistency in the raw feature space with consistency in a physics-aware latent space learned by a GNN, thereby improving equation recovery without modifying the underlying SR solver.
• We provide a comprehensive evaluation across ten SR models and ten physics-derived benchmark systems, including robustness to observational noise and data-efficiency analyses.

Figure 1 presents the schematic view of the advantage of our dual-representation approach.

Related work

SR is a modeling paradigm aimed at uncovering explicit mathematical expressions that capture the relationship between input variables and an output variable^46–48. Unlike traditional regression approaches such as linear regression, which require a predefined functional form⁴⁹, SR explores a vast space of candidate expressions to identify the one that best fits the observed data⁵⁰. This flexibility makes SR especially valuable in scenarios where the underlying functional form is unknown or highly nonlinear yet required to be expressed via a human-readable equation^51–53.

Roughly, SR methodologies can be grouped into four categories, depending on the computational strategy employed: brute force, sparse regression, deep learning, and genetic algorithms⁵⁴. Brute-force methods exhaustively enumerate candidate equations, theoretically guaranteeing a solution. However, in practice, they are computationally intractable and prone to severe overfitting⁵⁵. Their main distinction lies in how they partition and traverse the search space⁵⁶. In contrast, sparse regression approaches reduce the search space by exploiting optimization techniques that promote parsimony. A representative example is SINDy, which uses Lasso regression to identify compact models of nonlinear dynamical systems from time series data⁵⁷. SINDy alternates between least-squares fitting and thresholding steps to enforce sparsity⁵⁸. Another family of approaches is based on deep learning. These methods leverage neural networks to generate symbolic expressions, benefiting from robustness against noise and outliers, but often struggle with generalization beyond the fitting domain⁵⁹. For example, the Deep Symbolic Regression (DSR) framework⁶⁰ employs reinforcement learning to train a recurrent neural network (RNN) that generates symbolic formulas, using a risk-seeking policy gradient variant to improve expression discovery. Finally, genetic algorithm–based methods treat candidate equations as individuals in an evolving population⁶¹. Through selection, crossover, and mutation, the population gradually refines toward expressions that better explain the data. This evolutionary strategy supports interpretability and flexibility, since it does not assume a fixed model class^62,63. The Python package gplearn, for instance, implements this approach and has been shown to rediscover known physical laws⁶⁴. In gplearn, equations are represented as trees whose substructures evolve via stochastic optimization guided by fitness evaluation. In recent years, several SR toolkits have been proposed, making the methodology increasingly accessible. Widely adopted libraries include DEAP (Distributed Evolutionary Algorithm in Python)^65,66, gplearn^67,68, and PySR^69,70, all of which are based on genetic programming principles.

The application of SR in physics is especially appealing, as scientific domains often involve multivariate, noisy experimental data from nonlinear systems, yet demand interpretable analytical expressions to uncover governing laws^30–32. However, purely data-driven SR is insufficient in this context, since candidate equations must respect constraints such as dimensional consistency, conservation principles, and symmetries. To address these challenges, Physics-Informed Symbolic Regression (PiSR) has been proposed^32,71. PiSR integrates prior physical knowledge directly into the search process, guiding model discovery by embedding principles such as conservation laws, scale invariance, or sparsity assumptions³².

Several notable approaches exemplify this trend. AI-Feynman combines physics-inspired heuristics with machine learning and deep learning methods, demonstrating strong performance across diverse physics problems⁷². Similarly, the Equation Learner (EQL) framework⁷³ uses a specialized feedforward neural network designed to capture dynamical equations and extrapolate beyond observed data. Its architecture integrates linear mappings, nonlinear unary and binary operators, and a linear readout layer, enabling efficient gradient-based fitting. Another contribution is SPRINT (Scalable Pruning for Rapid Identification of Null vecTors), a sparse regression method that generalizes exhaustive search via iterative singular value decomposition (SVD) and pruning, enabling efficient equation discovery without prior model assumptions⁷⁴.

While PISR is tailored to enforce symbolic structure and physical constraints, Physics-Informed Neural Networks (PINNs) represent a parallel and highly active vein of research in scientific machine learning^75–77. PINNs embed known ordinary/partial differential equations (ODEs/PDEs) into the loss function of neural networks, forcing the learned approximation to satisfy the governing equations (via automatic differentiation) while fitting observed data⁷⁸. Recent reviews highlight that PINNs have become a central pillar in data-driven scientific modeling and inverse problems^79,80. They often combine boundary/initial condition terms, residual losses from PDE operators, and data mismatch terms into a unified fitting objective⁷⁸.

Recently, several enhancements have been proposed to improve the flexibility, robustness, and interpretability of PINNs. Variants in network architecture, domain decomposition, adaptive collocation sampling, adaptive loss weighting, and hybrid symbolic–neural formulations have all been introduced to mitigate fitting difficulties (e.g. stiffness, imbalance of loss terms) and improve generalization⁸¹. For example, GPT-PINN frames a meta-learning mechanism, where a meta-network composes pre-trained PINNs as basis units to generalize across parameterized PDEs efficiently⁸². TGPT-PINN extends this to nonlinear model reduction by combining transformed layers and shock-capturing losses⁸³. Another line of work, GINN-KAN, seeks to build interpretability into the PINN architecture itself, combining ideas from interpretable networks with classical Kolmogorov–Arnold (KAN) representations so that the resulting PINN is less of a black box⁸⁴. Despite their promising performance, they suffer from key drawbacks. First, fitting is often unstable due to imbalanced loss terms, and scalability to high-dimensional or complex systems is limited by computational cost. Second, PINNs also struggle with sharp features such as shocks or turbulence. Third, they remain largely black-box models without interpretable outputs, and finally, they are sensitive to noisy or sparse boundary conditions. These limitations motivate hybrid approaches with PISR, which can improve interpretability and efficiency while retaining physical consistency.

Data representation and symbolic regression fitting procedure

The motivation behind the proposed dual representation arises from the fundamental structure of physical reality itself. In physics, systems evolve according to local interactions in space and time — quantities such as velocity, temperature, or concentration at a given point are influenced primarily by their immediate spatial and temporal neighbors. This locality principle underpins the mathematical form of differential equations, where derivatives capture rates of change over infinitesimal distances or times. The tabular representation, while convenient for machine learning, inherently breaks these spatial and temporal linkages by treating each observation as an independent sample. In contrast, a graph representation naturally encodes these dependencies: nodes represent local states, and edges express their physical couplings. By embedding this graph-based structure alongside the traditional tabular view, we hypothesize that the SR process can infer differential relationships more efficiently and therefore, effectively reconstruct the latent operators (e.g., gradients, fluxes, interactions) that govern the system. Thus, the dual representation bridges two complementary perspectives: the table captures measurable data fidelity, while the graph encapsulates physical continuity and causality, together forming a more faithful computational mirror of how nature operates.

Formally, the proposed method aims to allow SR models to better discover the underlying governing equations of dynamical systems directly from spatio-temporal data by integrating Graph Neural Networks (GNNs) with SR. The methodology is organized into three principal stages: (i) construction of a spatio-temporal grid graph, (ii) graph-based feature representation learning using GNN, and (iii) SR optimization with a dual representation-based loss function. This architecture enables the recovery of interpretable, physically consistent equations from complex, high-dimensional data. A schematic overview of the in**li ne-eq-IEq1SRworkflow is illustrated in Fig. 2.

Graph representation construction

We assume an experiment in which measurements of one or more physical quantities are collected over time. These quantities may describe the state of a single object (e.g., the position and velocity of a falling ball) or the evolution of a distributed field (e.g., temperature, pressure, or velocity over a fluid domain). In the former case, spatial information may be fixed or limited, whereas in the latter, measurements vary continuously across both space and time. Regardless of the experimental setup, the collected data are reorganized into a unified structure that treats every observation as a node embedded in a four-dimensional space-time manifold ((x, y, z, t)), associated with a set of measured features. This graph-based formulation allows the learning model to capture two key forms of dependency:

• Temporal continuity, which ensures that consecutive measurements influence one another in accordance with physical causality.
• Spatial coherence, which enforces locality and smoothness of interactions among neighboring spatial points.

As a result, the learned representation respects the geometry of the physical process and provides a structured substrate upon which SR can operate to infer governing equations. To this end, let the dataset of measurements be denoted by i n line- eq- IEq2 where in l ine- eq- IEq 3 are the spatial coordinates of the i-th observation, in l i^ne-eq-IEq4 is the associated time coordinate, and in l ine-eq -IEq5 is the vector of measured or derived physical quantities.

From inline-eq-IEq6, we construct an undirected graph i n lin e-eq-IEq7, where each node in l ine-eq-IEq8 represents a unique space-time observation and is defined as i_n l ine_-e q_-I E_q9 and the edge set inline-eq-IEq10 is composed of spatial and temporal links i n li n e-eq-IEq11 such that spatial edges (in l ine-e q-I E q12 ) connect nodes corresponding to nearby spatial coordinates at the same time step, while temporal edges (in l ine-e q-I E q1 3 ) connect nodes corresponding to the same spatial location at consecutive times, where inline-eq-IEq14 and inline-eq-IEq15 are distance thresholds that determine the radius of spatial and temporal connectivity, respectively. This construction yields a four-dimensional spatio-temporal graph in which the geometry of the physical process is preserved and accessible to graph-based learning models.

Technically, this construction assumes that each observation can be associated with a meaningful spatio-temporal index, either through explicit spatial coordinates and timestamps inli ne-eq-IEq16, or through a consistent sensor/mesh indexing that can be mapped to inli ne-eq-IEq17. In consequence, the adjacency is derived from a user-defined neighborhood rule (e.g., radius or kNN) via the thresholds inl ine-eq-IEq18 . When the topology changes over time but coordinates remain available (e.g., particle systems or moving sensors), the same pipeline applies by recomputing inline-eq-IEq19 and inline-eq-IEq20 per time step (or within sliding windows), yielding a dynamic graph sequence. Thus, in such settings, dynamic-graph GNN encoders can be used as drop-in replacements for inline-eq-IEq21^85,86. To this end, in problems where spatial relations are unknown or latent, an additional graph-structure-learning component is required to infer the adjacency from data⁸⁷. For example, interaction graphs can be inferred directly from trajectories via latent-relation models⁸⁸.

Importantly, this graph-based structure explicitly defines “local differential relationships” among the components of in both space and time. For instance, spatial edges approximate partial derivatives ⁱⁿ _{line-eq-IEq23}, ⁱⁿ _{line-eq-IEq24} , and ⁱⁿ _{line-eq-IEq25} through finite differences between neighboring nodes at constant t, while temporal edges approximate temporal derivatives ⁱⁿ _{line-eq-IEq26} between successive time steps at fixed spatial coordinates. Hence, the graph implicitly encodes the discretized structure of a partial differential operator acting on : disp l ^ay -e q - (1) Equ1 where inline-eq-IEq28 denotes the latent differential operator governing the physical process, and inl ine-eq-IEq29 are coefficients to be discovered. An example of this phase, where a mass is falling under the influence of gravity, is provided in the Appendix.

Latent space learning

Once the spatio-temporal graph (G) is constructed, it is provided as input to the GNN component. The GNN operates directly on this structured graph representation, transforming node-level physical measurements into a set of learned latent embeddings that encode both spatial and temporal dependencies. These embeddings serve as physics-aware features that the SR module later uses to infer explicit analytical expressions of the underlying dynamics. Technically, each node in l ine-eq-IEq30 corresponds to a unique space-time observation inli ne-eq-IEq31, with an associated feature vector in l i^ne-eq-IEq32 representing the observed or derived physical quantities (e.g., position, velocity, acceleration, temperature, pressure). The initial node embedding is therefore i^nli e -eq-I E q33 n where inline-eq-IEq34 and inline-eq-IEq35 are learnable input projection parameters that map the raw feature space into a latent dimension inline-eq-IEq36. In addition, each edge inli ne-eq-IEq37 carries attributes inline-eq-IEq38 that may include spatial distance inl i ne-eq-IEq39 , temporal interval in l ine - eq-IEq40, or learned positional encoding. In addition, each edge inli ne-eq-IEq41 is associated with attributes inline-eq-IEq42 that may include spatial distance inl i ne-eq-IEq43 temporal interval in l ine - eq-IEq44, or learned positional encoding. In the general formulation, such attributes may be incorporated into the message function i^{nline-eq-IEq45}. In our implementation, however, they are used to construct and weight the adjacency matrix, as described in Graph embedding implementation Section.

The GNN updates node embeddings through iterative message-passing operations⁸⁹. At each layer l, the node inline-eq-IEq46 aggregates messages from its neighborhood i nli n e-eq-IEq47 according to:

d^isp _la y -_eq-Eq_u2 where i^{nline-eq-IEq48} is a message function, i^{nline-eq-IEq49} is a node update function, and inline-eq-IEq50 is a permutation-invariant aggregation operator (e.g., sum, mean, or attention-weighted average). In practice, i^{nline-eq-IEq51} and i^{nline-eq-IEq52} are implemented as multi-layer perceptrons (MLPs) or attention kernels parameterized by i^{nline-eq-IEq53}, yielding a differentiable mapping in l in^e-eq -I E q⁵⁴ . The formulation above describes a general message-passing framework. In this work, we instantiate it using a Graph Convolutional Network (GCN) where the GCN corresponds to a special case in which messages are aggregated via normalized adjacency weights without explicit edge-conditioned transformations.

Notably, after L layers of message passing, each node inline-eq-IEq55 obtains a latent embedding i^nli that encodes not ne-eq-IEq56 only its instantaneous physical state but also the aggregated influence of its spatio-temporal neighborhood: d^isp a y^-eq qu1 0 _{l -} E

This embedding space differs fundamentally from the raw measurement space. In the raw data domain, nearby points may appear similar due to sensor noise or limited resolution, obscuring meaningful dynamical variations. In contrast, the GNN-generated latent space emphasizes consistent dynamical relationships: message passing integrates directional information over temporal edges and spatial gradients, effectively filtering out incoherent noise while amplifying physically consistent signals. Consequently, nodes corresponding to regions or time steps governed by the same physical law form coherent clusters in the latent space, while those associated with different regimes or boundary conditions become separable. This property allows the learned representation to be interpreted as a physics-aware manifold inline-eq-IEq57, where distances reflect similarity in underlying dynamics rather than raw measurements. Within this manifold, trajectories that share identical governing equations follow geometrically smooth paths, whereas perturbations in physical behavior correspond to measurable divergences in latent geometry.

Symbolic regression fitting procedure

Let us assume K independent experiments indexed by i n lin e - e q -IEq58, each producing a spatio-temporal dataset i^{nline-eq-IEq59} and its graph i^nli n e- ^eq-I E^q60 . In addition, let i^nli n e-^eq- E_q61 denote the node features I (observed/derived physical quantities) on i^{nline-eq-IEq62}. A trained GNN encoder inline-eq-IEq63 maps node features on a graph to latent embeddings i^nli n e-eq_-IEq 6₄ . An SR candidate is an expression i n lin e-eq-IEq65 with structure s (an expression tree over a grammar inline-eq-IEq66 of primitives inl in e- eq -IEq 67 , etc.) and constants i n l^ine-eq-IEq68. The expression defines a right-hand side (RHS) operator inline-eq-IEq69 acting on the state inline-eq-IEq70 (scalar or vector field), e.g., d_is p l_a y-e q-E q_u11 Given experiment-specific initial/boundary conditions i^{nline-eq-IEq71} and exogenous inputs i^{nline-eq-IEq72}, we define a numerical time-stepping simulator inline-eq-IEq73 (e.g., explicit/implicit Runge–Kutta for ODEs, finite-difference/volume/element for PDEs) that yields a synthetic dataset on the same grid/graph: i^spl a y- d eq^-Equ 1² We assume inline-eq-IEq74 produces outputs on the same node set i ^{nline-eq-IEq75} (or a known interpolation/projection aligns nodes). The training objective is defined over the set of candidate expressions i n lin e-eq-IEq76, where each expression is evaluated through the simulator inline-eq-IEq77 and the GNN encoder inline-eq-IEq78. For a given candidate, the simulator produces ^{line-eq-IEq80}, which are then passed through the GNN to synthetic trajectories i obtain latent embeddings in n^{line-eq-IEq79} and corresponding node features in ^lin e -eq-^IEq8 1 . The quality of the candidate is quantified by the following components.

The first term penalizes discrepancies between simulated and observed data in the measurable (raw) feature space. Optionally, finite-difference derivatives can be included to enforce local smoothness and temporal or spatial consistency:

e d_isplay - q u (3) -Eq 3 lin where inl ine-eq-IEq82 are feature-wise weights, and i^nli denote numerical derivatives derived from i^{nline-eq-IEq85} and _ne-eq-IEq83 , in_e-eq-IEq84 ^{line-eq-IEq86}, respectively. in Beyond the raw space, we require that the simulated trajectories produce GNN embeddings consistent with those extracted from real data. This ensures that the candidate equation reproduces the encoded physical dependencies captured by inline-eq-IEq87:

q d^ispl _ay- eq- E u (4) 4 where MMD stands for Maximum mean discrepancy, and the total latent-space loss is expressed as d_isplay-eq - E q^u5 (5) This dual-level consistency enforces that the learned equation preserves both explicit signal patterns and implicit physics-aware relations. To favor concise, interpretable, and physically meaningful solutions, two regularization families are applied: model complexity and physics-based constraints:

display - eq-Equ 6 where inline-eq-IEq88 denotes the symbolic tree length, inline-eq-IEq89 encodes penalties for symmetry or invariance violations (e.g., conservation or Galilean invariance), and inline-eq-IEq90 represents residuals on boundary/initial conditions and other physics constraints (e.g., conservation or multi-field coupling residuals) when available. Beyond data fidelity, we explicitly encourage discovered equations to satisfy known physical structure. In Eq. (6), we use two complementary mechanisms: (i) an invariance/symmetry penalty inline-eq-IEq91 and (ii) a physics-residual term inline-eq-IEq92 that measures violations of hard constraints such as boundary/initial conditions and conserved quantities. Combining all components yields the total objective function:

d isp l ay-eq-Equ7 where the coefficients in l ine-eq-IEq93 balance the trade-offs among reconstruction accuracy, latent consistency, model parsimony, and physical validity. Minimizing inline-eq-IEq94 guides the discovered expression toward simultaneously fitting experimental observations and aligning with the manifold of physically consistent dynamics learned by the GNN.

Importantly, while the dual tabular-graph losses enforce data and representation fidelity, the above equations also allow injecting domain knowledge to prevent physically implausible equations. In particular, inline-eq-IEq95 can encode invariances by penalizing violations of known transformation properties (e.g., invariance to time translation, spatial translation/rotation when applicable, or Galilean transformations in dynamical systems). In a similar manner, inline-eq-IEq96 can be defined to enforce boundary and initial conditions by measuring the mismatch of the simulated trajectory induced by E at known boundary/initial samples (or their discrete approximations). These constraint terms reduce the effective search space of SR and act as a physically meaningful regularizer, which is especially beneficial in nonlinear and high-dimensional settings where many algebraically valid but physically spurious expressions can fit the data.

Computational complexity

We analyze the additional computational overhead introduced by the proposed dual tabular-graph representation. Formally, let i nl in e-eq-IEq97 be the number of space–time nodes, i nl ine-eq-IEq98 the number of edges (spatial+temporal), M the number of raw features per node, inline-eq-IEq99 the latent embedding dimension, and L the number of GNN message-passing layers.

Given the edge definitions based on spatial and temporal neighborhoods (thresholding via inl ine-eq-IEq100 ), the resulting graph has i n ^{line-eq-I E} edges, where ⁱnli ⁿe-eq-IEq102 are the average spatial/temporal degrees. A naive q101 all-pairs construction is inl ⁱne-eq-IEq103, but using standard spatial indexing (e.g., grid hashing or kNN structures) yields inl ine - e q-IEq104 time and inl i ne-eq-IEq105 space to store the adjacency. Moreover, each message-passing layer aggregates neighborhood information. With sparse adjacency, the dominant per-layer cost is inl ine-eq-IEq106 for aggregation added to inl iⁿ _e-eq-IEq107 for applying learned linear maps, resulting in:

In addition, for a candidate expression E, evaluation requires running the simulator inline-eq-IEq108 to generate trajectories on the same node set, and then encoding them through the frozen GNN encoder. Thus, per candidate:

where inline-eq-IEq109 depends on the numerical scheme and grid resolution. Compared to a tabular-only SR objective (which primarily incurs inline- e q-IEq110), the proposed method adds an explicit latent-consistency term whose overhead scales with sparse message passing.

Experiments

To systematically evaluate the effectiveness of the proposed dual tabular–graph data representation, we conducted a comprehensive experimental study involving 10 SR models and 10 synthetic datasets. Each dataset was generated from a distinct known physical equation, representing diverse dynamical behaviors and mathematical structures. The evaluation was organized into three experiments. In the first, we compared the performance of all SR models before and after applying the proposed dual representation to assess its overall impact on equation recovery accuracy. In the second, we tested the resilience of each SR method to increasing levels of observational noise, examining the robustness of the discovered expressions under data perturbations. Finally, in the third experiment, we evaluated data efficiency by progressively reducing the amount of available training data, measuring how the proposed approach improves generalization and recovery fidelity in data-scarce regimes. Figure 3 provides a schematic view of the experimental design.

Baseline models

Table 1 summarizes the ten SR models evaluated in this study. The selected methods represent a broad spectrum of algorithmic paradigms, including genetic programming, deep learning, sparse optimization, Bayesian inference, and physics-informed neural modeling. Specifically, evolutionary approaches such as gplearn and PySR serve as strong baselines for stochastic tree-based optimization, while sparse regression methods like SINDy and FLEXPDE-SR capture the benefits of parsimonious model identification under physically structured feature libraries. Neural and hybrid architectures, including SciMED, EQL, DSR, and AI-Feynman, provide a complementary perspective, leveraging differentiable symbolic operators and deep policy-based exploration.

SR model	Paradigm	Description	Source
SciMED	Multi-Stage Hybrid (Symbolic + GNN)	Integrates machine learning, symbolic search, and physics priors to recover interpretable governing equations from spatio-temporal data.	90
AI-Feynman	Hybrid (Physics-Informed Deep Learning)	Combines neural networks, dimensional analysis, and symbolic simplification to identify analytical physical laws from data.	72
Deep Symbolic Regression (DSR)	Deep Reinforcement Learning	Uses a recurrent neural network trained via policy gradients to generate symbolic expressions in a data-driven manner.	91
PySR	Genetic Programming (Gradient-Assisted)	Employs an evolutionary search over symbolic trees with gradient-based mutations to efficiently find compact equations.	92
gplearn	Genetic Programming	Classical SR implementation based on stochastic expression-tree evolution with selection, mutation, and crossover operations.	64
SINDy	Sparse Optimization	Identifies governing equations using sparse regression over a library of candidate nonlinear terms with Lasso regularization.	93
EQL (Equation Learner)	Neural–Symbolic Network	Neural network architecture with symbolic activation functions designed to learn analytic relations and extrapolate beyond training data.	94
GP-SR	Bayesian / Monte Carlo	Integrates symbolic expression structures into Gaussian Process kernels to perform probabilistic equation discovery.	95
Symbolic PINN (S-PINN)	Physics-Informed Neural Network	Embeds differential equation constraints in a neural network while maintaining explicit symbolic form for discovered terms.	96
FLEXPDE-SR	PDE-Constrained Sparse Regression	Discovers PDE forms by combining finite-difference approximations with sparse regression to recover differential operators.	97

Finally, physics-informed and probabilistic formulations, represented by Symbolic PINN and GP-SR, integrate prior knowledge or uncertainty quantification into the discovery process.

Datasets

Table 2 presents the set of ten benchmark equations used to generate the synthetic datasets in this study. The selected equations span a diverse range of physical domains and mathematical structures. In particular, we picked two examples from five equation types: functional, ODE, PDE, Integral, and ODE with delay. Together, these datasets form a balanced and interpretable benchmark suite that reflects the core mathematical structures encountered in physical modeling and provides a robust foundation for evaluating SR performance.

For every equation listed in Table 2, synthetic datasets were generated through controlled numerical simulation, producing both the features (independent variables) and targets (dependent variables) required for SR fitting. Namely, each dataset was constructed to include informative variables directly involved in the underlying physical law as well as non-informative distractors, allowing for the assessment of each SR model’s ability to identify relevant relationships under realistic conditions. The exact feature list for each equation is provided in the Appendix.

For each equation, a set of input features corresponding to the true physical variables was generated according to the mathematical form of the governing equation. For instance, in the case of Hooke’s Law, the displacement inline-eq-IEq121 and spring’s coefficient inline-eq-IEq122 were sampled from a uniform distribution within a predefined physical range, while the corresponding force inline-eq-IEq123 was computed analytically as i n line-eq-IEq124. In addition to the correct physical variables, each dataset was augmented with three extraneous features that were statistically independent of the governing equation. These features were generated deterministically using smooth mathematical functions of time or space to mimic structured yet irrelevant signals. Specifically, for each dataset we created three auxiliary variables defined as inlin e -e q-IEq12 5 where the coefficients inl ine -eq-IEq126 were sampled uniformly within predefined ranges (in l ine-e q-IEq127, in l ine-e q-IEq128, in l ine -eq-IEq129). For PDE-based datasets, the same formulation was extended to include both spatial and temporal dependencies, such that inlin e-eq-IEq130 combined sinusoidal or polynomial variations over the simulation grid. These extraneous features were included in all datasets but were not part of the underlying governing equations. Their role was to emulate structured measurement redundancy and confounding effects that frequently arise in physical experiments, thereby allowing the evaluation of each SR model’s capacity to isolate the truly relevant variables from correlated yet non-causal signals.

The samples were produced consistently across all experiments. For each equation, multiple independent runs (indexed by i n line-eq-IEq13 1 ) were generated, each representing a separate realization of the system with slightly perturbed initial conditions and parameters. Within each experiment, data were collected over a finite time horizon or spatial domain discretized into in l ine-eq-IEq132 time steps and in l ine-eq-IEq133 spatial points, yielding spatio-temporal tuples of the form inli ne- eq-IEq134. This procedure ensured both intra-experiment temporal consistency and inter-experiment variability, supporting generalization and cross-validation analyses of the SR models.

To introduce stochasticity and realism, two types of noise were added. First, a dedicated noise feature inline-eq-IEq135 was appended to every dataset, consisting of a random variable drawn independently from a zero-mean normal distribution i nli ne-eq-IEq136. This feature, while unrelated to the physical process, was included as an explicit distractor to test the robustness of the SR models to irrelevant information. Second, controlled Gaussian noise was injected directly into both the source (input) and target (output) variables of the governing equations. Specifically, for a target variable inline-eq-IEq137 and an input feature inline-eq-IEq138, noisy versions were generated as in l i n e-e q-I E q1 3 9 where inline-eq-IEq140 is a user-defined noise level determining the signal-to-noise ratio of the experiment. By default, we set iⁿ _{line-eq-IEq141} to be one percent of the average value of inline-eq-IEq142 and inline-eq-IEq143 respectively. Importantly, the datasets are generated once and used across all SR models to ensure consistency.

Experiments

Three complementary experiments were designed to rigorously evaluate the effect of the proposed dual tabular– graph data representation on SR performance. Each experiment targeted a distinct aspect of model quality— accuracy, robustness, and data efficiency, as these are the main properties required from a useful physically-relevant SR^108–110.

The first experiment assesses the baseline and improved performance of each SR method before and after the application of the proposed dual representation. Namely, for every SR model listed in Table 1, training and inference were performed twice: once using the standard tabular data representation and once using the combined tabular–graph representation. The comparison was based on three metrics: mean absolute error (MAE), coefficient of determination, and expression tree agreement⁴⁵. In addition, we report a post-hoc physical consistency check for each discovered equation by measuring its symmetry/invariance violation and boundary/initial residual in the Appendix. The goal of this experiment was to quantify the extent to which the proposed representation enhances the SR performance in reconstructing physical equations in a realistic while data-rich configuration.

The second experiment evaluates the robustness of the SR models to observational noise. For each dataset, Gaussian noise of increasing magnitude was independently added to both the input (source) and output (target) variables, following the formulation described in Section Datasets. The noise standard deviation inline-eq-IEq144 was systematically varied across six levels in the range in l ine -eq-IEq145 relative to the signal amplitude. Each SR model was retrained under each noise condition using both representations, and the resulting degradation in equation recovery accuracy was recorded for the same three metrics. This experiment quantified the stability and reliability of the symbolic discovery process under measurement uncertainty, as well as the degree to which the proposed dual representation mitigates overfitting to noisy observations.

The third experiment examines the performance of the SR models when trained on increasingly limited subsets of the available data. For each dataset, the total number of training samples was progressively reduced from inline-eq-IEq146 to inline-eq-IEq147 of the original set in ten linearly spaced intervals. At each data availability level, all SR models were trained and evaluated under both the standard and dual representations for all three metrics. This experiment aimed to measure the ability of each SR model—and particularly the proposed representation—to infer correct governing laws under data-scarce conditions, reflecting practical scientific scenarios where data collection is costly or constrained.

_q -IE

Graph embedding implementation

In all experiments, we employed a Graph Convolutional Network (GCN) architecture¹¹¹ as the encoder, chosen for its balance between interpretability, computational efficiency, and expressive capacity for relational learning. While more expressive edge-conditioned or attention-based GNNs could be employed, we adopt a GCN in this work to balance interpretability, efficiency, and stability. The GNN encoder consisted of two graph convolutional layers followed by a nonlinear activation (ReLU) and a final linear projection to obtain the latent representation i^nli. Specifically, the encoder inline-eq-IEq149 was defined as i^nli n e-eq^-IEq 1⁵⁰ where ne-eq-IEq148 i^{nline-eq-IEq151} is the normalized adjacency matrix of the graph, and inl i2 are trainable weight matrices. The output embeddings i^nli n e-^eq-E_q15 3 provide a physics-aware latent representation that reflects both the local I interactions and global spatio-temporal structure of the system. This model was trained in a self-supervised manner using the ground-truth simulated data prior to the SR fitting stage. The objective was to reconstruct local physical relationships between neighboring nodes and to preserve temporal consistency. To this end, the GNN was optimized with a composite loss function: i_nli n e_- e q15 4 ^{line-eq-IEq155} denotes reconstructed node features from the decoder, inline-eq-IEq156 denotes temporal finite differences, and inline-eq-IEq157, inline-eq-IEq158 are weighting coefficients controlling reconstruction and smoothness terms, respectively. This pretraining where in stage ensured that the GNN learned to embed local differential structure and temporal evolution without any explicit symbolic supervision. Once trained, the encoder parameters inline-eq-IEq159 were frozen and reused for all SR models and experiments to ensure a consistent and fair comparison.

This formulation corresponds to a specific instance of the general message-passing scheme in Eq. (2), where messages are defined as weighted neighbor embeddings and aggregated via the normalized adjacency matrix. Edge attributes (e.g., spatial and temporal proximity) are incorporated implicitly through the construction of the adjacency matrix i^{nline-eq-IEq160}. Specifically, edges are defined based on spatio-temporal neighborhoods, and their corresponding weights reflect these relations prior to normalization. Thus, edge information influences message passing through i^{nline-eq-IEq161}, rather than via explicit edge-conditioned transformations. Importantly, the decoder is implemented as a multi-layer perceptron (MLP) in l i^ne - e^q-IEq162, applied ^lin I^Eq1 - eq . The decoder parameters inline-eq-IEq164 are trained jointly with the encoder independently to each node in _e - ₆ 3 during the self-supervised pretraining phase. After training, the decoder is discarded, and only the encoder inline-eq-IEq165 is retained and frozen for the symbolic regression stage.

Results

Table 3 presents the comparative results of the ten SR models evaluated under two configurations: the standard tabular representation and the proposed dual tabular–graph representation. Each metric (MAE, i^{nline-eq-IEq166}, and ETD) was averaged across all ten benchmark datasets, with the reported uncertainty corresponding to the standard deviation across runs. As shown, the dual-representation configuration consistently improves or maintains performance across all models and metrics, with relative gains generally ranging between 1 and 3%. From the perspective of numerical accuracy (MAE), Symbolic PINN achieves the lowest mean error, followed closely by SciMED and DSR, while gplearn exhibits the highest error among the evaluated methods. In terms of explanatory power (i^{nline-eq-IEq167}), EQL performs best, with GP-SR and SciMED closely trailing, while DSR and AI-Feynman are less consistent. For symbolic fidelity (ETD), which measures the structural similarity of the discovered expressions to the ground-truth equations, SciMED and GP-SR outperform all others, followed closely by Symbolic PINN. Notably, even though the absolute numerical improvements appear modest, they are systematic across models and metrics, indicating that the dual representation improves convergence stability and promotes the discovery of physically consistent symbolic forms.

Figure 4 presents the robustness of the SR models to noise in both input and output variables. For each of the ten physical systems, increasing levels of zero-mean noise with standard deviation in l ine -eq-IEq173 were applied uniformly to all features. The figure is divided into the three metrics: MAE, i^{nline-eq-IEq174}, and ETD, averaged across all datasets. As expected, the numerical accuracy (MAE) and symbolic fidelity (ETD) deteriorate as the noise amplitude increases, while the explanatory power (i^{nline-eq-IEq175}) systematically declines. However, the rate of degradation differs substantially between models: deep-learning–based and physics-informed approaches (e.g., SciMED, Symbolic PINN, and FLEXPDE-SR) exhibit higher resilience to perturbations, while purely evolutionary or sparse-regression methods (e.g., gplearn, PySR, and SINDy) are more sensitive to noise. Across all metrics, the proposed dual tabular–graph representation consistently mitigates performance loss, maintaining up to 10–15% smaller error growth rates compared with the baseline tabular setting. Similar analysis for noise robustness for three more types of noise (heavy-tailed additive noise, impulsive outlier noise, and systematic bias/drift noise) is provided in the Appendix.

Figure 5 quantifies data efficiency (i.e., how performance scales with the number of independent experiments available for symbolic discovery). The number of experiments, K, was varied from 200 to 1000 while keeping all other settings fixed. The results show a clear sublinear trend: while early additions of new experiments substantially improve accuracy and stability, marginal gains diminish beyond approximately inline-eq-IEq176. The proposed dual-representation approach accelerates convergence toward the asymptotic regime, achieving equivalent accuracy and structural consistency with roughly 30–40% fewer experiments than the baseline models. Notably, the gains are most pronounced for methods that rely on strong inductive biases (e.g., SciMED, FLEXPDE-SR, and EQL), suggesting that the graph-informed encoding provides complementary contextual information that promotes generalizable equation recovery.

Figure 6 shows a two-dimensional performance landscape over the noise intensity inline-eq-IEq177 and the number of experiments K. The top row reports heatmaps of the aggregated metric values (mean across benchmarks and SR solvers) for MAE, i^{nline-eq-IEq178}, and ETD under the dual tabular-graph representation and the bottom row reports the corresponding delta heatmaps relative to the tabular baseline (MAE as relative reduction in %, and i^{nline-eq-IEq179}/ETD as absolute improvements), thereby visualizing how performance changes jointly as the input/output perturbation level increases and as the available experimental budget grows. Across the grid, the delta maps highlight larger improvements in the more challenging regimes (higher inline-eq-IEq180 and smaller K), while remaining consistently positive as K increases and the landscape saturates.

To further examine the smoothing effect attributed to the GNN encoder, we analyzed the learned latent representations under both clean and noisy conditions. Specifically, we compared low-dimensional projections of the raw feature space and the GNN latent space after training. Table 4 summarizes the geometric properties of the raw feature space and the learned GNN latent space under both clean and noisy settings using neighborhood compactness, computed as the average distance to the k nearest neighbors, and cluster separation, measured using a standard clustering-quality criterion such as the silhouette score. Notably, a smoothing effect induced by the GNN is indicated by both metrics when the latent space exhibits lower compactness values and higher separation values than the raw space, especially under noisy observations.

Discussion

In this study, we introduced a dual data representation framework for SR that integrates both tabular and graph-based encodings of spatio-temporal data. The proposed approach aims to overcome two major limitations of conventional SR methods - their assumption of independent samples and their lack of explicit physical structure awareness. Namely, by coupling standard feature-based representations with GNN embeddings that capture spatial and temporal dependencies, the method enables SR models to infer governing equations that are both numerically accurate and physically interpretable.

In order to evaluate the proposed method, we conducted three experiments exploring the performance comparison, noise resilience, and data efficiency, using ten physics-derived benchmark equations and ten representative SR algorithms spanning multiple methodological paradigms. Across all three experiments (Table 3 and Figs. 4, 5 and 6), the results consistently demonstrate that the proposed dual tabular–graph representation improves the robustness, data efficiency, and symbolic accuracy of SR models, on average. Specifically, Table 3) outlines that both genetic programming and sparse regression approaches, such as PySR and SINDy, exhibited the largest relative gains, with up to 25–40% improvement in symbolic match scores and 15–30% lower MAE. These improvements suggest that embedding relational priors through GNN-based encoding substantially guides the symbolic search process, reducing overfitting and improving physical consistency. Deep learning– based and hybrid frameworks (e.g., DSR, EQL, and GP-SR) displayed moderate but consistent gains, typically between 5–10%, indicating that they already internalize some spatio-temporal dependencies implicitly through their architectures. Notably, physics-informed models such as Symbolic PINN, SciMED, and FLEXPDE-SR achieved the lowest overall error levels, showing that coupling symbolic discovery with physical constraints remains a highly effective paradigm, consistent with findings by^72,82,84. In the noise resilience experiment, captured by Fig. 4, the dual-representation models maintained stable performance under Gaussian perturbations up to i_n l ine-eq-IEq184, whereas tabular-only baselines showed substantial degradation, especially among evolutionary and sparse methods. The relative robustness of dual models—evident from 10–15% slower deterioration in MAE and ETD—aligns with prior observations in physics-informed learning, where structural or inductive priors enhance generalization under uncertainty^78,79. Finally, the data-efficiency experiment, presented in Figs. 5 and 6, further revealed that dual-representation models required only 30–40% of the training data used by their tabular counterparts to reach comparable accuracy and symbolic fidelity. This efficiency gain was particularly pronounced for SciMED, FLEXPDE-SR, and EQL, which achieved near-asymptotic accuracy at fewer than six experiments, suggesting that the relational encoding effectively captures shared structure across spatially correlated observations. These findings mirror recent results in graph-based and operator-learning frameworks^112,113, which emphasize that leveraging domain structure dramatically improves sample efficiency in physics-driven inference. To this end, physical constraints improve interpretability not only by filtering out unphysical candidates, but also by guiding SR toward compact expressions whose terms admit direct physical meaning. In practice, invariance and boundary/initial-condition penalties discourage equations that rely on accidental correlations or dataset-specific artifacts, which become more prevalent as nonlinearity and dimensionality grow. As a result, the discovered symbolic trees tend to be shorter and more stable across noise/data-scarcity regimes, and the retained operators are more likely to correspond to physically meaningful mechanisms (e.g., transport, diffusion, restoring forces), rather than brittle high-order combinations of features. From a practical standpoint, even a lightweight graph representation, constructed from spatial or temporal adjacency, can significantly enhance symbolic model discovery without altering the underlying SR algorithm. Thus, the proposed dual loss function can be easily integrated into existing SR frameworks as an auxiliary constraint, providing an accessible pathway for upgrading current SR pipelines to physics-informed variants. Furthermore, the observed improved robustness to both noise and data sparsity makes this method particularly suitable for experimental settings where measurements are limited or uncertain, and indicates that the proposed method allows SR models to capture the physical dynamics even if its signal is slightly corrupted.

This study is not without limitations. First, the current implementation relies on a fixed GNN encoder trained separately from the SR optimization; joint training could further enhance representation alignment but would increase computational cost. Second, the graph construction procedure assumed a known spatial and temporal topology, which may not be available in all experimental settings. More precisely, our current implementation presumes that the neighborhood structure can be derived from coordinates (or known indices) using fixed rules such as radius or KNN connectivity. If the neighborhood evolves over time but coordinates remain observable, a straightforward extension is to rebuild the edge sets inl ine-eq-IEq185 per time step (or within windows) and employ a dynamic-graph GNN encoder. If the topology is latent (i.e., relations are not directly observable), then graph structure learning is required to infer the adjacency from data. Moreover, enforcing physical consistency requires specifying the relevant constraint family (e.g., conservation/symmetry/coupling constraints), which is straightforward in our controlled benchmarks but may require domain expertise in complex multi-physics settings. Third, the benchmark datasets were synthetic and noise-controlled, and thus may not fully reflect the complexity and irregularities of real-world measurements. As such, future work should evaluate the proposed methods’ capabilities in discovering new equations from real-world experimental data. Fourth, for high-dimensional problems, SR may face greater challenges to produce decent results. The proposed method is not inherently helpful for this use case. Future work may aim to combine dimensionality reduction methods as a pre- or post- analysis to GNN to provide a solution for such cases. Fifth, the SR solvers used in this study were off-the-shelf implementations; custom SR methods specifically designed to exploit graph embeddings could further improve interpretability and efficiency. Sixth, we focus exclusively on a GNN-based latent encoder and do not compare it against alternative embedding strategies such as autoencoders or neural operators. While GNNs are a natural choice in our setting because the data are explicitly organized as spatio-temporal graphs and the target dependencies are local, a rigorous comparison with other latent-space learning paradigms would require substantial architecture-specific development, tuning, and evaluation, and is therefore beyond the scope of the present work. Such a comparison is an important direction for future research, as it may help clarify the trade-offs between relational inductive bias, denoising ability, computational cost, and mesh-independence in physics-informed symbolic regression. Finally, we note that the current framework relies on standard finite-difference approximations for temporal and spatial derivatives. Although higher-order schemes may improve local derivative accuracy, they also require broader neighborhoods and tend to be more sensitive to noise, which is less aligned with the localized graph construction adopted in this work. In our preliminary analysis, these alternatives did not yield substantial downstream improvements in SR performance. Likewise, automatic differentiation is not directly applicable here because derivatives are estimated from discrete observations rather than from a differentiable surrogate model^114,115. Future work may therefore explore hybrid formulations that combine the proposed representation with neural surrogate models, allowing automatic-differentiation-based derivative estimation in a physics-aware symbolic discovery pipeline.

Taken jointly, the findings of this study demonstrate that enriching SR with structured, physics-informed graph representations substantially improves its capacity to recover governing equations from limited and noisy data. The dual representation bridges the gap between purely data-driven and physics-guided modeling by allowing SR models to leverage implicit spatial and temporal structure while retaining full analytical interpretability. As scientific datasets continue to grow in both complexity and scale, the proposed framework offers a promising and generalizable pathway toward more robust, explainable, and physically consistent equation discovery. To this end, future work should extend this approach to real experimental systems, explore joint SR–GNN optimization strategies, and investigate its applicability to multi-field coupled dynamics.

Data availability

The code and data that have been used in this study are available from the corresponding author

Appendix

Graph data transformation example

To illustrate this graph transformation, let us consider a simple, one-dimensional experiment in which a mass is released from rest and allowed to fall freely under the influence of gravity. The motion of the mass is recorded over time, yielding a discrete set of temporal measurements. The primary measured quantity is the height of the mass h(t), while its instantaneous velocity f(t) (corresponding to the vertical component of the velocity field) is either directly measured using sensors or derived numerically from the height data. To this end, the raw dataset can be written as i n line- eq- IEq1 86 , where inline-eq-IEq187 denotes the fixed horizontal coordinate of the ball, inline-eq-IEq188 is the i-th time sample, in l ine-eq-IEq189 is the observed height, in l ine-e q ^-I ₀ is the vertical velocity, and _Eq 19

in l ^ine _Eq191 is the estimated acceleration. _-eq -I As such, for a one-dimensional trajectory, we construct a two-dimensional grid graph over height and time. Let the measurement times be inliⁿ _e-eq-IEq192 with smallest nonzero spacing in li ne-eq-IEq 1 93 (robustly, the instrument time resolution or the mode of spacings under noise). Let the observed heights be inliⁿ _e-eq-IEq194 with smallest nonzero spacing in li ne-eq-I Eq1 9 5 (robustly, the height/position quantization or a percentile-based minimum under noise). Define grid axes in l ine- e q -IE q 1 96 , spanning the data ranges inline -eq-IEq197 and inline -eq-IEq198. The node set is the Cartesian grid i n li_ne- eq -IE_q1 9₉ , where i_nli n e_-eq-IEq200 holds the features available at that grid point. Given the raw samples i n line- eq-I Eq2 01 each record is mapped to its nearest grid node via inli n e-eq-I E_q202 and we set i_nline-eq- I Eq20 3 . If multiple samples map to the same (p, q), we aggregate (e.g., mean or median). Unobserved nodes retain a missing value token and can be imputed if desired. Edges encode local adjacency on the grid. We use 4-neighborhood (von Neumann) connectivity in l ine-eq- IEq204 so that i n li n e-eq-IEq205 .

Features for each equation

Table 5 details the complete set of features generated for each of the ten benchmark equations used in this study. For every dataset, both the true physical variables required by the governing equation and a set of additional synthetic features were included. The physical variables correspond directly to the quantities appearing in the mathematical formulation of each equation. For instance, displacement and stiffness in Hooke’s Law. In addition, to emulate the complexity and redundancy typical of real experimental measurements, three extra non-informative features (inline-eq-IEq206, inline-eq-IEq207, inline-eq-IEq208) were added to each dataset. These distractor features were generated from smooth random functions or stochastic processes unrelated to the target variable. In addition, every dataset included a dedicated Gaussian noise feature (inline-eq-IEq209), independently drawn from a normal distribution i nline-eq-IEq21 0 , serving as an explicit irrelevant input.

Additional noise types for robustness evaluation

The main text evaluates robustness under additive Gaussian perturbations applied to both input and output variables. Here, we extend the same protocol (six noise levels up to in l ine-eq-IEq223) to three additional noise families that frequently occur in physical measurements: heavy-tailed noise, impulsive outliers, and systematic drift. All results are averaged across the ten benchmark systems and reported using the same three metrics (MAE, i^{nline-eq-IEq224}, ETD). The noise profiles are:

• We replace the Gaussian perturbation with a Student-t noise model: inl i ne - eq- IE q 2 2 5 , where inline-eq-IEq226 controls tail heaviness (smaller inline-eq-IEq227 yields more outliers).
• We use a mixture model to emulate sporadic sensor glitches: i n li n e-e q-I E^q ₂2 8 , with inli n e-eq-IEq229 and small contamination probability p.
• We add a slowly varying bias to the target: inlin e -eq- I Eq23 0 , where b(t) is a low-frequency drift term (e.g., linear or sinusoidal), and inline-eq-IEq231 is small additive noise.Table 6 summarizes the statistics at the maximum noise level in l ine-eq-IEq232 for each SR model, comparing the tabular (baseline) and the proposed dual tabular-graph representation.

Physical consistency check

To complement fit-based metrics (MAE/i^{nline-eq-IEq238}/ETD), we quantify whether discovered equations satisfy known physical constraints. For each final discovered expression E, we report two scalar diagnostics computed on the simulated trajectory i^nli: (i) a symmetry/invariance violation score, and (ii) a boundary/initial-condition rene-eq-IEq239 sidual score. We compute inline-eq-IEq240 by sampling a small set of transformations i n line-eq-IEq241 that are valid for the target system (e.g., time translation for autonomous systems; Galilean boosts for applicable dynamical systems). We report the average score across experiments:

Moreover, using inline-eq-IEq242, we report a normalized constraint residual:

where inline-eq-IEq243 is the number of scalar constraint entries (nodes inline-eq-IEq244 times inline-eq-IEq245 fields).

Table 7 summarizes these results. These scores correspond to the same physical-consistency mechanisms used in the SR objective via inline-eq-IEq246 and inline-eq-IEq247 in Eq. (6). Across the ten benchmark systems, the dual tabular-graph representation reduces both types of violations for all SR solvers, with the largest relative gains typically observed for classical genetic-programming and sparse-regression baselines, indicating that the physics-aware graph/latent structure helps steer the search away from algebraically plausible but physically inconsistent expressions.

Performance landscapes over noise and data availability

To provide an intuitive, multi-dimensional view of how robustness and data efficiency interact, we visualize the dual-representation performance over the joint grid of noise intensity inline-eq-IEq251 and number of experiments K. Each point in Fig. 7 corresponds to one inli ne-eq-IEq252 setting, and the vertical axis shows the aggregated metric value (mean across SR solvers and benchmark systems). These 3D views complement the 1D sensitivity curves (Figs. 4–5) by making it easier to identify regimes where performance degrades most strongly (high noise / low data) and how quickly it recovers as more experiments are available.

Received: 17 March 2026; Accepted: 14 May 2026

References

1. Garrison, J. W. Newton and the relation of mathematics to natural philosophy. J. Hist. Ideas 48(4), 609–627 (1987).
2. Guicciardini, N. Isaac Newton on Mathematical Certainty and Method (Mit Press, 2009).
3. Shmuel, A., Glickman, O. & Lazebnik, T. Symbolic regression as a feature engineering method for machine and deep learning regression tasks. Mach. Learn.: Sci. Technol. 5(2), 025065 (2024).
4. Grabiner, J. V. Newton, maclaurin, and the authority of mathematics. Am. Math. Mon. 111(10), 841–852 (2004).
5. Kvasz, L. Revisiting the mathematisation thesis: Galileo, descartes, newton, and the language of nature. Int. Stud. Philos. Sci. 30(4), 399–406 (2016).
6. Redish, E. F. Using math in physics: Overview. Phys. Teach. 59(5), 314–318 (2021).
7. Yeargers, E. K., Herod, J. V. & Shonkweiler, R. W. An Introduction to the Mathematics of Biology: With Computer Algebra Models (Springer Science & Business Media, 2013).
8. Debreu, G. The mathematization of economic theory. Am. Econ. Rev. 81(1), 1–7 (1991).
9. Gingras, Y. What did mathematics do to physics? Hist. Sci. 39(4), 383–416 (2001).
10. National Research Council, Division on Engineering, Physical Sciences, Board on Mathematical Sciences, Their Applications, Committee on the Mathematical Sciences in, and 2025. Fueling Innovation and Discovery: The Mathematical Sciences in the 21st Century (National Academies Press, 2012).
11. Polya, G. Mathematical Discovery, 1962 (Wiley, 1962).
12. Abrahamson, D., Gutiérrez, J., Charoenying, T., Negrete, A. & Bumbacher, E. Fostering hooks and shifts: Tutorial tactics for guided mathematical discovery. Technol. Knowl. Learn. 17, 61–86 (2012).
13. Stam, J., Stokhof, M. & Van Lambalgen, M. Naturalising mathematics? A wittgensteinian perspective. Philosophies 7(4), 85 (2022).
14. Holton, G., Chang, H. & Jurkowitz, E. How a scientific discovery is made: A case history. Am. Sci. 84(4), 364–375 (1996).
15. Nickles, T. Scientific Discovery, Logic, and Rationality (Springer Science & Business Media, 2012).
16. Merton, R. K. The Sociology of Science: Theoretical and Empirical Investigations (University of Chicago press, UK, 1973).
17. De Solla Price, D. J. Little Science, Big Science (Columbia University Press, 1963).
18. Stevens, R., Taylor, V., Nichols, J., Maccabe, A.B., Yelick, K. & Brown, D. Ai for science: Report on the department of energy (doe) town halls on artificial intelligence (ai) for science. Technical report, (Argonne National Lab.(ANL), Argonne, IL (United States), 2020).
19. Van Noorden, R. & Perkel, J. M. Ai and science: what 1,600 researchers think. Nature 621(7980), 672–675 (2023).
20. Krenn, M. et al. On scientific understanding with artificial intelligence. Nat. Rev. Phys. 4(12), 761–769 (2022).
21. Grossmann, I. et al. Ai and the transformation of social science research. Science 380(6650), 1108–1109 (2023).
22. Alyass, A., Turcotte, M. & Meyre, D. From big data analysis to personalized medicine for all: Challenges and opportunities. BMC Med. Genomics 8, 1–12 (2015).
23. Einav, L. & Levin, J. Economics in the age of big data. Science 346(6210), 1243089 (2014).
24. Makke, N. & Chawla, S. Interpretable scientific discovery with symbolic regression: A review. Artif. Intell. Rev. 57(1), 2 (2024).
25. Kronberger, G., Burlacu, B., Kommenda, M., Winkler, S. M. & Affenzeller, M. Symbolic Regression (CRC Press, 2024).
26. Wang, Y., Wagner, N. & Rondinelli, J. M. Symbolic regression in materials science. MRS Commun. 9(3), 793–805 (2019).
27. Tian, Y. et al. Interactive symbolic regression with co-design mechanism through offline reinforcement learning. Nat. Commun. 16(1), 3930 (2025).
28. Dong, J. & Zhong, J. Recent advances in symbolic regression. ACM Comput. Surv. 57(11), 1–37 (2025).
29. Angelis, D., Sofos, F. & Karakasidis, T. E. Artificial intelligence in physical sciences: Symbolic regression trends and perspectives. Arch. Comput. Methods Eng. 30(6), 3845–3865 (2023).
30. Makke, N. & Chawla, S. Data-driven discovery of tsallis-like distribution using symbolic regression in high-energy physics. PNAS Nexus 3, (2024).
31. Tenachi, W., Ibata, R. & Diakogiannis, F. I. Deep symbolic regression for physics guided by units Constraints: Toward the automated discovery of physical laws. Astrophys. J. 959, 99 (2023).
32. Udrescu, S. & Tegmark, M. Ai feynman: a physics-inspired method for symbolic regression. Sci. Adv. 6, (2020).
33. Langley, P. Bacon: A production system that discovers empirical laws. In IJCAI, page 344. Citeseer, (1977).
34. Langley, P. Scientific Discovery: Computational Explorations of the Creative Processes (MIT press, 1987).
35. Petersen, B. K., Landajuela, M., Mundhenk, T. N., Santiago, C. P., Kim, S. K. & Kim, J. T. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. (2021).
36. Kubalík, J., Derner, E. & Babuška, R. Toward physically plausible data-driven models: A novel neural network approach to symbolic regression. IEEE Access 11, 61481–61501 (2023).
37. Kubalík, J., Derner, E. & Babuška, R. Multi-objective symbolic regression for physics-aware dynamic modeling. Expert Syst. Appl. 182, 115210 (2021).
38. Zhang, R., Liu, Y. & Sun, H. Physics-informed multi-lstm networks for metamodeling of nonlinear structures. Comput. Methods Appl. Mech. Eng. 369, 113226 (2020).
39. Cuomo, S. et al. Scientific machine learning through physics-informed neural networks: Where we are and what’s next. J. Sci. Comput. 92(3), 88 (2022).
40. Huang, B. & Wang, J. Applications of physics-informed neural networks in power systems-a review. IEEE Trans. Power Syst. 38(1), 572–588 (2022).
41. Miao, H., Zhao, Y., Guo, C., Yang, B., Zheng, K., Huang, F., Xie, J. & Jensen, C.S. A unified replay-based continuous learning framework for spatio-temporal prediction on streaming data. In 2024 IEEE 40th International Conference on Data Engineering (ICDE) 1050–1062 (IEEE, 2024).
42. Zhang, H., Zhang, W., Miao, H., Jiang, X., Fang, Y. & Zhang, Y. Strap: Spatio-temporal pattern retrieval for out-of-distribution generalization. arXiv preprint arXiv:2505.19547 (2025). link
43. Liu, C., Yang, S., Xu, Q., Li, Z., Long, C., Li, Z. & Zhao, R. Spatial-temporal large language model for traffic prediction. In 2024 25th IEEE International Conference on Mobile Data Management (MDM) 31–40 (IEEE, 2024). link
44. Liu, C., Hettige, K. H., Xu, Q., Long, C., Xiang, S., Cong, G., Li, Z. & Zhao, R. St-llm+: Graph enhanced spatio-temporal large language models for traffic prediction. IEEE Trans. Knowl. Data Eng. (2025).
45. Taskin, B., Xie, W. & Lazebnik, T. Knowledge integration for physics-informed symbolic regression using pre-trained large language models. Sci. Rep. (2026).
46. Icke, I. & Bongard, J. Modeling hierarchy using symbolic regression. In 2013 IEEE Congress on Evolutionary Computation 2980– 2987 (2013).
47. Kammerer, L., Kronberger, G., Burlacu, B., Winkler, S., Kommenda, M. & Affenzeller, M. Symbolic regression by exhaustive search: reducing the search space using syntactical constraints and efficient semantic structure deduplication. Genetic and Evolutionary Computation 79–99 (2020).
48. Valipour, M., You, B., Panju, M. & Ghodsi, A. Symbolicgpt: a generative transformer model for symbolic regression. arXiv preprint arXiv:2106.14131 (2021). link
49. Weisberg, S. Applied Linear Regression (Wiley, 2005). link
50. Makke, N. & Chawla, S. Interpretable scientific discovery with symbolic regression: A review. Artif. Intell. Rev. 57, (2024).
51. Chen, C., Luo, C. & Jiang, Z. Block building programming for symbolic regression. Neurocomputing 275, 1973–1980 (2018).
52. Jin, Y., Fu, W., Kang, J., Guo, J. & Guo, J. Bayesian symbolic regression. (2019).
53. Tohme, T., Liu, D. & Youcef-Toumi, K. Gsr: A generalized symbolic regression approach. (2022).
54. Korns, M. F. A baseline symbolic regression algorithm. In Genetic Programming Theory and Practice X 117–137 (Springer, 2013).
55. Riolo, R. Genetic Programming Theory and Practice X (Springer, 2013).
56. Heule, M. J. H. & Kullmann, O. The science of brute force. Commun. ACM 60(8), 70–79 (2017).
57. Brunton, S. L., Proctor, J. L. & Kutz, J. N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. 113(15), 3932–3937 (2016).
58. Kaptanoglu, A. A., de Silva, B. M., Fasel, U., Kaheman, K., Goldschmidt, A. J., Callaham, J. L., Delahunt, C. B., Nicolaou, Z. G., Champion, K., Loiseau, J.-C. et al. Pysindy: A comprehensive python package for robust sparse system identification. arXiv preprint arXiv:2111.08481 (2021). link
59. Cava, W. L., Orzechowski, P., Burlacu, B., de França, F. O., Virgolin, M., Jin, Y., Kommenda, M. & Moore, J. H. Contemporary symbolic regression methods and their relative performance. In Thirty-fifth Conference on Neural Information Processing Systems (NeurIPS 2021) Datasets and Benchmarks Track, (2021). link
60. Petersen, B. K., Landajuela, M., Mundhenk, T. N., Santiago, C., Kim, S. K. & Kim, J. T. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. (2019).
61. Koza, J. R. Genetic Programming: On the Programming of Computers by Means of Natural Selection (MIT Press, 1992).
62. Orzechowski, P., Cava, W. L. & Moore, J. H. Where are we now? In Proceedings of the Genetic and Evolutionary Computation Conference (2018).
63. Virgolin, M., Alderliesten, T., Witteveen, C. & Bosman, P. A. N. Improving model-based genetic programming for symbolic regression of small expressions. Evol. Comput. 29, 211–237 (2021).
64. Stephens, T. et al. Genetic programming in python, with a scikit-learn inspired api: gplearn. Documentation at https://gplearn.readthedocs.io/en/stable/intro.html, (2016). link
65. Neumann, P., Cao, L., Russo, D., Vassiliadis, V. S. & Lapkin, A. A. A new formulation for symbolic regression to identify physicochemical laws from experimental data. Chem. Eng. J. 387, 123412 (2020). link
66. Kim, S. et al. Integration of neural network-based symbolic regression in deep learning for scientific discovery. (2019).
67. Jiang, N. & Xue, Y. Racing control variable genetic programming for symbolic regression. In Proceedings of the AAAI Conference on Artificial Intelligence vol. 38, pp. 12901–12909 (2024).
68. Kubalik, J., Derner, E. & Babuska, R. Multi-objective symbolic regression for physics-aware dynamic modeling. Expert Syst. Appl. 182, 115210 (2021).
69. Gordin, S. A., Ermakov, A. N., Zakharov, A. Y. & Wang, J. Determining the maximum linear mass of a suspended conveyor belt using pysr symbolic regression. Mining 5(4), 83 (2025).
70. Eichhorn, S., Mohapatra, A. & Goebel, C. Pisr: Physics-informed symbolic regression for predicting power system voltage. In Proceedings of the 16th ACM International Conference on Future and Sustainable Energy Systems 92–107 (2025).
71. Chen, Z., Liu, Y. & Sun, H. Physics-informed learning of governing equations from scarce data. Nat. Commun. 12(1), 6136 (2021).
72. Udrescu, S.-M. & Tegmark, M. Ai feynman: A physics-inspired method for symbolic regression. Sci. Adv. 6(16), eaay2631 (2020). link
73. Martius, G. & Lampert, C. H. Extrapolation and learning equations. arXiv preprint arXiv:1610.02995 (2016). link · link
74. Golden, M. Scalable sparse regression for model discovery: The fast lane to insight. arXiv preprint arXiv:2405.09579 (2024). link · link
75. Cai, S., Mao, Z., Wang, Z., Yin, M. & Karniadakis, G. E. Physics-informed neural networks (pinns) for fluid mechanics: A review. Acta. Mech. Sin. 37(12), 1727–1738 (2021). link
76. Cuomo, S. et al. Scientific machine learning through physics-informed neural networks: Where we are and what’s next. J. Sci. Comput. 92(3), 88 (2022).
77. Cai, S., Wang, Z., Wang, S., Perdikaris, P. & Karniadakis, G. E. Physics-informed neural networks for heat transfer problems. J. Heat Transfer 143(6), 060801 (2021).
78. Hanna, J. M., Talbot, H. & Vignon-Clementel, I. E. Improved physics-informed neural networks loss function regularization with a variance-based term. arXiv preprint arXiv:2412.13993 (2024). link
79. Farea, A., Yli-Harja, O. & Emmert-Streib, F. Understanding physics-informed neural networks: Techniques, applications, trends, and challenges. AI 5(3), 1534–1557 (2024). link
80. Lawal, Z. K., Yassin, H., Lai, D. T. C. & Idris, A. C. Physics-informed neural network (pinn) evolution and beyond: A systematic literature review and bibliometric analysis. Big Data Cogn. Comput. 6(4), 140 (2022).
81. Luo, K. et al. Physics-informed neural networks for pde problems: A comprehensive review. Artif. Intell. Rev. 58(10), 1–43 (2025).
82. Chen, Y. & Koohy, S. Gpt-pinn: Generative pre-trained physics-informed neural networks toward non-intrusive meta-learning of parametric pdes. Finite Elem. Anal. Des. 228, 104047 (2024).
83. Chen, Y., Ji, Y., Narayan, A. & Zhenli, X. Tgpt-pinn: Nonlinear model reduction with transformed gpt-pinns. Comput. Methods Appl. Mech. Eng. 430, 117198 (2024).
84. Ranasinghe, N., Xia, Y., Seneviratne, S. & Halgamuge, S. Ginn-kan: Interpretability pipelining with applications in physics informed neural networks. arXiv preprint arXiv:2408.14780 (2024). link
85. Pareja, A. et al. Evolvegcn: Evolving graph convolutional networks for dynamic graphs. In Proceedings of the AAAI Conference on Artificial Intelligence (2020). link
86. Zheng, Y., Yi, L. & Wei, Z. A survey of dynamic graph neural networks. arXiv preprintarXiv:2404.18211 (2024). link
87. Zhu, Y., Xu, W., Zhang, J., Du, Y., Zhang, J., Liu, Q., Yang, C. & Wu, S. A survey on graph structure learning: Progress and opportunities. arXiv preprint arXiv:2103.03036 (2021). link · link
88. Kipf, T., Fetaya, E., Wang, K.-C., Welling, M. & Zemel, R. Neural relational inference for interacting systems. In Proceedings of the 35th International Conference on Machine Learning (ICML) (2018). link
89. He, H., Yu, X., Zhang, J., Song, S. & Letaief, K. B. Message passing meets graph neural networks: A new paradigm for massive mimo systems. IEEE Trans. Wireless Commun. 23(5), 4709–4723 (2023).
90. Keren, L. S., Liberzon, A. & Lazebnik, T. A computational framework for physics-informed symbolic regression with straightforward integration of domain knowledge. Sci. Rep. 13(1), 1249 (2023).
91. Tenachi, W., Ibata, R. & Diakogiannis, F. I. Deep symbolic regression for physics guided by units constraints: Toward the automated discovery of physical laws. Astrophys. J. 959(2), 99 (2023).
92. Cranmer, M. Pysr: High-performance symbolic regression in python and Julia. ascl-2409 (Astrophysics Source Code Library, 2024).
93. Brunton, S. L., Proctor, J. L. & Kutz, J. N. Discovering governing equations from data by sparse identification of nonlinear dynamics (sindy). Proc. Natl. Acad. Sci. 113(15), 3932–3937 (2016).
94. Kim, S. et al. Integration of neural network-based symbolic regression in deep learning for scientific discovery. IEEE Trans. Neural Netw. Learn. Syst. 32(9), 4166–4177 (2020).
95. Icke, I. & Bongard, J. C. Improving genetic programming based symbolic regression using deterministic machine learning. In 2013 IEEE Congress on Evolutionary Computation 1763–1770 (IEEE, 2013).
96. Liang, X. et al. Physics-informed neural network for chiller plant optimal control with structure-type and trend-type prior knowledge. Appl. Energy 390, 125857 (2025).
97. Rudy, S. H., Brunton, S. L., Proctor, J. L. & Kutz, J. N. Data-driven discovery of partial differential equations. Sci. Adv. 3(4), e1602614 (2017).
98. Rychlewski, J. On hooke’s law. J. Appl. Math. Mech. 48(3), 303–314 (1984).
99. De Lima, J. A. S. & Santos, J. Generalized stefan-boltzmann law. Int. J. Theor. Phys. 34(1), 127–134 (1995).
100. Rushka, M. & Freericks, J. K. A completely algebraic solution of the simple harmonic oscillator. Am. J. Phys. 88(11), 976–985 (2020).
101. Sunday, J. The duffing oscillator: Applications and computational simulations. Asian Res. J. Math. 2(3), 1–13 (2017).
102. Cannon, J. R. The One-Dimensional Heat Equation (Cambridge University Press, 1984).
103. Hon, Y.-C. & Mao, X. Z. An efficient numerical scheme for burgers’ equation. Appl. Math. Comput. 95(1), 37–50 (1998).
104. Taur, Y. Analytic solutions of charge and capacitance in symmetric and asymmetric double-gate mosfets. IEEE Trans. Electron Devices 48(12), 2861–2869 (2002).
105. Swinehart, D. F. The beer-lambert law. J. Chem. Educ. 39(7), 333 (1962).
106. Arino, J., Wang, L. & Wolkowicz, G. S. K. An alternative formulation for a delayed logistic equation. J. Theor. Biol. 241(1), 109–119 (2006).
107. Berezansky, L. & Braverman, E. Mackey-glass equation with variable coefficients. Comput. Math. Appl. 51(1), 1–16 (2006).
108. Wang, G., Wang, E., Li, Z., Zhou, J. & Sun, Z. Exploring the mathematic equations behind the materials science data using interpretable symbolic regression. Interdiscip. Mater. 3(5), 637–657 (2024).
109. Guo, J. & Yin, W.-J. Harnessing data using symbolic regression methods for discovering novel paradigms in physics. Sci. China Phys. Mech. Astron. 67(6), 267301 (2024).
110. AbdusSalam, S., Abel, S. & Romão, M. C. Symbolic regression for beyond the standard model physics. Phys. Rev. D 111(1), 015022 (2025).
111. Zhang, S., Tong, H., Jiejun, X. & Maciejewski, R. Graph convolutional networks: A comprehensive review. Comput. Soc. Netw. 6(1), 1–23 (2019).
112. Karniadakis, G. E. et al. Physics-informed machine learning. Nat. Rev. Phys. 3(6), 422–440 (2021). link
113. Brandstetter, J., Worrall, D. & Welling, M. Message passing neural pde solvers. arXiv preprint arXiv:2202.03376 (2022). link
114. Recktenwald, G. W. Finite-difference approximations to the heat equation. Mech. Eng. 10(01), (2004). link
115. Margossian, C. C. A review of automatic differentiation and its efficient implementation. Wiley Interdiscip. Rev.: Data Mining Knowl. Discov. 9(4), e1305 (2019). Author contributions Teddy Lazebnik: Conceptualization, Software, Formal analysis, Methodology, Investigation, Writing - Original Draft, Visualization. Alex Liberzon: Conceptualization, Validation, Writing - Review & Editing. Funding Open access funding provided by Jönköping University. Declarations Competing interests The authors declare no competing interests. Additional information Correspondence and requests for materials should be addressed to T.L. Reprints and permissions information is available at www.nature.com/reprints.© The Author(s) 2026