Vertical Symbolic Regression via Deep Policy Gradient

1 Introduction

Exciting progress has been made to accelerate scientific discovery using Artificial Intelligence (AI) [Langley et al., 1987; Kulkarni and Simon, 1988; Wang et al., 2023]. Symbolic regression, as an important task in AI-driven scientific discovery, distills physics models in the form of symbolic equations from experiment data [Schmidt and Lipson, 2009]. Notable progress in symbolic regression encompasses diverse methodologies, includes search-based methods [Lenat, 1977], genetic programming [Schmidt and Lipson, 2009; Virgolin et al., 2019], Monte Carlo tree search [Todorovski and Dzeroski, 1997; Sun et al., 2023; Kamienny et al., 2023], and deep reinforcement learning [Petersen et al., 2021; Mundhenk et al., 2021].

Recently, Vertical Symbolic Regression (VSR) [Jiang and Xue, 2023; Jiang et al., 2023] has been proposed to expedite the discovery of symbolic equations with many independent variables. Unlike previous approaches, VSR reduces the search spaces following a vertical path - it extends from reduced-form equations involving a subset of independent variables to full-fledged ones, adding one variable into the equation at a time. Figure 1 gives an example discovery of Joule's law Q ∝ I 2RT [Joule, 1843], where Q is heat, I is current, R is resistance, and T is time. VSR first holds I and R as constants and finds Q ∝ T. Next, I is introduced with targeted experiments that study the effect of I on Q. Such rounds repeat until all factors are considered. Compared with the horizontal paths, that modeling all independent variables simultaneously, vertical discovery can be significantly cheaper because the search spaces of the first few steps are exponentially smaller than the full search space.

Meanwhile, deep learning, especially deep policy gradient [Petersen et al., 2021], has greatly boosted the performance of symbolic regression approaches. However, VSR was implemented using genetic programming. Although we hypothesize deep neural nets should boost VSR, directly integrating deep neural nets to predict the symbolic equation tree in each vertical expansion step encounters difficulties. This will result in (1) difficulty passing gradients from trees to deep neural nets and (2) complications concatenating deep networks for predictions in each vertical expansion step. We provide a detailed analysis of the difficulty in Appendix A.

In this work, we propose Vertical Symbolic Regression using Deep Policy Gradient (VSR-DPG). We demonstrate that VSR-DPG can recover ground-truth equations involving 50 variables, which is beyond both deep reinforcement learningbased approaches and the previous VSR variant (i.e., VSRGP). Our VSR-DPG solves the above difficulty based on the following key idea: each symbolic expression can be treated as repeated applications of grammar rules. Hence, discovering the best symbolic equations in the space of all candidate expressions is viewed as the sequential decision-making of predicting the optimal sequence of grammar rules.

In Figure 1(right), the expansion from C1T to C2I 2T is to replace constant C1 with a sub-expression C2I 2 . We define a context-free grammar on symbolic expression, denoting the rules that certain constants can be replaced with other variables, constants, or sub-expressions. All candidate expressions that are compatible with C1T can be generated by repetitively applying the defined grammar rules in different order. VSR-DPG generates many sub-expressions (including C2I 2 ) by sequentially sampling rules from the Recurrent Neural Networks (RNN). A vertical discovery path is built on top of this sequential decision-making process of reducedform symbolic expressions. The RNN is trained to generate expansions that lead to high fitness scores on the dataset. In this regard, we train the RNN to maximize a policy gradient objective, similar to that proposed in Petersen et al.. The main difference is the RNN in our VSR-DPG predicts the next rules in the vertical discovery space, while the model from Petersen et al. predicts the next math token in the expression tree in the horizontal discovery space.

In experiments, we consider several challenging multivariable datasets of algebraic equations and of ordinary differential equations in material science and biology. (1) In Table 1, our VSR-DPG attains the smallest median NMSE values in 7 out of 8 datasets, against current popular baselines including VSR-GP. The main reason is deep networks offer more parameters than GP, which can better adapt to different datasets and sample higher-quality expressions from the networks. (2) Further analysis on the best-discovered equation (in Table 3) shows that VSR-DPG uncovers up to 50% of the exact governing equations with 5 input variables, where the baselines only attain 0%. (3) In Table 2, our VSR-DPG can find high-quality expressions on datasets with up to 50 variables, because of the vertical discovery idea. (4) On discovery of ordinary differential equations in Table 4, our VSR-DPG also improves over current baselines. 1

1 The code is at https://github.com/jiangnanhugo/VSR-DPG.

2 Preliminaries

Symbolic Regression aims to discover governing equations from the experimental data. An example of such mathematical expression is Joule's first law: Q = I 2RT, which quantifies the amount of heat Q generated when electric current I flows through a conductor with resistance R for time T . Formally, a mathematical expression ϕ connects a set of input variables x and a set of constant coefficients c by mathematical operators. The possible mathematical operators include addition, subtraction, multiplication, division, trigonometric functions, etc. The meaning of these mathematical expressions follows their standard arithmetic definition.

Given a dataset D = {(xi, yi)Ni=1 |xi ∈ Rn , yi ∈ R} with N samples, symbolic regression searches for the optimal expression ϕ∗ , such that ϕ∗ (xi, c) ≈ yi. From an optimization perspective, ϕ∗ minimizes the averaged loss on the dataset:

Here Π is the set of all candidate mathematical expressions; c denotes the constant coefficients in the expression; ℓ(·, ·) denotes a loss function that measures the difference between the output of the candidate expression ϕ(xi, c) and the ground truth yi. The set of all possible expressions, i.e., the hypothesis space Π, can be exponentially large. As a result, finding the optimal expression is challenging and is shown to be NPhard [Virgolin and Pissis, 2022].

Deep Policy Gradient for Symbolic Regression. Recently, a line of work proposes the use of deep reinforcement learning (RL) for searching the governing equations [Petersen et al., 2021; Abolafia et al., 2018; Mundhenk et al., 2021]. They represent expressions as binary trees, where the interior nodes of the tree correspond to mathematical operators and leaf nodes of the tree correspond to variables or constants. The key idea is to model the search of different expressions, as a sequential decision-making process for the preorder traversal sequence of the expression trees, using an RL algorithm. A reward function is defined to measure how well a predicted expression can fit the dataset. The deep RNN is used as the RL agent for predicting the next possible node in the expression tree at every step of decision-making. The parameters of the RNN are trained using policy gradient.

Control Variable Experiment studies the relationship between a few input variables and the output in the regression problem, with the remaining input variables fixed to be the same [Lehman et al., 2004]. In the controlled setting, the ground-truth equation behaves the same after setting those controlled variables as constants, which is noted as the reduced-form equation. For example, the ground-truth equation ϕ = x1 × x3 - x2 × x4 in Figure 2(a) is reduced to x1 × C1 - C2 when controlling x2, x3, x4. Figure 2(b,c) presents other reduced-form equations when the control variables are changed. For the corresponding dataset D, the controlled variables are fixed to one value and the remaining variables are randomly assigned. See Figure 2(a,b,c) for example datasets generated under different controlling variables.

Please refer to https://arxiv.org/abs/2402.00254 for the Appendix.

Vertical Symbolic Regression starts by finding a symbolic equation involving a small subset of the input variables and iteratively expands the discovered expression by introducing more variables. VSR relies on the control variable experiments introduced above. VSR-GP was the first implementation of vertical symbolic regression using genetic programming (GP) [Jiang and Xue, 2023]. To fit an expression of n variables, VSR-GP initially only allows variable x1 to vary and controls the values of all the rest of the variables. Using GP as a subroutine, VSR-GP finds a pool of expressions {ϕ1, . . . , ϕm} which best fit the data from this controlled experiment. Notice {ϕ1, . . . , ϕm} are restricted to contain only one free variable x1 and m is the pool size. A small error indicates ϕi is close to the ground truth reduced to the one free variable and thus is marked unmutable by the genetic operations in the following rounds. In the 2nd round, VSR-GP adds a second free variable x2 and starts fitting {ϕ′1 , . . . , ϕ′m} using the data from control variable experiments involving the two free variables x1, x2. After n rounds, the expressions in the VSR-GP pool consider all n variables. Overall, VSR expedites the discovery process because the first few rounds of VSR are significantly cheaper than the traditional horizontal discovery process, which searches for optimal expression involving all input variables.

3 Methodology

Motivation. The prior work of VSR-GP uses genetic programming to edit the expression tree. GP is not allowed to edit internal nodes in the best-discovered expression trees from the previous vertical discovery step, to ensure later genetic operations do not delete the prior knowledge on the governing equation. However, this idea cannot be easily integrated with deep RL-based symbolic regressors.

We explored using deep neural networks in vertical symbolic regression to predict the symbolic equation tree in each vertical expansion step, as detailed in Appendix A. However, we encountered two main issues: (1) difficulty in passing gradients from trees to deep neural nets. We need different constraints to embed into the output of RNN to ensure the whole model follows the vertical discovery idea. (2) complexity in concatenating deep networks for predictions at each vertical expansion step. The underlying issue is the expression tree representation.

Our solution is to adopt a new representation for expressions. We extend the context-free grammar definition for symbolic expressions, where a sequence of grammar rules uniquely corresponds to an expression. We view the prediction of symbolic expressions as a sequential decision-making process, selecting grammar rules step-by-step. The RNN predicts grammar rules instead of nodes in the expression tree. The best-discovered reduced-form equation is converted into the start symbol in the grammar, ensuring that the predicted expression from our RNN is always compatible with the prior knowledge of the governing equation. This reduces the hypothesis space and accelerates scientific discovery, as other non-reducible expressions are never sampled from the RNN.

Main Pipeline. Figure 1 shows our deep vertical symbolic regression (VSR-DPG) pipeline. The core idea is to construct increasingly complex symbolic expressions involving more input variables, based on control variable experiments with To fit an expression of n variables, we first hold n-1 variables as constant and allow only one variable to vary. We aim to find the best expression ϕ1, that fits the data in this controlled experiment. We use the RNN as the RL agent to search for the best possible expression, which is achieved by using the RNN to sample sequences of grammar rules defined for symbolic equations. Every sequence of rules is converted into an expression, where the constants in the expression are fitted with the dataset. The parameters of the RNN model will be trained through the policy gradient objective. The expression with the best fitness score is returned as the prediction of the RNN. A visualized process is in Figure 1(a, e).

Following the idea in VSR, the next step is to classify each

constant as either a summary constant or a standalone con-

stant. (1) A constant not relevant to any controlled variables is

considered standalone and is preserved in subsequent rounds.

(2) A constant that is a sub-expression involving controlled

variables is labeled as a summary type and will be expanded

in subsequent rounds. In our implementation, a constant with

high variance in fitted values across multiple control variable

experiments is classified as a summary type; otherwise, it is

fewer controlled variables. classified as standalone.

Assuming we find the correct reduced-from equation ϕ1 af-

ter several learning epochs. To ensure VSR-DPG does not forget this discovered knowledge of the first round, we want all the expressions to be discovered in the following rounds can be reduced to ϕ1. Therefore, we construct ϕ1 as the start symbol for the following round by replacing every summary constant in ϕ1 as a non-terminal symbol in the grammar (noted as A), indicating a sub-expression. For example, in Figure 2(a), both of them are summary constants so the 1st round bestpredicted expression x1 × C1 - C2 is converted to the start symbol x1 ×A-A for the second round.

In the second round, VSR-DPG adds one more free variable and starts fitting ϕ2 using the data from control variable experiments involving two free variables. Similar to the first round, we restrict our search to sub-expressions with the second variable. It is achieved by limiting the output vocabulary of the RNN model. In Figure 2(f), the RNN model finds an expression ϕ2 = x1×C3-(x2×C4). This means that the RL agent learns to expand C1 with another constant C3 and C2 with sub-expression (x2 × C4), based on the best-discovered result of the 1st round ϕ1 = x1 ×C1 -C2.

VSR-DPG introduces one free variable at a time and expands the equation learned in the previous round to include this newly added variable. This process continues until all the variables are considered. After n rounds, we return the equations with the best fitness scores. The predicted equation will be evaluated on data with no variable controlled. See the steps in Figure 1(b,c,d) for a visual demonstration. We summarize the whole process of VSR-DPG in Algorithm 1 in the appendix. The major difference of this approach from most state-of-the-art approaches is that those baselines learn to find the expressions in the full hypothesis space with all input variables, from a fixed dataset collected before training. Our VSR-DPG accelerates the discovery process, because of the small size of the reduced hypothesis space, i.e., the set of candidate expressions involving only a few variables. The task is much easier than fitting the expression in the full hypothesis space involving all input variables.

3.1 Expression Represented as Grammar Rules

We propose to represent symbolic expression by extending the existing context-free grammar [Todorovski and Dzeroski, 1997]. A context-free grammar is represented by a tuple (V,Σ,R, S), where V is a set of non-terminal symbols, Σ is a set of terminal symbols, R is a set of production rules and S is a start symbol. In our formulation, (1) Σ is the set of input variables and constants {x1, . . . , xn, const}. (2) The set of non-terminal symbols V represents sub-expressions, like {A}. Here A is a placeholder symbol. (3) The set of production rules R represents mathematical operations such as addition, subtraction, multiplication, and division. That is {A → (A + A), A → (A - A), A → A × A, A → A ÷ A}, where → indicates that the left-hand side is replaced with the right-hand side. (4) The start symbol S is extended to be A, x1 × A - A, or other symbols constructed from the bestpredicted expression under the controlled variables.

Starting with the start symbol, successively applying the grammar rules in different orders results in different expressions. Each rule expands the first non-terminal symbol in the start symbol. An expression with only terminal symbols is a valid mathematical expression, whereas an expression with a mixture of non-terminal and terminal symbols is invalid. The expression can also be represented as a binary tree. We chose the grammar representation over the binary tree representation because it simplifies the vertical symbolic regression process, avoiding gradient passing and heavy engineering issues.

Figure 3 presents two examples of constructing the expression ϕ from the start symbol using a given sequence of grammar rules. In Figure 3(a), we first use the subtraction rule A → (A - A), indicating that the symbol A in expression ϕ = A is expanded to ϕ = (A - A). By repeatedly using grammar rules to expand the non-terminal symbols, we eventually arrive at our desired expression ϕ = x1 × C1 - C2. In Figure 3(b), the start symbol is x1 × A - A. The rule A → const replaces the first non-terminal symbol with a constant C1, resulting in x1 ×C1 -A. Finally, we obtain a valid expression x1 × C1 - C2. Figure 3 provides another example conversion with the start symbol x1 ×A-A.

3.2 Expression Sampling fromRecurrent Network

In our vertical symbolic regression setting, the input and output are in the set of grammar rules that cover each input variable, constants, and math operations. We create an embedding layer for the input, noted as Embd. For each input rule r ∈ R, its d-dimensional embedding vector is Embd(r) ∈ Rd .

Sampling Procedure. The RNN module samples an expression by sampling the sequence of grammar rules in a sequential decision-making process. Denote the sampled sequence of rules as τ = (τ1, τ2, . . .). Initially, the RNN takes in the start symbol τ1 = S and computes the first step hidden state vector h1. At t-th time step, RNN uses the predicted output from the previous step as the input of current step τt. RNN computes its hidden state vector ht using the embedding vector of input token τt and the previous time-step hidden state vector ht-1. The linear layer and softmax function are applied to emit a categorical distribution p(τt+1|τt, ht) over every token in the output vocabulary, which represents the probability of the next possible rule in the half-completed expression pθ(τt+1|τt, . . . , τ1). The RNN samples one token from the categorical distribution τt+1 ∼ p(τt+1|τt, ht) as the prediction of the next possible rule. To conclude, the computational pipeline at the t-th step is shown below:

The weight matrix W ∈ Rd×|R| and bias vector b ∈ Rd are the parameters of the linear layer. The last row in Equation 2 is the softmax layer. The sampled rule rt+1 will be the input for the t + 1-th step. We denote θ as the total parameters of the embedding layer, the RNN, and the linear layer.

After L steps, we obtain the sequence τ = (τ1, . . . , τL) with probability pθ(τ) = ∏L-1 t=1 pθ(τt+1|τ1, . . . , τt). We convert this sequence into an expression by following the procedure described in Section 3.1. If we reach the end of the sequence while there are still non-terminal symbols in the converted expression, we randomly add rules containing only terminal symbols to complete the expression. Conversely, if we obtain a valid expression before the sequence ends, we disregard the remaining elements of the sequence and return the valid expression.

Policy Gradient-based Training. We follow the reinforcement learning formulation to train the parameters of the RNN module [Wierstra et al., 2010]. The sampled rules before the current step t, i.e., (τ1, . . . , τt), is viewed as the state of the t-th step for the RL agent. Those rules in the output vocabulary are the available actions for the RL agent. In the formulated decision-making process, the RNN takes in the current state and outputs a distribution over next-step possible actions. The objective of the RL agent is to learn to pick the optimal sequences of grammar rules to maximize the expected rewards. Denote the converted expression from τ as ϕ. A typical reward function is defined from the fitness score of the expression reward(τ ) = 1/(1+NMSE(ϕ)). The objective that maximizes the expected reward from the RNN model is defined as J(θ) = Eτ∼pθ(τ)(reward(τ)), where pθ(τ) is the probability of sampling sequence τ from the RNN.

In the next step, the gradient with respect to the objective ∇θJ(θ) needs to be estimated. We follow the classic REINFORCE policy gradient algorithm [Williams, 1992]. We first sample several times from the RNN module and obtain N sequences (τ 1 , . . . , τ N ), an unbiased estimation of the gradient of the objective is computed as ∇θJ(θ) ≈ 1N ∑Ni=1 reward(τ i )∇θ log pθ(τ i ). The parameters of the deep network are updated by the gradient descent algorithm with the estimated policy gradient value. In the literature, several practical tricks are proposed to reduce the estimation variance of the policy gradient. A common choice is to subtract a baseline function b from the reward, as long as the baseline is not a function of the sample batch of expressions. Our implementation adopts this trick and the detailed derivation is presented in Appendix B.

Throughout the whole training process, the expression with optimal fitness score from all the sampled expressions is used as the prediction of VSR-DPG at the current round.

Start Symbol Construction in Vertical Discovery. Given the best-predicted equation ϕ and controlled variables xc, the following step is to construct the start symbol of the next rounds. This operation ensures all the future expressions can be reduced to any previously discovered equation thus all the discovered knowledge is remembered. It expedites the discovery of symbolic expression since other expressions that cannot be reduced to ϕ will never be sampled from the RNN. It requires first classifying the type of every constant in the expression into stand-alone or summary type, through multitrail control variable experiments. Then we replace each summary constant with a placeholder symbol (i.e., "A") indicating a sub-expression containing controlled variables.

Following the procedure proposed in [Jiang and Xue, 2023], we first query K data batches (D1, . . . , DK) with the same controlled variables xc. The controlled variables take the same value within each batch while taking different values across data batches. We fit open constants in the candidate expression ϕ with each data batch by the gradient-based optimizer, like BFGS [Fletcher, 2000]. We obtain multiple fitness scores (o1, . . . , oK) and multiple solutions to open constants (c1, . . . , cK). By examining the outcomes of K-trials control variable experiments, we have: (1) Consistent close-tozero fitness scores imply the fitted expression is close to the ground-truth equation in the reduced form. That is ok ≤ ε for all 1 ≤ k ≤ K, where ε is the threshold for the fitness scores. (2) Conditioning on the result in case (1), the j-th open constant is a standalone constant when the empirical variance of its fitted values across K trials is less than a threshold ε′ . In practice, if the best-predicted expression by the RNN module is not consistently close to zero, then all the constants in the expression are summary constants. Finally, the start symbol is obtained by replacing every summary constant with the symbol "A" according to our grammar.

4 RelatedWork

Recently AI has been highlighted to enable scientific discoveries in diverse domains [Kirkpatrick et al., 2021; Jumper et al., 2021; Wang et al., 2023]. Early work in this domain focuses on learning logic (symbolic) representations [Bradley et al., 2001]. Recently, there has been extensive research on learning algebraic equations from data [Petersen et al., 2021; Mundhenk et al., 2021] and learning differential equations from data [Dzeroski and Todorovski, 1995; Brunton et al., 2016; Wu and Tegmark, 2019; Zhang and Lin, 2018; Iten et al., 2020; Cranmer et al., 2020; Raissi et al., 2020; Raissi et al., 2019; Liu and Tegmark, 2021; Xue et al., 2021; Chen et al., 2018]. In this domain, a line of works develops robots that automatically refine the hypothesis space, some with human interactions [Langley et al., 1987; Vald´es-P´erez, 1994; King et al., 2004; King et al., 2009]. These works are related to ours because they also actively explore the hypothesis spaces, while they are in biology and chemistry.

Existing works on multivariate regression predominantly rely on pre-trained encoder-decoder structures with extensive training datasets [Biggio et al., 2021], and larger-scale deep models [Kamienny et al., 2022]. Our VSR-DPG method is tailored to solve multivariate symbolic regression.

The idea of using a control variable experiment tightly connects to the BACON system [Langley, 1977; Langley, 1979; Langley et al., 1981; King et al., 2004; King et al., 2009; Cerrato et al., 2023]. While their methods are rule-based systems due to historical limitations, our approach leverages contemporary deep recurrent neural networks.

Our method is also related to symbolic regression methods using probabilistic context-free grammar [Todorovski and Dzeroski, 1997; Brence et al., 2021; Gec et al., 2022]. While they use a fixed probability to sample rules, we use a deep neural network to learn the probability distribution.

Our VSR-DPG is also tightly connected to deep symbolic regression [Petersen et al., 2021; Mundhenk et al., 2021]. We both use deep recurrent networks to predict a sequence of tokens that can be composed into a symbolic expression. However, their method predicts the preorder traversal sequence for the expression tree while our method predicts the sequence of production rules for the expression.

5 Experiment

5.1 Regression on Algebraic Equations

Experiment Settings. For the dataset on algebraic expressions, we consider the 8 groups of expressions from the Trigonometric dataset [Jiang and Xue, 2023], where each group contains 10 randomly sampled expressions. In terms of baselines, we consider (1) evolutionary algorithms: Genetic Programming (GP), Eureqa [Dubc´akov´a, 2011] and Control Variable Genetic Programming (VSR-GP) [Jiang and Xue, 2023]. (2) deep reinforcement learning: Priority queue training (PQT) [Abolafia et al., 2018], Vanilla Policy Gradient (VPG) [Williams, 1992], Deep Symbolic Regression (DSR) [Petersen et al., 2021], and Neural-Guided Genetic Programming Population Seeding (GPMeld) [Mundhenk et al., 2021]. (3) Monte Carlo Tree Search: Symbolic Physics Learner (SPL) [Sun et al., 2023], (4) pretrained Transformer: end-to-end Transformer (E2ETransformer) [Kamienny et al., 2022]. In terms of evaluation metrics, we use the normalized mean-squared-error (NMSE) of the best-predicted expression by each algorithm, on a separately-generated testing dataset. We report the median instead of the mean value. Symbolic regression belongs to combinatorial optimization problems, which have no mean values and are sensitive to outliers. The detailed settings are in Appendix C.2.

Goodness-of-fit Comparison. We consider our VSR-DPG against several challenging datasets involving multiple variables. In Table 1, we report the median NMSE on the selected algebraic datasets. Our VSR-DPG attains the smallest median NMSE values in 7 out of 8 datasets, against a line of current popular baselines including the original VSR-GP. The main reason is deep networks offer many more parameters than the GP algorithm, which can better adapt to different datasets and sample higher-quality expressions from the deep networks.

Extended Large-scale Comparison. In the real world, scientists may collect all available variables that are more than needed into symbolic regression, where only part of the inputs will be included in the ground-truth expression. We randomly pick 5 variables from all the n variables and replace the appeared variable in expressions confined as (5, 5, 5) in Table 1. In Table 2, we collect the median NMSE values on this large-scale dataset setting. Our VSR-DPG scales well because it first detects all the contributing inputs using the control variable experiments. Notice that those baselines that are easily timeout in this setting are excluded for comparison.

Exact Recovery Comparison. We compare if each learning algorithm finds the exact equation, the result of which is collected in Table 3. The discovered equation by each algorithm is further collected in Appendix D.1. We can observe that our VSR-DPG has a higher rate of recovering the ground-truth expressions compared to baselines. This is because our method first use a control variable experiment to pick what are the contributing variables to the data and what are not.

5.2 Regression on Ordinary Differential Equations

Task Definition. The temporal evolution of the dynamic system is modeled by the time derivatives of the state variables. Let x be the n-dimensional vector of state variables, and ˙x is the vector of their time derivatives. The ordinary differential equation (ODE) is of the form ˙x = ϕ(x, c), where constant vector c ∈ Rm are parameters of the dynamic system. Following the definition of symbolic regression on ODE [Gec et al., 2022; Sun et al., 2023], given a trajectory dataset of state variable and its time derivatives {(x(ti), ˙x(ti))}Ni=1 , the symbolic regression task is to predict the best expression ϕ(x, c) that minimizes the average loss on trajectory data. Other formulations of this problem [d'Ascoli et al., 2024] assume we have no access to its time derivatives, i.e., {(ti, x(ti))}Ni=1 .

Experiment Setting. We consider recent popular baselines for differential equations, including (1) SINDy [Brunton et al., 2016], (2) ODEFormer [d'Ascoli et al., 2024], (3) Symbolic Physics Learner (SPL) [Sun et al., 2023]. (4) Probabilistic grammar for equation discovery (ProGED) [Brence et al., 2021]. In terms of the dataset, we consider the Lorenz Attractor with n = 3 variables, Magnetohydrodynamic (MHD) turbulence with n = 6 variables, and Glycolysis Oscillation

-20

Lorenz MHD Glycolysis

SPL 100% 50% 14.2%

SINDy 100% 0% 0%

ProGED 0% 0% 0%

ODEFormer 0% 0% NA

VSR-DPG (ours) 100% 100% 87%

with n = 7 variables. All of them are collected from [Brunton et al., 2016]. To evaluate whether the algorithm identifies the ground-truth expression, we use the Accuracy metric based on the coefficient of determination (R2 ). The detailed experiment configurations are in Appendix C.4.

Result Analysis. The results are summarized in Table 4. Our proposed VSR-DPG discovers a set of differential expressions with much higher quality than the considered baselines. We further provide a visual understanding of the proposed VSR-DPG method in Figure 4. The data of our VSR-DPG are drawn from the intersection of the mesh plane and the curve on the Lorenz attractor. In comparison, the current baselines draw data by picking a random trajectory or many random points on the curve. We notice the ODEFormer is pre-trained on differential equations up to two variables, and thus does not scale well with more variable settings. The predicted differential equations by each algorithm are in Appendix D.2.

6 Conclusion

In this research, we propose VSR-DPG to accelerate the discovery of governing equations involving many variables, which is beyond the capabilities of current state-of-the-art approaches. VSR-DPG follows a vertical discovery path, building equations involving more and more input variables using control variable experiments. VSR-DPG has the potential to supercharge current popular approaches because the first few steps following the vertical discovery route are much cheaper than discovering the equation in the full hypothesis space. Experimental results show VSR-DPG can uncover complex equations with more variables than what current approaches can handle.

60

40

20

0

x3

Acknowledgments

We thank all the reviewers for their constructive comments. This research was supported by NSF grant CCF-1918327 and DOE - Fusion Energy Science grant: DE-SC0024583.