arXiv CSopen access

Discovering Multiscale Deep Formulas in Complex Systems via Neural-Guided Lambda Calculus

Discovering Multiscale Deep Formulas in Complex Systems via Neural-Guided Lambda Calculus Hanqiao Yu1,2† , Shusen Yang1,2*† , Xuebin Ren1,3*, Cong Zhao1,2

arXiv:2606.07426v1 [cs.LG] 5 Jun 2026

1*

National Engineering Laboratory for Big Data Analytics, Xi’an Jiaotong University, Xi’an, 710049, Shaanxi, China. 2* School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, 710049, Shaanxi, China. 3* School of Computer Science and Technology, Faculty of Electronic and Information Engineering, Xi’an Jiaotong University, Xi’an, 710049, Shaanxi, China.

*Corresponding author(s). E-mail(s): [email protected]; [email protected]; † These authors contributed equally to this work. Abstract A fundamental problem in science is identifying underlying patterns of complex systems in the form of concise mathematical formulas. Current Artificial Intelligence (AI)-based methods have shown strong performance in single-scale systems, yet remain limited in identifying scale-specific formulas in multiscale complex systems. We present Deflex, an end-to-end AI method to automatically extract multiscale formulas with potentially different forms, including invariants and distributions, from complex systems. Deflex consists of two subsystems named Deflexformer and Deflexpressor. Deflexpressor is a lambda-calculus symbolic regression model for higher-order formulas. Deflexformer is a decomposable deep energy model for learning unified representations across scales. Deflexpressor generates synthetic data to pre-train Deflexformer, which then guides formula discovery by decoupling multiscale latent relationships. Across six representative complex systems with diverse behaviors, Deflex achieves up to 7-fold higher efficiency than the state-of-the-art methods while enabling automated multiscale discovery. Our work could be a useful tool for scientific discovery across disciplines.

1

1 Introduction Discovering rules in mathematical forms, such as energy conservation and Boltzmann distribution, in nature and society has always been an important endeavor in scientific research [1, 2]. In recent years, artificial intelligence (AI) has played an increasingly significant role in identifying mathematical formulas from different systems [2–5]. From computational and cognitive perspectives, AI-based automatic formula discovery for complex systems has to solve three fundamental scientific problems: 1) the volume challenge, indicating that the huge number of interacting elements would result in computational and representational prohibitiveness of formula findings, 2) search-space explosion, meaning that all possible free-form formulas increase exponentially as the number of variables and potential expression depth increase, due to the volume and complex nonlinear interactions in complex systems [3, 4, 6–8], and 3) the scale gap issue, suggesting that the mathematical rules emerge at different scales in complex systems and automatic discovery requires unified representations and analysis at different levels [9, 10]. Facing the above challenges, discovering genuine mathematical formulas in a fully automated fashion, without pre-specified equation structures, remains largely limited to single-scale and few-body systems of lower complexity [1, 6, 11–15]. Meanwhile, increasing effort has been put into research on complex systems, such as fluid dynamics [4], crowd movement [16], climate evolution [10, 17], and social systems [8], which explain many complex phenomena and behaviors in the real world. Many distinct properties of complex systems come from the transition to chaos and the emergence of order from disorder across different spatial and temporal scales [18, 19]. Such chaos and order emergence also introduces different levels of regularity, requiring scientific analysis to account for distinct mathematical forms: deterministic descriptions in conservation laws and probability distributions for stochastic behaviors [20]. AI methods have demonstrated great potential for discovering the governing rules behind the interactions of massive elements at multiple scales in these complex systems. For instance, their capacity to efficiently approximate solutions for complex known equations in physics [21], their ability to discern intricate patterns from largescale observational data in domains such as meteorology or climate science [22], and the emerging success of large models in learning or representing complex mathematical relationships from scientific data [23] all underscore this potential. Unfortunately, no existing AI methods can achieve automatic discovery of multiscale deep formulas in complex systems. While methods like sparse regression can effectively learn mathematical patterns from large datasets [13, 24], their reliance on pre-defined functional forms limits the exploration of unknown regularities. Symbolic regression (SR), based on symbolic trees and genetic programming [2, 6, 7, 13, 25], is a powerful tool for automatically discovering mathematical formulas with unknown forms from data. However, traditional SR methods suffer from low success rates and computational efficiency for identifying complex relationships. In particular, discovering probability distributions from raw observational data—rather than merely fitting predefined distribution functions—represents a core challenge for multiscale formula discovery, as stochastic patterns naturally emerge at coarse-grained scales [20]. While some advanced methods (e.g., AI Feynman [6]) can handle distribution fitting, they

2

require substantial manual specification of candidate functional forms. Recently, substantial progress has been made in symbolic regression, which has greatly enhanced the discovery of laws in real-world data. For instance, physics-inspired symbolic regression has explored incorporating physical knowledge (e.g., higher-order symmetries) as priors [26, 27]. Also, neural network-integrated SR discovers novel patterns from representations learned from deep models (e.g., GNNs) [13]. Besides, many studies have enabled symbolic regression to find complex relationships on benchmark datasets by improving the searching strategies with evolutionary algorithms [6, 28], reinforcement learning [29], or Monte Carlo tree search [27]. Nonetheless, without higher-order representations and tailored design, the automatic discovery of the above studies can hardly work for multiscale complex systems. In this article, we propose a novel framework of Deep Formula discovery for Complex systems (Deflex) by incorporating new symbolic regression systems with strong expressiveness, unified formulation methodology, and decomposable neural network architecture. The main contributions of our work are summarized as follows: We first regard all mathematical formulas as probability distributions [30], which can unify the system patterns at different scales. Specifically, equations and invariants are treated as distributions with very narrow peaks, which indicates their tiny errors or uncertainty. Original distribution-based formulas are treated as distributions with wider peaks. We then formulate these distribution patterns through energy-based models (EBMs), where each pattern is represented as an energy function. This formulation facilitates both mathematical modeling and algorithmic processing of these patterns. Particularly, by measuring the likelihood of EBMs over the observation data, we build a tractable learning framework for finding the optimal energy functions. With these unified representations and formulations, the overall workflow of Deflex is shown in Fig. 1. We next develop Deflexformer, the neural network subsystem of Deflex. Deflexformer is a deep neural network (DNN) for learning the energy functions. It is designed as a self-attention-based network (Fig. 2a,b) to leverage the demonstrated efficacy of self-attention architectures in processing set-structured data and capturing the patterns for both element-wise and global variables. Different from existing neural-integrated methods, Deflexformer can further decompose the learned functions into elementary components that are tractable for symbolic regression. Deflexformer is composed of multiple Transformer-like blocks, each of which represents a distinct elementary mathematical transformation. Thus, any intricate mathematical relationships learned by Deflexformer, can be collectively parameterized by the modular transformations of all Deflexformer blocks. The neural network with decomposable blocks can efficiently guide symbolic regression algorithms for tractable learning of mathematical formulas from each component, rather than from the whole. However, traditional tree-based symbolic regression may have significant limitations in both expression ability and discovery efficiency for higher-order mathematical relationships (such as summation or mapping) over a large volume of interacting elements in complex systems. The former limitation is lacking native support for variable bindings, so a standard expression tree must either rely on a rigid, predefined summation primitive or force the structure to unroll into an intractably long arithmetic chain of pairwise additions. The latter is due to the combinatorial search

3

space of potentially massive arguments for complex systems. To address these limitations, we present Deflexpressor, a symbolic regression subsystem that leverages lambda-calculus expressions. The inherent capabilities of lambda calculus for function abstraction, variable binding, and higher-order functions provide a powerful and flexible representational framework. This allows Deflexpressor to express and discover free-form higher-order mathematical relationships, such as mapping and aggregating interactions over distinct pairs without pre-specifying the operators. In addition, Deflexpressor can efficiently navigate the search space for complex symbolic structures via two novel designs. One is to guide lambda expression construction by incorporating a rule-based generator and an auto-regressive neural network. The other optimizes the efficiency-diversity trade-offs of expression generation and mutation by extending genetic programming mutations with transformations specified for lambda calculus (e.g., inspired by beta-reduction or exploring variations in function arity and application). The workflow of our approach proceeds as follows. First, Deflexpressor generates a large number of formulas for sampling a synthetic dataset (Fig. 1a). The synthetic dataset is then used to pre-train a Deflexformer block, which is a Transformer-like neural network (Fig. 1b). This pre-training serves as prior knowledge for the subsequent step, where we cascade the multiple copies of the pre-trained block as the whole Deflexformer to learn patterns from the observation dataset of a certain complex system (Fig. 1c). During learning, Deflexformer employs element-wise and time-wise self-attention mechanisms to capture the interactions among different entities and time frames. Each trained block then essentially represents a distinct component of the latent formulas of the patterns. In this phase, the observation dataset, together with some random data, is then fed to the trained Deflexformer for inference, which generates a large volume of data representations output by all blocks (Fig. 1d). Finally, the representation dataset output by the elementary blocks in Deflexformer is fitted using Deflexpressor to obtain the desired formulas from each component, through iterative formula generation and mutation in the form of lambda expressions (Fig. 1e). Thus far, we have outlined the design principles and general processes of the proposed method, which aim to address the aforementioned challenges for AI-based formula discovery from complex systems. Concerning the volume problem, Deflex uses a symbolic regression system based on lambda calculus to compress the mathematical rules into concise forms, and utilizes scalable deep learning and symbolic regression algorithms to handle large-scale datasets. For the complexity of interactions, we design a self-attention-based DNN for Deflex to decompose the complex systems into different scales and sub-problems, making the approach efficient for highly complex formulas. For the scale gap problem, Deflex uses energy-based models to unify invariant and distributional patterns across scales. A detailed illustration of the system modelling, the subsystems design, and the specific workflow will be presented in the Method Section.

2 Results We evaluate the effectiveness and efficiency of Deflex in six representative complex systems with different behaviors. Particularly, these systems can be categorized into

4

simulation experiments, including particle motion, two-phase water particles, and fluid dynamics, and the real-world data of human and bird collective movements. More details about the experiments can be found in the SI.

2.1 Compared methods We benchmark our approach against several leading AI methods for formula discovery, including state-of-the-art SR method Operon [31], typical physics-inspired SR method SciMED [27], leading sparse regression method SINDy [24], and AI Feynman 2.0 [6], as well as classic SR methods gplearn [32], PySR [33] and DEAP [34] in terms of accuracy, efficiency, and adaptability to multiscale systems. Different from Deflex with energy models and lambda calculus, existing AI methods (mainly based on symbolic regression) often support limited forms of problems and formulas. To maximize comparability, we made several accommodations for these compared methods. For the comparatives (PySR, AI Feynman, and DEAP) lacking native support for higher-order operators (e.g., map or sum), we added these operators to their symbol set manually. For the rest ones (Operon, SciMED, SINDy, and gplearn) that are incompatible with these operators, we provided an alternative. That is, if such an operator was at the top level of the target formula, we applied the method for the sub-expressions only and then handled the top-level operators within the loss function manually. Similarly, for methods unable to approximate distributions, we provided custom energy-based loss functions (see SI for comparison details). Note that some of these frameworks have variants that would gain better performance with scientists in the loop (e.g., manually providing feature sets or constraints on the expression tree for SciMED and AI Feynman), while we mainly used the raw version for simplicity despite the above adaptations. A method was considered to have failed a task if it still could not discover the rule after these best-effort accommodations.

2.2 Evaluation metrics Since uncovering formulas within a complex system may represent an open-ended challenge, we incorporate qualitative assessments, basic quantitative analyses, renormalization and coarse-graining techniques [35] in the performance evaluation. For established or empirical formulas, we first qualitatively assess whether the structure and parameters of the discovered formulas align with the ground truth. The formulas of invariants are considered successful recovery if they can be expressed as the correct formula with an additional stochastic residual term. Table 1 lists all discovered formulas by Deflex in the studied six systems, compared with the ground truths. Note that, for undocumented discoveries, we provide justifications for their possible real-world interpretations. For reference, we also summarize the capability of the compared methods for discovering these formulas in Table 2. To quantitatively evaluate the accuracy of discovered formulas, we consider two cases based on the different forms of formulas separately. • For the formulas in the form of equations and invariants, we employ the relative mean square error (RMSE) between the prediction and the true value. Specifically, all compared methods are given the same convergence threshold. Note that, in such 5

a case, the RMSEs of different methods are hardly comparable since there is almost no structural difference (except for residual terms with negligible errors) in the commonly identified formulas in these forms. • For distribution-based formulas, we use the negative log-likelihood (NLL) loss defined in Eq. (8) in the method section, which measures the likelihood of the distribution fitting real data. The feasible methods were also set with the same convergence constraint. Note that we also report the RMSE to measure the fitting error of simple distributions for reference. Fig. 3a compares the NLL loss of Deflex’s distribution discovery in comparison with other methods, and reports the RMSE on the commonly identified equations and invariants. Besides, we compare the recovery rate to empirically measure the percentage of times certain formulas can be consistently recovered from the data. To evaluate the efficiency of Deflex, we run Deflex and other comparative algorithms on a workstation with 32 CPU cores working in parallel, and two Tesla V100 GPUs for GPU-accelerated computing (details shown in SI Sec. 3). We compare the time efficiency of Deflex with that of other methods by recording the run-time of finding the same or equivalent formulas from the same dataset, within the same tolerable error in NLL or RMSE. Note that the run-time of Deflex is measured from applying Deflex on the observation dataset, including post-training and inference of Deflexformer, and hierarchical symbolic regression for formula discovery. It excludes the time cost of the Deflexformer block pre-training, which was conducted offline over 10,000 sampled pairs and independent of any studied system. A discovery is considered successful when its error remains within a small margin of the residual error, defined as the minimum achievable discrepancy when fitting the ground truth formula to observations [6]. This baseline accounts for irreducible noise from measurements and stochastic fluctuations in the system (details in SI). Fig. 3b shows the comparison results. In our experiments, we employ different observation scales by varying the spatiotemporal scope of data sampling. Small-scale observations use localized sampling (e.g., selecting spatiotemporally neighboring data points) to capture fine-grained interactions and facilitate the discovery of governing equations, while large-scale observations employ uniform random sampling across the entire domain to reveal macroscopic statistical patterns such as distributions (see SI for detailed sampling protocols). One of the most important characteristics of complex systems is the relevance between scales and behaviors. For evaluating the effectiveness of multiscale formula discovery, we introduce the techniques of coarse-graining and renormalization [35] in our experiments. Coarse-graining refers to grouping the small (fine-grained) system elements {ef } into a set of large (coarse-grained) entities {ec }, e.g., aggregating adjacent particles into clusters or merging smaller fluid grid cells into larger ones. Renormalization means transforming the formula P ({ef }, J) with parameters J describing a system at fine-grained level elements {ef }, into the formula P ({ec }, J ′ ) in the same form, with the only difference in parameters J ′ describing coarse-grained level entities {ec }. Scale-invariant and self-similar systems often satisfy formula renormalization across multiple different scales. In our experiments, we also set up the 6

experimental environments for water molecules and fluid dynamics to evaluate the performance of our Deflex framework on such systems. Specifically, we measure the distance between our discovered formulas at different scales with the normalized ground-truth formulas at corresponding scales. Here, since all the formulas discovered by Deflex can be represented in the probabilistic form, we use Earth Mover’s Distance (EMD) [36] as an alternative error measure. EMD is particularly useful for distributions with very sharp peaks and small overlaps. EM D(µ, ν) =

inf

γ∈Γ(µ,ν)

E(x,y)∼γ d(x, y)p ,

(1)

where µ and ν are two distributions, γ is the set of all possible joint distributions of µ and ν, and d(x, y) is the distance between x and y. Fig. 4a compares the EMD value across scales in the scenarios of water particle motion and fluid dynamics.

2.3 Particle motion The study of the motion of particles, like pollen, molecules, and fundamental particles, has contributed to some breakthrough discoveries, such as the law of energy conservation, the Maxwell distribution, the Boltzmann distribution, and the Langevin equation. Many of these have helped reveal fundamental understandings of how patterns emerge in systems across scales [2, 37]. For instance, the Langevin equation manifests that small-scale dynamics can be approximated as noise at larger scales. We first applied Deflex to the particle motion systems of rare gases and pollen in water [37]. In the rare gas particle motion experiment, we simulated the motion of argon particles at three temperature points near the critical point (Tc , Tc ± 1K, where Tc = 150.72K) within a cube, which allowed us to capture the system’s behavior across different thermodynamic states. When the particles approached the cube edge, they experienced a repulsive force analogous to intermolecular forces. See SI 2.1 for a detailed experimental setup. In the pollen experiment in water, we simulated water and pollen within a microtube. This represents a typical two-scale system, involving microscopic water molecules and mesoscopic pollen. We sampled 16 to 2048 particles for 50 time frames as a data item and trained our network with a dataset of 65,536 sampled data items. For each particle at a time frame, the mass, position, and velocity are embedded into the embedding layer of 64 dimensions. As shown in Table 1, Deflex successfully identifies a series of formulas, including the energy conservation, Boltzmann distribution, Maxwell-Boltzmann distribution (MBD), MBD with Lennard-Jones potential, and the Langevin equation, with only one training pass of Deflexformer in each of these particle motion systems. Notably, since temperature is provided as an input variable in the argon system, the discovered formulas naturally capture the temperature-dependent behaviors across different thermodynamic states, with temperature T appearing explicitly in the energy expressions. In particular, MBD with Lennard-Jones potential and the Langevin equation cannot be discovered by other compared methods in our experiments (shown in Table 2). And those tree-based symbolic regression method requires manual specification MBD with Lennard-Jones potential identified by Deflex in the dynamics of rare gas particles is

7

a nonlinear equation that consists of up to 23 nodes, demonstrating Deflex’s ability to discover complex nonlinear equations. On the commonly identified Maxwell-Boltzmann distribution, Deflex manages to achieve the lowest NLL loss among the compared methods in both Argon and Pollen particles (shown in Fig. 3a). Meanwhile, the RMSEs of all compared methods are nearly the same on both energy conservation and momentum conservation in the two systems. This is because the structures of correctly identified formulas are almost identical in these methods, as mentioned before. As shown in Fig. 3b, Deflex incurs significantly less computational time than most of the comparatives. SINDy and Operon are quite close to or sometimes faster than Deflex, but only limited to some lowerorder formulas like momentum conservation. This is due to their limitations in formula structures. Furthermore, ablation results indicate that our Langevin sampling strategy is essential for efficiently identifying these complex equilibrium distributions (Fig. 5g).

2.4 Human and bird mobility The movement of human crowds and bird flocks exemplifies significant randomness and multiscale characteristics in complex systems. In the experiment, we used data from the Zurich Carnival dataset [38] and the eBird bird migration dataset [39] to simulate the collective movement of humans and birds, respectively. These datasets consist of the moving trajectories of humans and birds over a specific period. Specifically, the former dataset contains 24M locations and time points of 29,000 visitors, with 2 minutes as a time frame. We sampled the trajectory data of 16 to 2048 people in 60 time frames (2 hours) for a data snippet in the training stage. We used the reanalysis of 7 species from eBird as the dataset, containing 176M data points. For each data snippet, we sampled 16 individual birds over 256 to 1024 time frames. As presented in Table 1, Deflex identified multiple formulas, including power law equations [40], Lévy flight patterns [41], and vector navigation laws [42] in human crowds, and the Lévy flight patterns in bird migrations. Besides Deflex, AI Feynman, PySR, and DEAP can also discover the power law distribution and Lévy flight pattern (shown in Table 2). But no other method can identify the vector navigation law. It is worth noting that Deflex also uncovered an undocumented power law distribution, P



(i) ∆xt



α −γ

∝r 0

P j̸=i

I

(j)

(i)

xt −xt

<2.46

.

(2)

This rule may indicate a trend superimposed on the baseline power-law distribution and describe the collective behavior of human crowds: individuals exhibit a higher probability of making a long-distance move as their neighbors do. For the commonly found rules of the power law distribution and Lévy flight pattern, we first compare the NLL loss and RMSE of Deflex with other methods in Fig. 3a. It shows that Deflex’s discovered distributions have consistently high accuracy (both smaller NLL loss and RMSE value) than other capable methods. We also compare their time efficiency for discovering these distributions under the same accuracy (i.e., by setting an identical precision threshold for convergence). Fig. 3b shows that Deflex also leads to the least running time among the compared methods.

8

2.5 Water particles To further validate the performance of Deflex in multiscale systems, we conducted experiments on cross-scale water particles at different temperatures (300K and 273.15K) and scales. In particular, ice particles can be regarded as either clusters of water molecules or larger particles. This experiment was conducted using the same microtube as the pollen experiment. For scale-independent characteristics, Deflex demonstrates similar performance in water particles as on Argon and pollen particles, as shown in Table 1, Figs. 3a and 3b. To validate the cross-scale properties, we conducted coarse-graining of water particles by aggregating groups of water molecules that form lattice structures into macro-particles (i.e., ice crystals). This process was iteratively repeated to construct a hierarchically coarse-grained model of water. Theoretically, as the coarse-graining progresses, the dynamics of ice crystals were expected to transition from a particlebased description to a Langevin dynamics framework. In our analysis, we focused on tracking the short-term dynamics of the water system, typically over a span of 50 time frames. This constraint arises because the coarse-grained model cannot capture the formation or melting of ice crystals, limiting its applicability to scenarios where these processes are not significant. As shown in Fig. 4a, Deflex manages to capture the cross-scale emergence of Langevin dynamics in water particles. In contrast, other methods are limited to observing molecular dynamics at smaller scales and the dynamics of ice crystals at larger scales, failing to bridge the transition across scales. It also shows that Deflex manages to detect the scale invariance of the system at the phase transition point of water (i.e., at 273.15K).

2.6 Fluid dynamics We conducted experiments on fluids (Fig. 3a, b), one of the most typical complex systems in physics [43]. Specifically, we focused on the wake behind a cylinder in a two-dimensional flow. The specific experimental setup is detailed in SI Document S2.2. We built a 1024×1024 mesh for simulation, and each one of 1M grids was regarded as an element, containing the position, pressure, and velocity of the grid as the data fields, and 2048 time frames were collected in total. We collected 9 to 256 adjacent or random grids and 128 time frames as a data snippet for training. To test the generalizability of our approach across different flow regimes and spatial dimensions, we conducted a second fluid dynamics experiment using the Johns Hopkins Turbulence Database (JHTDB) [44]. Unlike the confined cylinder wake flow, the JHTDB provides Direct Numerical Simulation of three-dimensional forced isotropic turbulence, representing a canonical case of fully developed turbulence without boundary effects. The dataset features a 10243 computational grid with a Taylor microscale Reynolds number Reλ ≈ 433. We extracted volumetric subsamples from the database, with each grid point containing the three-dimensional velocity field (u, v, w), pressure p, and their spatial derivatives. Following the same training protocol as the cylinder experiment, we sampled 9 to 256 spatial locations and 128 temporal snapshots as

9

input data. Remarkably, Deflex successfully identified the three-dimensional NavierStokes equation from this turbulent dataset (Table 1, Figs. 3 and 4), capturing the full three-dimensional convective and diffusive terms. The complementary nature of these two experiments—2D laminar wake versus 3D turbulent flow—demonstrates Deflex’s capability to discover governing equations across different flow regimes, spatial dimensions, and Reynolds numbers. As listed in Table 1, Deflex detects the Navier-Stokes equation of fluid dynamics through small-scale observations, which cannot be discovered by other methods. In larger-scale random observations, Deflex extracts the velocity distribution of the fluid within the vortex flow. These formulas correspond to different manifestations of longstanding laws in physics. For complex distribution formulas like velocity distribution, most compared methods are infeasible (shown in Table 2), and only AI Feynman and Deflex manage to discover. Notably, from the discovered Navier-Stokes equation, we can further identify the term corresponding to Newton’s viscosity law, which is also listed among our discoveries in Tables 1 and 2. Meanwhile, SINDy and Operon can recover a simplified form of 2D Navier–Stokes without pressure-field terms (i.e., ∇p), which leads to noticeably larger errors than Deflex, even though the simpler structure can make their discovery faster. Deflex is the only method that can recover the full 3D Navier–Stokes equation in the JHTDB setting. Fig. 3a also shows that the equation-based formulas discovered by Deflex can achieve a similar RMSE as other methods. For the velocity distribution, Deflex’s discovery has a much smaller NLL loss than AI Feynman’s, showing its better accuracy. Fig. 3b, most comparatives, except for SINDy and Operon, become much slower than Deflex with the increase of formula complexity. Among them, AI Feynman is around 10 times slower than Deflex. In addition, we compared the EMD across different coarse-graining scales of Deflex with that of AI Feynman and the renormalization of known formulas. Fig. 4a shows that predictions from Deflex keep the lowest EMD across different scales, indicating the Deflex’s ability to overcome the scale gap issue in complex systems. Meanwhile, other methods with a single setting can only capture the patterns at a single scale, and the EMD increases significantly when the scale changes.

2.7 Comparison on benchmark dataset We finally evaluated Deflex on the public benchmark dataset of Feynman Symbolic Regression Database [7]. This benchmark dataset is generated by sampling from 120 classical formulas spanning different fields of physics, with added noise to simulate real-world conditions. It captures essential characteristics of real-world mathematical formulas, such as symmetry, while also being a classic, well-established set of problems for symbolic regression. This dual characteristic facilitates a feasible performance comparison against other SR-based AI methods. In the experiment, Deflex was compared against Operon, SciMED, gplearn, AI Feynman 2.0, and PySR. We compared their performance by measuring the trade-off between the average RMSE of discovered formulas and the total runtime for the discovery. Note that SINDy was not included as it is sparse regression-based and not designed for symbolic regression tasks, while DEAP was omitted since it works as the backend of gplearn in the current task. 10

As shown in Fig. 4b, Deflex demonstrates quite competitive performance among all compared methods on the benchmark dataset in terms of the accuracy-efficiency tradeoff, although it is not specialized for these traditional, closed-form equations. In particular, Deflex is only slightly outperformed by Operon, which relies on a fixed symbol set and is optimized for runtime. In our experiment, Deflex managed to discover the formula with an RMSE of 0.075 in just 146.2 seconds, which was quite close to Operon with an RMSE of 0.0663 in 135.8 seconds. And Deflex and Operon show significant advantages over other methods. This high efficiency is largely attributed to the pre-training mechanism, which accelerates formula discovery by incorporating mathematical priors (Fig. 5a, d).

2.8 Ablation analysis To validate the contribution of each component in Deflex, we conducted systematic ablation studies across multiple dimensions (Fig. 5). We first examined the impact of pre-training on discovery efficiency by varying the synthetic dataset size from 0 to 106 formulas. Results demonstrate that larger pre-training volumes lead to substantially faster convergence and higher formula recovery rates (Fig. 5a, d). Pre-training with 106 examples achieved over 80% recovery rate within 5,000 steps, while training from scratch required more than 20,000 steps to reach similar performance, confirming that incorporating mathematical priors through pre-training is crucial for efficient exploration of the formula space. Architectural ablations reveal that both network depth and embedding dimensions critically affect the model’s capacity to capture complex system dynamics. Increasing the number of Deflexformer blocks from 1 to 8 progressively improves performance, with NLL loss decreasing from approximately 1.8 to 0.6 (Fig. 5b, e). Similarly, expanding embedding dimensions from 32 to 128 enhances the model’s representational power (Fig. 5c, f). However, performance saturates beyond 6 blocks and 128 dimensions, indicating that these configurations provide sufficient capacity for the tested systems. We also conducted ablation studies on the derivative (gradient) features from the datasets, as detailed in the Supplementary Information, finding that accuracy is mostly preserved but discovery efficiency degrades, and the complicated 3D Navier–Stokes case (JHTDB) does not reliably converge under the same time budget. Moreover, the Deflexformer effectively internalizes physical constraints, autonomously learning energy and momentum conservation manifolds (Fig. 5h). Finally, comparing data sampling strategies demonstrates that Langevin sampling is essential for accurately identifying complex equilibrium distributions (Fig. 5g).

3 Discussion We present an AI-based framework Deflex for scientific rules discovery from complex systems, which can extract multiscale mathematical formulas through neural-networkguided decomposition and symbolic regression. In Deflex, we propose a unified energy model and distribution-based representation of formulas, which enables multiscale probability modeling, while leveraging lambda calculus extends symbolic regression’s capacity to handle massively parallel entity interactions without prior structural 11

assumptions. Deflex first trains a neural network to decompose intricate system interactions into elementary components, then interprets these into compact mathematical formulas using lambda-calculus-based symbolic regression. By breaking down formula learning into tractable sub-problems, this decomposition strategy solves the issue of exponential complexity in traditional symbolic regression. Extensive experiments on six typical complex systems across diverse domains, ranging from particle dynamics and fluid mechanics to human and bird mobility, demonstrate the effectiveness and efficiency of Deflex in multiscale formula discovery. Experimental results demonstrate that Deflex can manage to overcome the scale gap issue by achieving the least differences across scales. Compared with state-of-theart single-scale methods, it can also achieve better performance in both accuracy and computational efficiency in most cases, where the runtime reduction can achieve up to orders of magnitude. In addition to reconstructing established principles (e.g., Langevin equations, energy conservation), Deflex also shows the potential to reveal undocumented system relationships from the observed data. The critical advancement of Deflex lies in the multiscale characterization. As evidenced by the experiments on ice crystals and fluid dynamics, Deflex can capture scale-consistent patterns without renormalization procedures, surpassing traditional multiscale description methods. This advantage stems from the synergy between the invariant modeling of EBMs and the compositional flexibility of lambda calculus, eliminating the prerequisite for explicit scale separation or prior formula knowledge. In contrast, conventional methods require substantial prior specification of the potential formula structures, severely limiting their autonomous discovery capability. Nonetheless, on the benchmark AI Feynman dataset, we observe that Deflex has a lower recovery rate on formulas of quantum transition probability derived from time-dependent perturbation theory [45] (see SI for detailed introduction and the recovery comparison). We guess that this phenomenon exposes a potential limitation of Deflex’s multiscale discovery. That is, when candidate structures of markedly different magnitudes coexist, low-frequency approximations may overwhelm the original high-frequency but complex structures. This suggests an open problem in balancing competing scales so that small yet physically meaningful terms are not suppressed. From a practical application perspective, Deflex requires that the observation of complex systems can be organized into learning samples. In a typical physical setting, each sample combines global context with evolving states of a set of n interacting entities over consecutive time steps. In practice, one can obtain a long continuous record and then construct training samples by extracting spatiotemporal sub-trajectories, allowing the model to learn an energy-form expression that aggregates variablesize interactions. Successful symbolic recovery further demands sufficiently diverse and broad coverage of the system’s state space. This can be done by collecting multiple trajectories under varied initial conditions. If data collection is limited, sampling over multiple windows and entity subsets may help to increase the number of training samples without altering the underlying physics. In real applications, the scientist-in-the-loop approach can also be incorporated into data collection and feature engineering when domain knowledge is needed.

12

While demonstrating significant progress, several limitations in the current implementations still warrant further investigation. First, the computational efficiency of deep formula synthesis remains constrained by excessive invalid candidate generation in symbolic regression. Second, the current framework assumes pre-processed observational data, whereas real-world applications may require integrated solutions for automated data acquisition and hypothesis generation. Therefore, future work may explore improving symbolic search efficiency via hybrid inference strategies combining neural guidance with type-theoretic constraints. Third, regarding interferences of the approximated patterns among multiple scales, a promising direction is to separate scales explicitly in the discovery procedures, including problem formulation, rather than extracting implicitly within a single joint run. Also, it is important to establish theoretical connections between the learned energy landscapes and renormalization group theory. Another practical problem is handling noisy, incomplete observational data to enhance the robustness for real-world complex systems. In this work, we have presented Deflex, a framework for discovering governing equations from multiscale complex systems. Our empirical results demonstrate its ability to identify symbolic rules across different scales and, in some cases, reveal previously undocumented relationships from observational data. The proposed methodology for representing interactions among high-dimensional entities suggests potential applicability to domains beyond those tested in our experiments. Future work could explore the practical meaning of the undocumented formulas with area experts, and Deflex’s further use in areas such as atmospheric dynamics modeling [22] and materials science [46], particularly for systems where strong scale interdependencies pose challenges to conventional methods. We hope that this work serves as a step toward more capable automated scientific discovery in a range of disciplines.

4 Methods We present the details of the Deflex framework for discovering mathematical laws in complex systems, which can automatically learn concise mathematical formulas from the observed data of the systems. Throughout this section, we use bold symbols to denote vectors and matrices (e.g., x for individual samples, θ for model parameters). Specifically, X = {x1 , . . . , xN } represents the observation dataset with N samples, each xi characterizing the dynamics of element i over time. We denote theoretical energy functions as E(·) or E ∗ (·) for the optimal one, and use Eθ (·) to refer to the parameterized energy function approximated by the Deflexformer neural network. The partition function is denoted as Zθ , and probability distributions as p(·) for data distributions and Qθ (·) for model distributions.

4.1 Problem formulation We consider describing a complex system by a set of measurable features x, which normally include the element-level ones s and optionally the global-level ones v, i.e., x = [s; v]. Specifically, the element-level features s are state variables that describe

13

the local, dynamic condition of individual components (or elements) within the system (e.g., a particle’s position and velocity), which are time-varying. The global-level features v may characterize the system from a macroscopic level (e.g., a system-wide temperature) and is considered independent of both element amount and time. In particular, Suppose we can obtain an observation dataset X = [S; V] from the system, where S ∈ RN T Ds denotes the Ds -dimensional state s of N elements within a time period of T , and the optional V ∈ RDv represents the Dv dimensional global information v. Note that, except for the category of global or elementary features, we have no prior knowledge about the relationships among the individual features. Let F denote the set of all possible mathematical formulas, including equations, invariants, and distributions. The behaviors in the system can be described by a set of formulas f ⊆ F on X. Normally, f = {f1 , · · · , fm } (where ∀fi ∈ f , fi ∈ F) with m formulas is expected to satisfy two properties: 1) f should fit the observed data well, i.e., have a large likelihood. Specifically, we can define a certain loss function L(fi , X) to measure how well fi ∈ F fits the datasets X, where lower loss indicates better fit. 2) f should be concise, that is, have low complexity in mathematical formulas, denoted as R(f ). Then, the goal of Deflex is to find the mathematical formulas f ∗ ⊆ F that balance the trade-off between the complexity of the expression of the function and the fitness of the data (in terms of the expected loss over the observed data X). This can be formulated as the Pareto optimal set of patterns f ∗ ⊆ F that satisfies: ∀g ⊆ F,

X f ∈f

L(f, X) ≤ ∗

X g∈g

L(g, X) ⇒

X f ∈f

R(f ) ≤

X g∈g∗

R(g).

(3)

This is essentially a symbolic regression (SR) problem. Existing SR methods explore to solve the above problem by evaluating large-scale candidate expressions. However, facing complex systems with complicated observation data and multiscale rules, existing studies would suffer from three major challenges: • Limited expressibility for higher-order relationships: traditional expression trees struggle to represent mappings and higher-order functions autonomously, often resorting to rigid, pre-specified operators or suffering from explosive tree growth when modeling collective behaviors. • Disjoint representations across scales: existing approaches lack a unified framework to characterize and learn both deterministic laws (invariants) and stochastic patterns (distributions) that emerge at different scales. • Inefficiency in complex search spaces: the search space for complex nonlinear formulas grows exponentially with the number of variables, making standard search algorithms computationally prohibitive for high-dimensional systems.

4.2 Distribution-based formula representation We first unify formulas in the form of probability distributions, which can be represented as energy functions that map the observation data X (with both element-level state information S and the optional global context V) to scalar-valued energy f1 (X), . . . , fm (X). Then the search space F of Deflex is unified as the set of energy functions. 14

4.2.1 Unified formula representation with distributions. For a unified representation, we can regard any formula f ∈ F as a probability distribution p(x) [30]. Specifically, a collective formula f = {f1 , · · · , fm }, describing the behaviors of a single system, can be seen as a set p = {p1 , · · · , pm } of m independent distributions pi (x) over x, i.e., p ≡ f . Among them, equations and invariants are treated as distributions with very narrow peaks, which indicates their tiny errors or uncertainty. The original distribution-based formulas are treated as distributions with wider peaks. This representation unifies different forms of formula expressions, thus simplifying the formula discovery. Furthermore, any two independent distributions pi (x), pj (x) ∈ p can be equivalently described as a single mixed distribution p′{i,j} (x) of them. Ultimately, p with m independent distributions can also be seen as a set p′ = {p′{1,...,m} } with a single mixed distribution p′{1,...,m} (x) of p1 (x), · · · , pm (x). Then, we also say p′ = p ≡ f . This aligns well with the characteristics of different scales in complex systems and naturally resolves the scale gap issue in multiscale pattern discovery. That is, probability distributions on a large scale can describe the behaviors on a small scale. In other words, narrow distributions correspond to precise small-scale patterns, whereas wider distributions correspond to large-scale patterns. Based on the unified representation, solving the problem in Eq. (3) requires discovering an optimal distribution p∗ (x) that best fits (with the maximum likelihood of) pdata (x), which denotes the underlying data distribution of the observation dataset.

4.2.2 Characterizing distributions with energy functions. The unified representation can tackle the scale gap issue in multiscale systems. However, without the ability to normalize probability density functions, SR cannot directly discover the distribution-based formulas. Therefore, we further characterize distributions as energy functions under the energy-based models (EBM [47–50]) framework, which are less restrictive in functional forms (i.e., without requiring a normalized probability). EBMs are probabilistic models inspired by the Boltzmann distribution in statistical physics. Let p(x) denote a distribution over features x. Then, p(x) can be characterized as an energy function E(x), which indicates different energy values over the features x. A probability distribution function (PDF) p(x) given by an EBM is: p(x) =

1 exp(−E(x)) Z

(4)

R where Z = exp(−E(x))dx is the normalization factor (i.e., partition function) and ensures the integral of the distribution equals 1. It can avoid complex normalization calculations and volume transformations for complex and multi-dimensional variables. Therefore, finding the optimal distribution p∗ (x) is essentially equivalent to finding the optimal energy function E ∗ (x). In our framework, we use the Deflexformer neural network Eθ (x) to approximate E ∗ (x), which can be achieved by training on observation data.

15

4.3 Deflexpressor: SR subsystem based on lambda calculus Traditional tree-based SR relies on first-order expression trees and manually specified higher-order primitives, which makes it difficult to learn relationships such as summation or mapping over massive interacting elements in complex systems. To this end, we present Deflexpressor, an SR system based on the extended lambda calculus, which supports various higher-order expressions while ensuring valid types and variable length of elementary inputs. In particular, Deflexpressor leverages a strict rule-based generator together with an efficient neural network-based one to generate a large number of valid candidate lambda expressions for mathematical relationships, which are then evaluated to identify Pareto-optimal expressions in terms of accuracy and conciseness.

4.3.1 Expressing relationships by lambda calculus Lambda calculus is a formal system in mathematical logic and computer science for expressing computation based on function abstraction and application [51]. A lambda term is built from three kinds of expressions: (i) Variable x denotes a placeholder for a value or parameter; (ii) Abstraction λx.M defines a function with parameter x and body M ; (iii) Application M N applies function M to argument N . Although there is no built-in data, the lambda calculus can encode data and control structures (e.g., recursion) as functions via Church encoding. The lambda calculus is Turing-complete and can simulate any computation based on symbol replacement. Specifically, for any expression made of variables, symbol λ, and parentheses, we can find a pattern (λx.M )N and replace it with M where every x is replaced by N . By repeating this process until no such pattern remains, we can obtain the final expression as the result. The basic grammar rules are provided in Section 3.1 in SI. Since the lambda calculus is Turing-complete and natively supports function abstraction and recursion, Deflexpressor can express operations such as summation and mapping directly rather than as fixed templates. Specifically, operations over collections can be written as recursively defined functions or higher-order maps, enabling the system to learn how to aggregate interactions from data without requiring preP specified structures. For instance, we consider a common formula E = i̸=j V (rij ) that describes the aggregation of a variable number of particle interactions. Here, V (rij ) denotes a function that calculates the pairwise interaction potential based on the distance rij between two particles indexed by i and j. Traditional SR either requires hard-coding symbols or extremely deep trees to express such a relationship. In contrast, suppose that x = [rij ]i̸=j denotes the list that enumerates these distances under any fixed ordering, Deflexpressor can naturally write the aggregation as the following Lambda expression.

16



 Fix λf λxs λa.

Application (i)

| | | |

{z

Variables

}



if Null(xs) then a else f (Tail(xs)) a + V (Head(xs)) {z } | {z } |

x0.

Application (ii)

|

{z

}

Function body describing the whole relationship (i and ii) among f, xs, a

{z

Abstraction defining a function G(f, xs, a) over variables f, xs, a

}

{z

Fixed-point simulating a recursive function F(xs, a) based on G(f, xs, a)

{z

Application denoting applying recursive function F(xs, a) over parameters x and 0

} }

(5) Eq. (5) defines a recursive aggregation function F via a fixed-point construction and applies it to a variable-length list x and an initial accumulator 0. The underbraces decompose the term: the innermost abstraction defines a one-step update rule G(f, xs, a) whose body implements (i) a base case that returns a when Null(xs) and (ii) an update-and-continue step that calls f on Tail(xs) with the updated accumulator a + V (Head(xs)). The operator Fix(·) turns this update rule into a self-referential recursion F (formal definition in the SI), and application associates to the left, so F x 0 means (F x) 0. See SI for a step-by-step expansion and an explicit walk-through of the bracketed decomposition. This demonstrates that the system constructs the algorithm of summation from scratch, covering cases where predefined summation templates are inapplicable or insufficient. Besides strong expressiveness, another merit of the lambda calculus is that its relatively unified structure enables intuitive comparison of the complexity of formulas. In particular, we can directly use the expression length as the complexity metric R(f ) (in Eq. 3) in Pareto selection, enabling an explicit trade-off between data fitness and conciseness.

4.3.2 Extended lambda calculus with types and arrays The pure lambda calculus is Turing complete, but it has two limitations in our complex system setting. The first is that it may generate meaningless expressions and even non-terminating reductions, which would lead to the low efficiency of the symbolic regression process in Deflexpressor. The second is that data structures like lists or arrays are often encoded as nested higher-order functions, which leads to significantly high complexity in data traversal. To this end, in Deflexpressor, we adopt an extended lambda calculus with types and arrays. A type specifies what kind of value a sub-expression denotes (e.g., a real number, a vector, or a function), and therefore constrains how expressions can be composed. Specifically, we enforce the following type constraints. (i) Type compliance in the applications. For example, an application (M N ) is valid only when M denotes a function whose input type matches the type of N . Instead, if M : R denotes a realvalued number, then (M N ) will be rejected. (ii) Type requirements for the operands of operators. For instance, expression M = 1 + (λx. x) will be identified as invalid because λx. x (i.e., f (x) = x) belongs to a function type, whereas + expects two operands with the real-valued type. Type constraints on input-output format. For 17

example, many tasks require a real-valued function on an n-dimensional state, i.e., f : Rn → R. In such cases, a well-typed function termP can be used, whereas a constant will be rejected in Deflexpressor. For instance, E = i̸=j V (rij ) requires a mapping from a variable-length list of pair distances x = [rij ]i̸=j to a scalar energy. Eq. (5) has exactly this form (a function that consumes x and returns E), thus being a typevalid expression. In Deflexpressor, we implement a type inference engine based on the classical type inference algorithm W [52] (see SI for details), which can automatically infer the type of each candidate expression, and then keep only type-valid candidates. Variable-sized data structures, such as lists, arrays, and vectors, are essential to represent various data in complex systems. However, encoding such structures as functions in pure lambda calculus would lead to extra workload of searching the encoding itself rather than the relationship of encoded data. To adapt the framework for scientific discovery, we predefine array Array[τ ] as a built-in primitive data type, and provide O(1)-time access and basic operators (e.g., indexing and size). For instance, P we treat the pair distances rij in E = i̸=j V (rij ) as an array x = [rij ]i̸=j . This allows the aggregation to be written in the same form as Eq. (5), letting the search focus on learning V (·) and the aggregation structure rather than on data representation details (formal definitions in SI).

4.3.3 Hybrid expression generation with rules and NNs Exploring the vast space of lambda expressions presents distinct challenges in terms of efficiency and coverage. Firstly, rule-based mutations (e.g., random subtree replacement) explore the space diffusively and often generate low-quality expressions. To boost search efficiency, we incorporate a causal autoregressive neural network based on the Transformer Decoder architecture (see SI for architecture details). This network is trained offline via self-supervised learning on a manually-prepared corpus of valid expression trees and their rule-based mutation pairs. By learning the probabilistic priors of valid structural transformations, the network can generate syntactically correct and type-consistent lambda expressions either from scratch or by mutating existing ones. Note that, the produced expressions are parsed and type-checked to ensure validity. For example, a subtree mutation can replace the local update term V (Head(xs)) in Eq. (5) with V ′ (Head(xs)) without changing the recursion structure. β-reduction further simplifies candidates by evaluating function applications, e.g., ((λx. x) y) → y (formal rules in SI). Secondly, relying solely on a neural network can lead to mode collapse, where the search gets stuck in familiar patterns and misses novel solutions. To expand search coverage, we implement a hybrid sampling strategy with simulated annealing. The process transitions from predominantly rule-based generation (for broad exploration) to neural-guided generation (for exploitation). Third, expressions after multiple mutations often become bloated with redundancy. To maintain compact and interpretable forms, we enforce lambda-calculus-specific mutations like β-reduction and common sub-expression elimination in the rule-based component (see SI for mutation details).

18

4.3.4 Initialization with global aggregation priors While Deflexpressor can generate expressions of any computable functions via the lambda calculus, standard random initialization tends to produce expression trees that operate on a fixed number of arguments. In such a case, it would be combinatorially expensive to evolve a structure that iteratively computes over a variable-length input from scratch. To accelerate this process without imposing hard constraints, we employ a heuristic initialization strategy that generates the seed expressions performing various aggregations over variable-length inputs. Specifically, these seed expressions are written in the form λx. Reduce(kernel, x), which denotes the recursive procedure of updating an aggregation state s by traversing a variable-length input list x = [x1 , . . . , xN ] as follows  x1 , i = 1 s← (6) kernel (s, xi ) , i > 1 where kernel(s, xi ) is a learnable binary operator that updates the aggregation state s using an element xi of the list x. In particular, the kernel is a randomly generated expression tree that defines the interaction logic between the aggregation P state s and each data element xi of a variable-length list. For example, in E = i̸=j V (rij ), we can define a kernel as kernel(s, r) = s + V (r) and take x = [rij ]i̸=j as the list. Then the function Reduce(kernel, x) can yield the desired aggregation, i.e., it stands for computing E by iteratively applying s ← s + V (r) over all pairwise distances rij in x (with the recursive procedures realized via Fix(·) as in Eq. (5)). This initialization serves as a soft inductive bias that prioritizes aggregation operations over massive elements, encouraging SR to focus immediately on discovering local interaction laws (encoded in the kernel) while fitting the data. Meanwhile, since λx. Reduce(kernel, x) is also a recursive lambda expression that can be implemented with the same recursion primitive Fix(·) in Eq. (5), it remains fully mutable under subsequent lambda calculus mutations. Subsequent evolutionary steps (e.g., subtree mutation, β-reduction) can alter the kernel, break the recursion, or discard this “Reduce-like” structure entirely if a simpler non-recursive law fits the data well.

4.4 Deflexformer: energy-based NNs with decomposable blocks The optimal energy function E ∗ (·) can theoretically be inferred via SR over the observation dataset. However, this can be computationally expensive and even infeasible for complex expressions. Generally, any complex E(·) can be expressed as the composition of a series of (e.g., k) base functions E i (·), e.g., E(·) = E 1 ◦ E 2 ◦ . . . ◦ E k (·). So, we propose to approximate E(·) by training a deep neural network Deflexformer Eθ (·) (parameterized by θ) with k decomposable blocks E i (·), each approximating a base function E i (·) (where i = 1, . . . , k). The decomposable structure enables penetrating each block E i (·) to derive explicit expressions of E i (·), making SR feasible to discover complex expressions of E(·).

19

4.4.1 Architecture of Deflexformer The Deflexformer network predicts the energy value for a given dataset using the aligned state sequences of elements and global context as inputs. The architecture comprises an input embedding layer, k stacked Deflexformer blocks, and an output projection module. The input layer embeds raw observations X into a high-dimensional feature space via Fourier multiscale embedding (see SI), producing element-level representations P0 ∈ Rn×t×Dmodel and global representations G0 ∈ Rc×Dmodel , where n is the number of elements, t is the number of time frames, c is the number of global tokens, and Dmodel = 64 is the model dimension. For some observations of a complex system, X may only include the element-level state information S without the global information v. In such a case, we set a dummy global feature v and randomly generate a dummy global context dataset V by drawing some random noise over v (see SI for details). Note that, designed for a much more complicated SR problem over observation data of complex systems with multiple scales, Deflexformer can also be directly used for classical SR tasks (e.g., applied to the AI Feynman dataset in Section 2.7), which correspond to a single scale without considering element-level dynamics. Here, the SR dataset is input as the global context, and the element-level representation is left empty at the input layer. Each Deflexformer block E i (·) (Fig. 2a) transforms representations (Pi−1 , Gi−1 ) into (Pi , Gi ), approximating a base mathematical transformation E i (·) in the energy decomposition E(·) = E 1 ◦ E 2 ◦ . . . ◦ E k (·). Each block comprises two sub-modules operating on the concatenated input of dimension (c + nt) × Dmodel : (1) A pointwise transformation network implemented as a four-layer fully connected network with architecture (Dmodel , 2Dmodel , 2Dmodel , Dmodel ) and ReLU activations, modeling element-wise evolution (e.g., the dynamics of an individual particle). (2) A multi-head cross-attention network with h = 8 heads for capturing the interactions among particles. The attention operates in two sequential stages (Fig. 2b): spatial mixing performs attention over the element dimension to capture inter-element interactions, while temporal mixing applies causal attention over the time dimension to model temporal dependencies. By factorizing attention into spatial and temporal axes, computational complexity is reduced from O(nt) to O(n + t). Residual connections and layer normalization are applied after each sub-module (see SI for detailed specifications). The output module projects the final global representation Gk back to a scalar energy via inverse Fourier embedding e = E −1 (Gk ) (see SI). For example, when modeling particle energy E = 12 mi vi2 + U (xi ), early blocks may capture kinetic energy contributions while later blocks model potential interactions, with their composition yielding the total energy.

4.4.2 Pre-training the Deflexformer block Directly training Deflexformer E(·) from scratch requires a vast amount of observation data. Yet, its blocks E i (·) are normally simple functions, which have a much smaller search space. To this end, we adopt a common two-stage training framework of pretraining and post-training in other transformer-based deep networks. In particular, a

20

single Deflexformer block is first pre-trained with a large amount of synthetic data based on some specified relationships. Each Deflexformer block E i (·) represents a basic mathematical transformation, similar to some base functions (e.g., sin(·)) with a simple expression form. Therefore, we can pre-train the Deflexformer block on a synthetic dataset sampled from some known base functions. Specifically, we first call Deflexpressor to generate a pool of expressions for functions f (·). Then, based on these functions, we can sample numerous data pairs (x, y) to synthesize a pre-training dataset D, where x is drawn from a Gaussian distribution and y = f (x). Here, pre-training aims to fit a block model E ′ (·) via supervised learning on the synthetic dataset. So the loss function is defined simply as the MSE between the predictions with labels, given the P|D| features of data pairs in the synthetic dataset. That is, LE ′ = i=1 MSE(y, E ′ (x)) (see SI for training details).

4.4.3 Post-training the full Deflexformer Then, the full Deflexformer, composed of multiple copies of the pre-trained block and an output module, is post-trained on the observation data. Note that Deflexformer Eθ (·) is an energy-based model for approximating the optimal energy function E(·). So, its PDF Qθ (·) should fit pdata (·) with a maximized likelihood L(Eθ , X). Here, our likelihood function adopts the common expected log-likelihood over the data distribution in the EBM framework as L(Eθ , X) = Ex∼pdata [log Qθ (x)] = E[log

1 exp(−Eθ (x))], Zθ

(7)

where θ denotes the network parameters, and Zθ is the partition function dependent on θ. Then, we can train Eθ to maximize (7). A common training methodology is to define a loss function as the negative log-likelihood (NLL): Lθ = −L(Eθ , X) =

1 XN (Eθ (xi ) + log(Zθ )), i=1 N

(8)

and adopt the Markov Chain Monte Carlo (MCMC) method with Langevin sampling [53, 54]. While computing the partition function Zθ itself is intractable, its gradient ∇θ log Zθ can be efficiently estimated via MCMC sampling, which is sufficient for gradient-based optimization (see SI for detailed methodology and parameter setups). Crucially, this approach discovers probability distributions directly from raw data points via likelihood maximization, rather than fitting predefined probability density functions (PDFs) to pre-computed statistics.

4.4.4 Extracting intermediate representations by inference The post-trained Deflexformer Eθ approximates an energy function E(·) while the post-trained block E i corresponds to the base functions E i (·) within E(·). As Pi and Gi denote the element-level and global-level representations output by the E i (P0 and G0 correspond to the representations of raw data X), there should be Pi ||Gi = E i (Pi−1 ||Gi−1 ) where Pi ||Gi means the concatenation of representations. 21

More generally, there is Pi ||Gi = E i−l ◦ · · · ◦ E i (Pi−l−1 ||Gi−l−1 ). Therefore, to identify any relationships E i−l ◦· · ·◦E i , we can perform SR on a large number of input-output representation pairs (Pi−l−1 ||Gi−l−1 , Pi ||Gi ) by inferring over the post-trained Deflexformer Eθ . Here, we can feed each sample in the observation data into Deflexformer to conduct forward computation over a block array E i−l ◦ · · · ◦ E i of multiple sequential blocks, producing the corresponding intermediate representations. After multiple inferences, we can collect a dataset to perform SR, thus identifying any relationship represented by an array with an arbitrary number of sequential blocks. However, the number of collected samples may be constrained by that of the observation dataset, thereby limiting SR accuracy. Fortunately, as an energy model, Deflexformer inherently supports generating more data samples with a similar distribution to the original observation data. Specifically, it can work by using Langevin sampling from random noise (see SI for sampling details).

4.4.5 Hierarchical SR from intermediate representations SR for a smaller block array with fewer Deflexformer blocks can reduce the computation complexity, but may converge to fragmented expressions with no meaning or oversimplified with larger decomposition error. Then their composition would lead to significantly biased energy functions. Nevertheless, these incomplete/inaccurate expressions can serve as base expressions to accelerate the SR for a larger block array, which has smaller decomposition errors. To this end, we propose a hierarchical SR strategy for recovering the whole energy function. Here, denote Ê i−l ◦. . .◦Ê k ← SR(Pi−l−1 ||Gi−l−1 , Pk ||Gk ) as SR for the l consecutive blocks indexed by k − l, k − l + 1, . . . , k. Let k = 1, 2, . . . , K denote the block index. We perform the SRs by building a binary tree T with the depth of ⌈log K⌉. For each node v in T , we denote T (v) as the sub-tree with node v as its root. In T , each node v represents an SR process over the consecutive blocks indexed by the leaf nodes of the subtree T (v). The data samples (input-output representation pairs) for SRv are constructed by extracting the pair of the input representation of T (v)’s leftmost leaf node and the corresponding output representation of T (v)’s rightmost leaf node. In particular, for the leaf layer, the base expression set for SR is user-specified. For other layers in the tree, the base expression set for any node v is constructed by taking the union of the SR-derived expression set of its two children. Up to the root, we can gradually refine the SR results and finally obtain an expression set {Ê1 , . . . , Êm } for the entire Deflexformer that approximates the complete energy function E(·). Although this incurs O(K log K) SRs, it can greatly reduce the searching space of SR for complex relationships (see SI for detailed algorithm).

Data availability The data generated or processed in this study are available in the open-source project repository at https://github.com/yhqjohn/deflex. For the simulation-based projects, the repository provides the generated datasets for rare gas particle motion, pollen-in-water particle motion, cross-scale water particles, and two-dimensional cylinder-wake fluid dynamics. The original third-party datasets used in this study 22

are publicly available from their source repositories or publications: the Zurich Carnival human-mobility dataset at https://doi.org/10.4108/icst.urb-iot.2014.257190; the eBird Status and Trends Data Version 2022 at https://doi.org/10.2173/ebirdst.2022; the JHTDB dataset at https://turbulence.pha.jhu.edu/; and the Feynman Symbolic Regression Database at https://space.mit.edu/home/tegmark/aifeynman.html, with associated publication DOI https://doi.org/10.1126/sciadv.aay2631. No additional access restrictions apply to the data made available through the project repository.

Code availability The code used in the study is available in the open-source development repository at https://github.com/yhqjohn/deflex.

References [1] Wigner, E.P.: The unreasonable effectiveness of mathematics in the natural sciences. In: Mathematics and Science, pp. 291–306. World Scientific, Singapore (1990) [2] Schmidt, M., Lipson, H.: Distilling free-form natural laws from experimental data. Science 324(5923), 81–85 (2009) https://doi.org/10.1126/science.1165893 [3] Irrgang, C., Boers, N., Sonnewald, M., Barnes, E.A., Kadow, C., Staneva, J., Saynisch-Wagner, J.: Towards neural earth system modelling by integrating artificial intelligence in earth system science. Nature Machine Intelligence 3(8), 667–674 (2021) [4] Qi, D., Majda, A.J.: Using machine learning to predict extreme events in complex systems. Proceedings of the National Academy of Sciences 117(1), 52–59 (2020) https://doi.org/10.1073/pnas.1917285117 [5] Varadi, M., Anyango, S., Deshpande, M., Nair, S., Natassia, C., Yordanova, G., Yuan, D., Stroe, O., Wood, G., Laydon, A., et al.: Alphafold protein structure database: massively expanding the structural coverage of protein-sequence space with high-accuracy models. Nucleic acids research 50(D1), 439–444 (2022) [6] Udrescu, S.-M., Tan, A., Feng, J., Neto, O., Wu, T., Tegmark, M.: Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Neural Information Processing Systems 33, 4860–4871 (2020) [7] Udrescu, S.-M., Tegmark, M.: Ai feynman: A physics-inspired method for symbolic regression. Science Advances 6(16), 2631 (2020) https://doi.org/10.1126/ sciadv.aay2631 [8] Bellomo, N., Dogbe, C.: On the modeling of traffic and crowds: A survey of models, speculations, and perspectives. SIAM review 53(3), 409–463 (2011) 23

[9] Bar-Yam, Y.: General features of complex systems. Encyclopedia of Life Support Systems (EOLSS), UNESCO, EOLSS Publishers, Oxford, UK 1 (2002) [10] Council, N.R., et al.: A National Strategy for Advancing Climate Modeling. MIT Press, Cambridge, MA (2012) [11] Fortunato, S., Bergstrom, C.T., Börner, K., Evans, J.A., Helbing, D., Milojević, S., Petersen, A.M., Radicchi, F., Sinatra, R., Uzzi, B., et al.: Science of science. Science 359(6379), 0185 (2018) [12] Pomeau, Y.: The long and winding road. Nature Physics 12(3), 198–199 (2016) [13] Cranmer, M., Sanchez-Gonzalez, A., Battaglia, P., Xu, R., Cranmer, K., Spergel, D., Ho, S.: Discovering symbolic models from deep learning with inductive biases. In: Proceedings of the 34th International Conference on Neural Information Processing Systems. NIPS’20. Curran Associates Inc., Red Hook, NY, USA (2020) [14] Rudy, S.H., Brunton, S.L., Proctor, J.L., Kutz, J.N.: Data-driven discovery of partial differential equations. Science Advances 3(4), 1602614 (2017) https://doi. org/10.1126/sciadv.1602614 [15] Gao, T.-T., Barzel, B., Yan, G.: Learning interpretable dynamics of stochastic complex systems from experimental data. Nature Communications 15(1), 6029 (2024) [16] Van Doren, B.M., Horton, K.G.: A continental system for forecasting bird migration. Science 361(6407), 1115–1118 (2018) [17] Hasselmann, K.: Stochastic climate models part i. theory. tellus 28(6), 473–485 (1976) [18] Vicsek, T.: Collective Motion. Elsevier (2012) [19] Kaneko, K., Tsuda, I.: Complex Systems: Chaos and Beyond: Chaos and Beyond: A Constructive Approach with Applications in Life Sciences. Springer, Berlin/Heidelberg, Germany (2001) [20] Weinan, E.: Principles of Multiscale Modeling. Cambridge University Press, Cambridge, UK (2011) [21] Pun, G.P., Batra, R., Ramprasad, R., Mishin, Y.: Physically informed artificial neural networks for atomistic modeling of materials. Nature communications 10(1), 2339 (2019) [22] Bi, K., Xie, L., Zhang, H., Chen, X., Gu, X., Tian, Q.: Accurate medium-range global weather forecasting with 3d neural networks. Nature 619(7970), 533–538 (2023) 24

[23] Romera-Paredes, B., Barekatain, M., Novikov, A., Balog, M., Kumar, M.P., Dupont, E., Ruiz, F.J., Ellenberg, J.S., Wang, P., Fawzi, O., et al.: Mathematical discoveries from program search with large language models. Nature 625(7995), 468–475 (2024) [24] Brunton, S.L., Proctor, J.L., Kutz, J.N.: Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences 113(15), 3932–3937 (2016) https://doi.org/10. 1073/pnas.1517384113 [25] Koza, J.R.: Genetic programming as a means for programming computers by natural selection. Statistics and computing 4, 87–112 (1994) [26] Keren, L.S., Liberzon, A., Lazebnik, T.: A computational framework for physics-informed symbolic regression with straightforward integration of domain knowledge. Scientific Reports 13(1), 1249 (2023) [27] Sun, F., Liu, Y., Wang, J.-X., Sun, H.: Symbolic physics learner: Discovering governing equations via monte carlo tree search. In: The Eleventh International Conference on Learning Representations (2023) [28] Jiang, N., Xue, Y.: Racing control variable genetic programming for symbolic regression. In: Proceedings of the 38th AAAI Conference on Artificial Intelligence (AAAI 2024). AAAI Press, Vancouver, Canada (2024). https://doi.org/10.1609/ aaai.v38i11.29187 [29] Mundhenk, T., Landajuela, M., Glatt, R., Santiago, C.P., faissol, D., Petersen, B.K.: Symbolic regression via deep reinforcement learning enhanced genetic programming seeding. In: Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P.S., Vaughan, J.W. (eds.) Advances in Neural Information Processing Systems, vol. 34, pp. 24912–24923. Curran Associates, Inc., Second Virtual Conference (2021) [30] Siegenfeld, A.F., Bar-Yam, Y.: An introduction to complex systems science and its applications. Complexity 2020, 1–16 (2020) [31] Burlacu, B., Kronberger, G., Kommenda, M.: Operon c++: An efficient genetic programming framework for symbolic regression. In: Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion. GECCO ’20, pp. 1562–1570. Association for Computing Machinery, New York, NY, USA (2020). https://doi.org/10.1145/3377929.3398099 [32] Stephens, T., Kemenade, H., Rai, A., Price, B., Watson, C., Bell, I., McDermott, J., Vos, N., Niculae, S., Di Caprio, U., sun ao: trevorstephens/gplearn. https://github.com/trevorstephens/gplearn (2025). https://github.com/ trevorstephens/gplearn [33] Cranmer, M.: Interpretable Machine Learning for Science with PySR and 25

SymbolicRegression.jl. arXiv (2023). https://doi.org/10.48550/arXiv.2305.01582 [34] Fortin, F.-A., De Rainville, F.-M., Gardner, M.-A., Parizeau, M., Gagné, C.: DEAP: Evolutionary algorithms made easy. Journal of Machine Learning Research 13, 2171–2175 (2012) [35] Villegas, P., Gili, T., Caldarelli, G., Gabrielli, A.: Laplacian renormalization group for heterogeneous networks. Nature Physics 19(3), 445–450 (2023) [36] Villani, C., et al.: Optimal Transport: Old and New. Grundlehren der mathematischen Wissenschaften, vol. 338. Springer, Berlin Heidelberg (2009) [37] Langevin, P., et al.: Sur la théorie du mouvement brownien. CR Acad. Sci. Paris 146(530-533), 530 (1908) [38] Blanke, U., Guldener, R., Feese, S., Tröster, G.: Crowdsourced pedestrian map construction for short-term city-scale events. In: Proceedings of the First International Conference on IoT in Urban Space. URB-IOT ’14, pp. 25–31. ICST (Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering), Brussels, BEL (2014). https://doi.org/10.4108/icst.urb-iot.2014. 257190 [39] Fink, D., Auer, T., Johnston, A., Strimas-Mackey, M., Ligocki, S., Robinson, O., Hochachka, W., Jaromczyk, L., Crowlye, C., Dunham, K., Stillman, A., Davies, I., Rodewald, A., Ruiz-Gutierrez, V., Wood, C.: eBird Status and Trends, Data Version: 2022; Released: 2023. Cornell Lab of Ornithology (2023). https://doi. org/10.2173/ebirdst.2022 [40] Gonzalez, M.C., Hidalgo, C.A., Barabasi, A.-L.: Understanding individual human mobility patterns. nature 453(7196), 779–782 (2008) [41] Schläpfer, M., Dong, L., O’Keeffe, K., Santi, P., Szell, M., Salat, H., Anklesaria, S., Vazifeh, M., Ratti, C., West, G.B.: The universal visitation law of human mobility. Nature 593(7860), 522–527 (2021) [42] Bongiorno, C., Zhou, Y., Kryven, M., Theurel, D., Rizzo, A., Santi, P., Tenenbaum, J., Ratti, C.: Vector-based pedestrian navigation in cities. Nature Computational Science 1(10), 678–685 (2021) [43] Solera-Rico, A., Sanmiguel Vila, C., Gómez-López, M., Wang, Y., Almashjary, A., Dawson, S.T., Vinuesa, R.: β-variational autoencoders and transformers for reduced-order modelling of fluid flows. Nature Communications 15(1), 1361 (2024) [44] Li, Y., Perlman, E., Wan, M., Yang, Y., Meneveau, C., Burns, R., Chen, S., Szalay, A., Eyink, G.: A public turbulence database cluster and applications to study lagrangian evolution of velocity increments in turbulence. Journal of

26

Turbulence (9), 31 (2008) [45] Sakurai, J.J., Napolitano, J.: Modern Quantum Mechanics, 2nd edn. Cambridge University Press, Cambridge, UK (2017) [46] Benzi, R., Paladin, G., Parisi, G., Vulpiani, A.: On the multifractal nature of fully developed turbulence and chaotic systems. Journal of Physics A: Mathematical and General 17(18), 3521 (1984) [47] Song, Y., Kingma, D.P.: How to train your energy-based models. arXiv preprint arXiv:2101.03288 (2021) [48] Wenliang, L., Sutherland, D.J., Strathmann, H., Gretton, A.: Learning deep kernels for exponential family densities. In: International Conference on Machine Learning, pp. 6737–6746 (2019). PMLR [49] Song, Y., Ermon, S.: Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems 32 (2019) [50] Roney, J.P., Ovchinnikov, S.: State-of-the-art estimation of protein model accuracy using alphafold. Physical Review Letters 129(23), 238101 (2022) [51] Barendregt, H.P.: The Lambda Calculus: Its Syntax and Semantics, Revised 3rd edn., pp. 23–315460. North-Holland, Amsterdam (1984) [52] Milner, R.: A theory of type polymorphism in programming. Journal of Computer and System Sciences 17(3), 348–375 (1978) https://doi.org/10.1016/ 0022-0000(78)90014-4 [53] Parisi, G.: Correlation functions and computer simulations. Nuclear Physics B 180(3), 378–384 (1981) [54] Grenander, U., Miller, M.I.: Representations of knowledge in complex systems. Journal of the Royal Statistical Society: Series B (Methodological) 56(4), 549–581 (1994)

Funding This work was supported in part by the National Key Research and Development Program of China (2022YFA1004100 [S.Y.]); in part by the National Natural Science Foundation of China (U21A6005 [S.Y.], 62172329 [X.R.]).

Author contributions S.Y. and X.R. conceived and supervised the project. H.Y., S.Y., and X.R. designed the study and algorithms. H.Y. implemented the code. H.Y., S.Y., X.R., and C.Z.

27

performed the experiments. H.Y., S.Y., and X.R. wrote the manuscript. All authors discussed the results and reviewed the manuscript.

Competing interests The authors declare no competing interests.

28

Generation Deflexpressor

𝒙𝒙

λ 𝒚𝒚 + −× 𝒙𝒙 𝒂𝒂

87.9

𝒙𝒙𝟏𝟏

42.6

𝒙𝒙𝟐𝟐

96.1

35.9

43.1

14.6

49.7

34.5

Generating equations

𝒚𝒚

Intermediate representation

Representation dataset

Deflexformer Block N

63.2

(27.3 26.1), (75.9 94.2), (95.1 51.6)

Deflexformer Block N

Deflexformer Block 2 Deflexformer Block 1 Input

(69.4 43.2), (31.4 22.1), (82.1 73.1)

SR

Training

(67.1 66.3), (15.5 84.8) (29.2 33.4), (51.6 52.7)

SR

�(⋅) + C

SR

𝒏𝒏

� 𝒊𝒊=𝟎𝟎

+ 𝑪𝑪

𝟏𝟏 𝒎𝒎 𝒗𝒗𝟐𝟐 + 𝑼𝑼(𝒙𝒙𝒊𝒊 ) 𝟐𝟐 𝒊𝒊 𝒊𝒊

𝟏𝟏 𝒎𝒎 𝒗𝒗𝟐𝟐 + 𝑼𝑼(𝒙𝒙𝒊𝒊 ) 𝟐𝟐 𝒊𝒊 𝒊𝒊

SR Generating 𝒔𝒔1 + 𝒔𝒔3

𝒔𝒔[1] + 𝒔𝒔[3]

Mutation & evaluation

𝒙𝒙

λ 𝒚𝒚 + −

𝒙𝒙

×

𝒂𝒂

λ λ x λ x . λ x . x

𝒏𝒏

�(⋅) + (⋅)1 (Fix (λ𝑓𝑓𝑓𝑓𝑓𝑓𝒔𝒔. null 𝒔𝒔 ? 𝑎𝑎 ∶ 𝑎𝑎 + f (tail 𝒔𝒔) g (head 𝒔𝒔)))𝒔𝒔[1] 𝒔𝒔 𝒊𝒊=𝟎𝟎

Mutation & evaluation

(⋅)𝟐𝟐 , 𝑼𝑼(⋅)

𝒏𝒏

�(⋅) + C 𝒊𝒊=𝟎𝟎

(67.1 66.3), (15.5

Inference

𝒊𝒊=𝟎𝟎

Observation dataset

𝒊𝒊=𝟎𝟎

SR

𝟏𝟏 𝒎𝒎 𝒗𝒗𝟐𝟐 + 𝑼𝑼(𝒙𝒙𝒊𝒊 ) + 𝑪𝑪 𝟐𝟐 𝒊𝒊 𝒊𝒊 … 𝟏𝟏 𝒎𝒎 𝒗𝒗𝟐𝟐 + 𝑼𝑼(𝒙𝒙𝒊𝒊 ) 𝟐𝟐 𝒊𝒊 𝒊𝒊 𝒗𝒗𝟐𝟐𝒊𝒊 𝑼𝑼(𝒙𝒙𝒊𝒊 )

�

…

𝒏𝒏

𝟏𝟏 𝒂𝒂 𝒃𝒃 + 𝒄𝒄i 𝟐𝟐 i i

Latent function ≈

𝒏𝒏

(⋅)𝟐𝟐 , 𝑼𝑼(⋅)

Input

e Hierarchical symbolic regression

(27.3 26.1), (75.9 94.2), (95.1 51.6)

𝟏𝟏 𝒂𝒂 𝒃𝒃 + 𝒄𝒄i 𝟐𝟐 i i

Deflexformer Block 1

Pre-trained block

…

𝒊𝒊=𝟎𝟎

Deflexformer Block 2

𝜽𝜽,𝒊𝒊 𝒙𝒙

SR

𝒏𝒏

�(⋅) + C

…

Deflexformer Block 𝓝𝓝

Output

Training

Output

Pre-Training

Sampling data

Post-training

Deflexformer

…

𝒚𝒚 = 𝒙𝒙𝟏𝟏 + 𝒔𝒔𝒔𝒔𝒔𝒔 𝒙𝒙𝟐𝟐 𝒚𝒚 = 𝟐𝟐𝒙𝒙𝟐𝟐𝟏𝟏 − 𝒍𝒍𝒍𝒍𝒍𝒍 𝒙𝒙𝟐𝟐 𝒚𝒚 = 𝒆𝒆𝒙𝒙𝟏𝟏 + 𝒙𝒙𝟑𝟑𝟐𝟐

Deflexformer

Synthetic dataset

λ λ x λ x . λ x . x

Formulas

d Inference

c

b Pre-training

…

a

84.8) …

λ𝑔𝑔(Fix (λ𝑓𝑓𝑓𝑓𝑓𝑓𝒔𝒔. null 𝒔𝒔 ? 𝑎𝑎 ∶ 𝑎𝑎 + f (tail 𝒔𝒔) g (head 𝒔𝒔))) 𝐶𝐶𝒔𝒔

Observation dataset

Fig. 1 Illustration of the workflow in Deflex. Deflex utilizes two subsystems: Deflexpressor and Deflexformer, where the former is a symbolic regression subsystem based on lambda calculus and the latter is an energy-based DNN subsystem. The overall process flows sequentially from a through b, c, d, and finally to e. a, Deflexpressor generates a large number of expressions, from which it samples a synthetic dataset. b, the synthetic dataset is used to pre-train a Deflexformer block, a Transformer-like neural network. c, the pre-trained block is repeatedly cascaded multiple times to build Deflexformer, which is then trained to learn latent formulas on the observation dataset in the target complex system. d, the trained Deflexformer is fed with the observation dataset and some random data for inference, which generates a large number of data representations output by all blocks. e, the multi-level representation dataset from all Deflexformer blocks is processed by Deflexpressor through hierarchical symbolic regression. The SR framework iteratively generates expression candidates, evaluates them across different representation levels, and refines them through mutation to discover the desired formulas as both mathematical expressions and executable lambdacalculus code.

29

a Deflexformer Block 1

Point-wise data �0 (� × � × �)

Block 1

Deflexformer Block k

Output

Deflexformer block× � +

Q KV

�1

… Block 2

Point-wise transform

+

Global context �0 (� × �)

Deflexformer Block 2

…

Input

�1

…

Block �

Global context �� (� × �) Point-wise data �� (� × � × �)

b Global context (� × �)

Point-wise data (� × � × �)

Elementwise self attention

+

+

Framewise causal self attention +

+

Global context (� × �) Point-wise data (� × � × �)

Fig. 2 Illustration of the neural network architecture of Deflex. a illustrates the overall model architecture of Deflexformer. Arrows show the information flow among the various components. The intermediate expressions are separated into element-wise and global variables and are processed by different neural networks. b illustrates the self-attention network for a Deflexformer block. The pointwise data consists of t time frames, each containing n elements with d dimensions. The point-wise data is first processed by the element-wise residual self-attention with the fusion of global variables and then processed by the frame-wise residual causal self-attention with the fusion of global variables.

30

Table 1 Discovered expressions by Deflex in comparison with the ground truths Rules Energy conservation Momentum conservation Maxwell-Boltzmann distribution (MBD) MBD with potential Continuity of quantity Newton viscosity law

Velocity distribution Navier-Stokes equation Langevin Lévy flight (Human) Lévy flight (Avian) Power law Crowd power law

Ground truths and discovered by Deflex True Discovered True Discovered True

E = v 2 · m/kB T E = 0.992v 2 · m/kB T P P mi v i = C mi v i = C mv 2 E = 2k T

Discovered

E = 0.511mv k T

True

mv E = 2k + T

Discovered True Discovered

0.97 ∥x⃗j − x⃗i ∥−12 − 1.89 ∥x⃗j − x⃗i ∥−6 E = 0.506mv + kB T j̸=i E = ∥∇ · v∥ × ∞ E = ∥∇ ·  v∥ × 1026 

True

B

2

B 2

B

τij = −µ

P

2

j̸=i

∥x⃗j − x⃗i ∥−12 − 2 ∥x⃗j − x⃗i ∥−6

P

∂u ∂ui + ∂xj ∂xj i



∂u ∂ui + ∂xj ∂xj i P



Discovered

τij = −1.037µ

True Discovered True Discovered True Discovered True Discovered True Discovered True Discovered True Discovered

E = 21 mv 2 + (δv) E = 0.501mv 2 + (δv)
3.15 + u · ∇u + ∇p − µ∇2 u∥ × ∞ E = ∥ρ ∂u ∂t
∂u E = ∥ρ ∂t + u · ∇u + ∇p − µ∇2 u∥ × 1.62 × 1021 E = ∥4.468 × 1016 mv̇ + v∥2 E = ∥4.438 × 1016 mv̇ + 1.003v∥2 E = D1 log ∆x E = 2.37 log ∆x E = D2 log ∆x E = 1.627 log ∆x E = α log l + γl E = 1.56 log l + 0.00027l Undocumented
P (j) (i) E = 1.63 log ri − 0.037 log ri I xt − x t < 2.46 j̸=i

E denotes the energy value of the corresponding energy function. kB (Boltzmann constant), δ, ρ, P , D1 , D2 , α, and γ are unprovided constants in the experiments. All other symbols (including temperature T in the argon system) denote the provided features of observation data in the experiments. Detailed explanations of these features are described in SI.

31

Table 2 Comparison of methods on discovery capability Rules Energy conservation Momentum conservation MBD MBD with potential Continuity of quantity Newton viscosity law Velocity distribution Navier-Stokes equation Langevin Lévy flight (Human) Lévy flight (Avian) Power law

Deflex

Operon

SciMED

SINDy

gplearn

AI Feynman

PySR

DEAP

✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓

✓ ✓ ✓ N.A. ✓ ✓ ✗ ✓1 ✗ ✗ ✗ ✗

✓ ✓ ✓ N.A. ✓ ✓ ✗ ✗ ✗ ✗ ✗ ✗

✓ ✓ N.A. N.A. ✓ ✓ N.A. ✓1 N.A. N.A. N.A. N.A.

✓ ✓ ✓ ✗ ✓ ✓ ✗ ✗ ✗ ✗ ✗ ✗

✓ ✓ ✓ N.A. ✓ ✓ ✓ ✗ ✗ ✓ ✓ ✓

✓ ✓ ✓ N.A. ✓ ✓ ✗ ✗ ✗ ✓ ✓ ✓

✓ ✓ ✓ ✗ ✓ ✓ ✗ ✗ ✗ ✓ ✓ ✓

Operon and SINDy can recover a 2D Navier–Stokes form (cylinder wake) but fail in the 3D JHTDB setting; moreover, their recovered 2D forms do not include pressure-field terms (i.e., ∇p). ✓ denotes a successful discovery, ✗ denotes a failure, and N.A. denotes an unsupported task (i.e., not applicable).

32

a

Me thod

De flex

SciMED

gple arn

SINDy

Ope ron

AI Feynman

PySR

DE AP

Argon Particles Pollen Particles Human Mobility Birds Migration Water Particles Fluid Dynamics

RMSE

100 10-1 10-2 10-3

NLL(log)

10 5 0 -5

b

Argon Particles Pollen Particles Human Mobility Birds Migration Water Particles Fluid Dynamics

Time (s)

104 103 102 101 100

Form ula

Fig. 3 Experimental results of Deflex in terms of accuracy (a) and efficiency (b). a compares the RMSE and NLL loss of the discovered formulas by Deflex and other methods in the experiments of particle motions, human and bird collective movement, water particles, and fluid dynamics. The abbreviations represent the Maxwell-Boltzmann distribution (MBD), energy conservation law (EC), momentum conservation (MC), Langevin dynamics (LD), continuity of quantity (CQ), Newton viscosity law (NV), velocity distribution in fluid dynamics(VD), Navier-Stokes equation (NS), power law in human mobility (Power), Lévy flight rule (Lévy), vector navigation law (Vec), and crowd power law, respectively. Among them, the formulas in the form of distributions are marked with ∗. The absence of the method indicates that the corresponding metrics are not available for the method. b depicts the running time and usability of Deflex compared with other methods in different experiments. The horizontal axis represents the running time of the methods and the absence of the method indicates it is infeasible to discover the formula. Only Deflex can discover all the formulas in the experiments, and the time consumption of Deflex is significantly lower than that of most other methods for complex formulas. In panels a and b, error bars show mean ± standard deviation.

33

a

Me thod

Molecular 300K

SR-I

LD

MBD

ND

VD

NS

b 0.4

Fluid Dynamics

0.4

1

3 10 Scale

0.2 0.1 0.0

1

3 10 Scale

0.6

RMSE

0.2 0.1

0.3

EMD/Scale

0.3

EMD/Scale

EMD/Scale

SR-D

Molecular Transition

0.4

0.0

De flex

0.4 0.2 0.0

Method De flex SciMED gple arn Ope ron AI Feynman PySR

0.3 0.2 0.1

1

10 Scale

100

200

300 400 Time (s)

Fig. 4 Experimental results of Deflex in terms of cross-scale performance (a) and benchmark performance (b). a shows the EMD values of known formulas and formulas discovered by Deflex and other methods across different coarse-graining scales for water particles’ motion and fluid dynamics, where “SR-I” and “SR-D” represent the symbolic regression methods with invariant and distributional form settings, respectively. b shows the running time and relative mean square error for Deflex and other methods on the public benchmark dataset of Feynman Symbolic Regression Database.

34

a

b

d

e

f

g

c

h

Fig. 5 Ablation studies and component analysis of Deflex. a, d, Impact of pre-training volume on formula discovery performance (recovery rate and relative speed) and NLL trend. Larger pretraining volumes significantly accelerate discovery and improve accuracy. b, e, Performance metrics and NLL loss with varying Deflexformer block numbers. Increasing depth improves representation power. c, f, Effect of embedding dimensions on performance. g, NLL loss comparison across different sampling schemes (Random, Random+Real, Langevin). Langevin sampling yields the best fitness. h, Verification of conservation laws (Energy and Momentum) learned by the Deflexformer, showing close alignment between sampled states (black dots) and ideal laws (red lines). In panels a–c, the error bars show mean value ± standard deviation.

35