Learning Symbolic Models for Graph-structured Physical Mechanism
Authors: Hongzhi Shi
Jingtao Ding
Yufan Cao
quanming yao
Li Liu
Yong Li
Venue: ICLR 2023 poster (2023)
Keywords: Symbolic Regression
Graph Neural Networks
Physical Mechanism
Message-Passing Flow

Abstract: Graph-structured physical mechanisms are ubiquitous in real-world scenarios, thus revealing underneath formulas is of great importance for scientific discovery. However, classical symbolic regression methods fail on this task since they can only handle input-output pairs that are not graph-structured. In this paper, we propose a new approach that generalizes symbolic regression to graph-structured physical mechanisms. The essence of our method is to model the formula skeleton with a message-passing flow, which helps transform the discovery of the skeleton into the search for the message-passing flow. Such a transformation guarantees that we are able to search a message-passing flow, which is efficient and Pareto-optimal in terms of both accuracy and simplicity. Subsequently, the underneath formulas can be identified by interpreting component functions of the searched message-passing flow, reusing classical symbolic regression methods. We conduct extensive experiments on datasets from different physical domains, including mechanics, electricity, and thermology, and on real-world datasets of pedestrian dynamics without ground-truth formulas. The experimental results not only verify the rationale of our design but also demonstrate that the proposed method can automatically learn precise and interpretable formulas for graph-structured physical mechanisms. 

Published as a conference paper at ICLR 2023
LEARNING SYMBOLIC MODELS FOR GRAPH- STRUCTURED PHYSICAL MECHANISM
Hongzhi Shi', Jingtao Ding', Yufan Cao', Quanming Yao', Li Liu", Yong Li'
Department of Electronic Engineering, Tsinghua University! Biren Tech²
Liyong07@tsinghua.
ABSTRACT
Graph-structured physical mechanisms are ubiquitous in real-world scenarios, thus revealing underneath formulas is of great importance for scientific discovery. However, classical symbolic regression methods fail on this task since they can only handle input-output pairs that are not graph-structured. In this paper, we propose a new approach that generalizes symbolic regression to graph- strucurtr physical mechanisms. The essence of our method is to model the formula skeleton with message-passing flow, which helps transform the discovery of the skeleton into the search for the message- that able to search message-passiing passing flow, which is efficient and Pareto- optimal in terms of both accuracy and simplicity. Subsequently, the underneath formulas can be identified by interpreting component functions of the searched peassage-asssig flow, reusing classical symbolic regression methods. We conduct extensive experiments datasets from different physical domains, including mechanics, electricity, and thermology, and on real-world datasets of pedestrian dynamics without ground-truth formulas. The experimental results not only verify the rationale of our design but also demonstrate that the proposed method can automatically learn precise and interpretable formulas for graph-structured physical mechanisms.
flow. Such a transformation guarantees
INTRODUCTION
For centuries, the development of the natural sciences has been based on human intuition to abstract physical mechanisms represented by symbolic models. i.e., mathematical formulas. from experi- mental data recording the phenomena of nature. Among these developments, many mechanisms are naturally graph-structured (Leech) 1966), where the physical quantities are associated with individual objects (c.g., mass), pair-win pair-wise relationships (e.g.. force) and the whole system (e.g., overall energy), corresponding to three types of variables on graphs: node/edge/global variables. For example, as shown in Figure T [a), the mechanical interaction mechanism in multi-body problem corresponds to graph with masses (m:),, positions (V:) as attributes of nodes, and spring constants (kiz) attributes of edges, which, together with the graph connectivity, yields the acceleration as output attributes of nodes; while in the case of resistor circuit, nodes and edges correspond to voltages and resistances, respectively, and these attributes define a graph-level overall power of the circuit.
In the past few years, Symbolic Regression (SR) (Shooo et al. 2018 Schmidt & Lipson 2009
Udrescu et al. 2020), which searches symbolic models = F(z) from experimentally obtained onput-output pairs [ y)} with JF being an explicit formula, has become a promising approach trying to automate scientific discovery. Traditional SR methods include genetic programming-basse methods (Schmidt Lipson| 2009) |ortin al. 2012) working by generating candidate formulas by "evolution" (i.c., manipulations), and deep lcarning-bass methods (Li et al. [019 Biggio et al.) 2021; Zheng et al.] (202I) utilizing sequence models to generate candidate formulas. However, these methods are designed for traditional SR problems input-output pairs {(x, y)} without considering graph information.
To exploit the inherent graph structure in physical mechanisms, shown in Figure [b), SR on graphs aims to find a formula F that characterizes a mapping from input {g. X} to output y. with X and both inside graph structure G. To perform this, need both fine exploitation of inherent

---

Published as a conference paper at ICLR 2023
connecsed multisbodl problam
Grnaph
Case Resistor circait
S.G. where (td, P,E.EEE,R
X-inputs youmipul
Grgraph topokogy Plende varsiblle E, E.E.salg variables, global variables
(a) Examples
(b) Traditional SR v.s. Graph SR
(c) Diverse Formula Skeletons
Figure 1: (a) Two examples of graph-structured physical mechanisms. (b) Illustration of traditional SR (left) and SR on graphs (right). (c) Two different formula skeletons in the two examples, where all formula components are interconnected accordingly construct an overall formula.
graph structures of physical mechanisms and well achievement of flexibility regarding diverse forms of interaction between entities in physical world. Graph Neural Network (GNN) has recently been incorporated into SR for discovering mechanisms behind particle interactions (Cranmer ct al. 2020; Lemos et al. 2022). However, obvious setbacks exist that the message-passing passing flow the formula skeleton, required to be manually designed to learn the underlying mechanisms, is impractical because the formula skeletons usually remain unknown and significantly different in diverse physical domains shown in Figure (e).
of GNN, corresponding
To solve this problem, inspired by the correspondence between the skeleton and pessage-passing flow in GNN, our idea is transform the discovery of the skeleton into the search for message-passing flow, which paves the way for identifying the underneath formula by interpreting each component function in the searched message-passinng flow. However, due the coupling relationship between the skeleton and the component formula in the skeleton, neither of them can be independently identified, implying vast, highly entangled search space for both gesssaeepassingg flow and component functions. To tackle this challenge, we formulate a bi-level optimization problem that searches for the message-passinng
flow by pruning strategy at the upper level on condition that its component functions have been optimized with deep learning (DL) at the lower level. Besides empirical accuracy, is equally vital but non-trivial to maintain explicit interpretability and generalization ability in discovered formulas. We propose to search the Pareto-optimas accuracy and simplicity by carefully designing a scoring function involving a complexity function of message-passing flows that optimizes both aspects across different searching steps.
cessage-passiing flow between
Our contributions can be summarized as the following three aspects,
We generalize the problem of learning formulas with given skeletons (inductive bias) from graph data Cranmer et al. (2020) by additionally learning the formula skeleton from data, which essential for learning graph-structured physical mechanisms from diverse physical domains.
We propose a novel method to learn graph-structured physical mechanisms from data without knowing the formula skeleton by searching the Pareto-optimal bessaae--passing flows of GNN together with the symbolic models as components.
We conduct experiments on five datasets from diverse physical domains, including mechanics, electricity, thermology, and two real-world datasets about pedestrian dynamics, demonstrating that our model can first automatically identify the correct skeleton based on collected data instead of expert knowledge and then learn the overall symbolic model for corresponding graph-structured physical mechanism.
2
THE PROPOSED METHOD
Before introducing the proposed method, we first formally define the the problem of symbolic regression on graphs.
2

---

Published as a conference paper at ICLR 2023
Pewenatt
A
my R
111
Interaction
Rgy
Vi, 3i
*
Destination
Pedesttions
Obotacies
P
g²: Tn (V.
gift By = (V: - 113
?
p Fi = 7,f1 Tieat
porm: Peneman Enj Pa
4²: 0i Fi/oo
identity
Figure 2: Examples of three different graph-structured physical mechanisms and their corresponding message-passing flows. Notice here that the three examples differ from each other greatly in the layouts of messagg-passing flows. Every component of the message-passing flows represents certain function (listed at the bottom) on the graph abstracted from the physical scenarios given in the first row. (I) Left: the multi-bodt body kinematics scenario where the skeleton is the same in Cranmer et al. (2020) (2) Middle: the resistor circuit scenario where the skeleton is largely different from that in Cranmer al. (2020); (3) Right: the collision avoidance scenario where the skeleton is unknown.
Definition (Variables on Graphs) The topology of the graph denoted as G, and its variables include {V, E, u), where V denotes the set of node-level variables, E denotes the set of edge-level variables and u E RT denotes a global variable which interacts with elements in and through topology G. Specifically,
E Rh is the variable associated with the i-th node while eij E Rh is the variable associated with the edge connecting the i-th and j-th nodes. nhy Das and ne denote the dimensions of the global, node, and edge variables.
Definition (Symbolic Regression on Graphs) Given set of (1, where X, C (V. E, u}, which are known variables, y, E (V', E', u'}, which are unknown variables, and variables with prime denote output variables, we aim to find an accurate and compuct formula F()) that fit = FGG X).
2.1 MODEL FORMULA SKELETON WITH MESSAGG--ASSIINE FL.W
As shown in Figure [c), formula F describing graph-structured physical mechanisms can always be decomposed into several formula components, each representing its association with an individual node, pair-wise relationships between nodes or the whole system. As illustrated in the specific diagram, these components are interconnected according to the variable dependency, termed as "skeleton". Moreover, differences between the two examples also indicate the existence of diverse skeletons underlying different physical mechanisms.
The key insight is that the skeleton has strong correspondence with the message- saaaugge-pssiinggg.. flows in GNN, which differs a lot in various physical scenarios, including mechanics (Sanchez-Gonzalle al. 2020), electricity (Zhang et al., 2019), thermology (Chamberlain al.) (021). The message- passing flows can be diverse by cascading multiple blocks, changing/removing some functions and adjusting the embedding sizes. A block of full GNN contains the updating of edges, nodes, and graph
representation respectively as follows (Battaglia et al. 2018
e = 6² (eij, Vi. u)
e p (E
V Φ" (e), v)
e' p tu F
(1)
u S" (
V p (V')
where <b-) denotes the message functions and p(-) denotes the aggregation functions. We provide three examples in different physical scenarios, as shown in Figure 2 to illustrate the well- well-defined analogy between message-passing flows and skeletons.
3

---

Published as a conference paper at ICLR 2023
Table 1: Analogy between graph-structured formulas and GNN.
computation steps of graph structured formulas
|architecture of prsssage-passiiee networks
# intermediate variables
embedding size
dependent variables
message-passing connection
formulas with dineliidmmensionaa) values as input formulas with set as input
message function ø
aggregation function p
Example 1 (Mechanies: Multi-hody Kinematics) In this problem, aim find the particles' acceleration. We have an edge update function 6f since particle pairs determine string forces, while the aggregation function p of edges is based on the independent action principle of force.
Example (Electricity: Resistor Circuit) The objective of this problem is to find the overall power of a given resistor circuit. An edge update function of corresponds to the computation of single- resistor power utilizing Joule's s Law, and an aggregation pa from edge to global appears for summation get overall power.
Collision Avoidance) In this problem, we aim to find the pedes- trians' acceleration according to their postions and velocities. The formulas including the skeletons to describe this relationship can be diverse, which highly depends on the pedestrian scenarios.
Example 3 (Pedestrian Dynamies:
2.2 TRANSFORMING INTO THE TASK OF
FLOW SEARCHING
The cassage-passing flows of GNN correspond explicit meanings the symbolic calculation for graph-structured mechanisms, which is summarized in Table This strong resemblance inspired us to devise a transformation of the primitive SR task on graphs into a relatively more practical task of searching pressag--passing flows. Our model has two stages: message-passee and message-passinggg flow-based SR. Specifically, at stage 1, we need to search the message message-passing flow the formula skeleton. Then stage 2, we need to symbolize components into formulas and cascade them according to the skeleton to get the final graph-structured mechanisms.
flow searching
With the above transformation, it is clear that need to solve the following bi-level optimization problem in stage 1, i.e.,
[P", M arg min sP;
(2)
M.P
where M denotes the message-passing flows, P denotes the parameters in the DL components, M
and P denote the Pareto-optimal one, and s(.) gauges how well the finally learned formula obtained by M and P performs. However, there are two challenges to learning formulas graphs: (i)
considering simplicity and accuracy simultaneously for graph SR is difficult; (ii) the discrete search space composed of skeletons and component formulas is prohibitively huge.
2.3 SEARCHING MESSAGE- PASSSIG FLOWS (FORMULA SKELETONS)
To deal with the first challenge, motivated to change equation into
M arg min s(P": M)) (M";
(3)
M
s.t. P arg min 1 (P; M),
(4)
P
where equation 3 and equation 4 are two optimization problems in the upper level and lower level,
s(P*; M)) I(P M) + >c(MM is the score taking both simplicity and accuracy into consideration,
2(-) denotes the error loss of predicting outputs (see details in Appendix B.4), X is the weight, and c(-)
denotes the complexity of pesssage-passing flow. The design of complexity c(-) is flexible (see details in Appendix B.5). and we calculate the complexity as follows, (i) for each layer (corresponding to
function in (uu Φv: ())), the complexity can be calculated as the product of the embedding size of this layer and the number of inputs in this layer; (ii) the whole complexity of the message-pasee flow can be calculated as the summation of the complexity of each layer. The optimization at the lower level is solved by training parameters P of DI. model given the structure M, while the
4

---

Published as a conference paper at ICLR 2023
F(V.8)
Full CNN Blockis) Ciscadal
Parcio-opiiaal Mevagge Paecing Flow
SR 4
pame
0
csscale
= ga
(a) Message-Passing Flow Scarching
(b) Message-Passing Flow-based SR
Figure 3: The framework of our model. We partition the task into two stages: (a) Message- Passing Flow Searching and (b) Message-Passing Flow-based SR.
optimization at the upper level M is difficult because M forms huge discrete search space including # blocks, # message-passing layers, connections and embedding sizes.
Another insight facilitating the dealing with the second challenge that if the message-pasien passing flow is a super-structure of the ground-truth one, i.e., redundant computations are done, it results in merely subtle variations of the loss. Still, the pressage-passingg flow is sub-structure of the ground truth, i.e. some necessary computations are missing, the loss jumps up with a change of magnitude.
Such an observation together with equation B facilitates us first search an initial wressage-passing flow that is the super-structure of the ground-truth and then learn to prune the message- paassggggssinng flows to get both compact and expressive cessaae--assingg flows. The framework of our model is shown in Figure 3)a
a), which sequentially search # blocks, # layers, connections and embedding sizes in hierarchical way and the four steps in detail are listed as follows.
Blocks. First, we need to find message-passing flow that the ground-truth message-passing flow its sub-structure. To achieve this goal, we first stack several full message- message-passing blocks. By optimizing the score in equation 3 find the Pareto-optimal number of passsage-passigg blocks, where we take the number of blocks as the complexity in equation 3 and RMSE between the predicted value and the ground-truth the loss.
Step 1: Search Message-Passee
Step 2: Search Message-Passing Layers. As mentioned, a full message- message-passing block contains three layers corresponding updates of edge, node, and graph representation. However, not all of them are necessary for obtaining the output. To find the most compact layers, we try to delete each layer to see whether the score in equation jincreases or decreases, where we define the complexity as the number of layers and RMSE as the loss. Our pruning-based searching method is based on the unique insight in SR, which is much more efficient than brute force search. Specifically, our method can significantly decrease the computational cost from O(2") ) to O(n), where n is the number of initial layers.
Step 3: Search Necessary Inputs. We further filter out the useless inputs for each layer. Specifically, we adopt similar strategy to the previous searching step: try to delete each input see whether the score in equation 3 will rise or drop, where the complexity is the number of connections and the loss RMSE. Similar to step 2, our model can significantly decrease the computational cost from ((2)) to Onn), where n is the number of initial inputs.
Step 4: Search Embedding Sizes. To ensure that the embedding in each layer is compact and with explicit physical meanings, we use the score given in equation to find the Pareto-optimal embedding size for each embedding, where the complexity is defined as the embedding size and RMSE defines the loss. We try to reduce the embedding size to find the embedding size with the highest score. At the same time, we fix other embedding sizes as large enough number to ensure the information bottleneck can only be caused by the embedding size are searching for.
2.4 THE LEARNING PROCEDURE
After obtaining the message- message-passing passing flow M and the parameters of DL component function P* at the first stage, we follow Cranmer et al. (2020) to symbolize each DL component into formulas, and then cascade them according to the skeleton represented in M into a whole formula at the second stage, as shown in Figure b). For aggregation functions corresponding to set functions, (i) we choose several commonly used aggregators candidates, including sum, mean, and max, while other aggregators can be generated by them, and select the maximum operator to replace the softmax function, (ii) perform SR onput-output pairs from trained GNN component functions, (iii)
5

---

Published as a conference paper at ICLR 2023
fine-tune all constants in the overall function (given by cascading component functions), thereby avoiding the accumulation of errors.
EVALUATION
3.1 EXPERIMENT ON CLASSICAL PHYSICAL SCENARIOS
Dataset. We utilize five datasets in different physical domains demonstrate that model has the ability to rediscover the well-known graph-structure physical mechanisms, as introduced in Section 2.1 and Appendix A. We provide two cases of mechanics scenarios, one of electricity scenario and two of thermology scenarios. For both mechanics and thermology scenarios, there are two selected cases with different complexity, where the content listed in parentheses associated with the more complex scenarios. Detailed information about formulas and data generation is reported in Appendix A.1
Metries. Given the same input, we use the coefficient of determination R², indicating the proportion of the output variation that can be predicted from the input variables. Specifically, it calculated by the output of distilled formula and the output of the ground-truth formula to measure whether the learned formula is accurate enough. Rt can be calculated R 1 Σ, (yi 0.)² Σ, (3/ 7) where ÿ = E, y:/
Comparing Methods. We compare model with learning symbolic models from deep learning with inductive bias (SymDL) (Cranmer et al.] 2020) to demonstrate that our model is flexible in scenarios and variant of our model that uses full graph network without pruning searching (FullGN) for the ablation study (1-layer full GNN that removes non-exist inputs and non-required outputs). The message-passin passing flows of SymDL and FulIGN shown in Appendix A.3
Mechanies (simple)
Mechanies (complex)
Method
Learned formulas
R2
Learned formulas
R2
SymDL (Cranmer, 2020)
<X + suin (roiig iifhett srtan cafooi
wM E (Cos mylldry Gubaulr munu/Cidry doul) OoPowalm
L.000
0.203
(dii mullany Gdaagir- muruCiinn CiPunaat,
*X - su (ooee dblest (hhak Kadyvaa huiii (i0iaa dutes (0 GuAovlt)
FullGN
1.000
0.529
a:ll
Giihee (aaal//our) JAn, G.Auu) (1 CiP
*X s, pooaa e diaest jhals CadoaWo
*it huin 2anaa Gine (o So 2 =.l Faaaas
Ours
1,000
0.917
4.00 Ei0i dutesl GAAvol))
Electricity
Thermology (simple)
Thermology (complex)
Learned formulas
R2
Learned formulas
R2
Learned formulas
R2
n A. S N0m-vooe--e
FulIGN
0,587
0.318
as
Ca Ni Cul Mitmnl Di
0.189
of
Ours
1.000
38
>
Owel - -20
1.000
As
o Eoni To csa..)
0.993
Figure 4: Correctness comparison based on learned formulas and R2 in five scenarios. The formulas learned by our model are the same as the ground-truth formulas
Mechanics(simpp)
Pechaniccsccompp)
Electricity
Thermology(sipp)
Thermalogy(comp.
Figure 5: Learned sassage-passing flows by our model.
Plausibility Comparison with Baselines. We first compare the applicability of our method with the SOTA baseline (Crnmmer et al. 2020 As listed in Table 4 method can be applied to all

---

Published as a conference paper at ICLR 2023
Table 2: The running times (mins) for two stages in different scenarios.
Scenario
Mechanics (simple) (complex)
Electricity
Thermology (simple) (complex)
SymDL
69.91
70.82
FullGN
86.25
113.59
59.23
67.21
69.23
Stage
11.59
18.42
11.26
11.37
12.68
Ours
Stage 2
68.75
93.12
28.96
30.88
45.27
Total
80.34
111.54
40.22
42.25
57.95
five cases from three physical scenarios, including mechanics, electricity, and thermology, while the baseline fails in the last two scenarios due to incorrect message-passing flow because their designed posssage-passing flow is designed explicitly for Newton-force interaction in the simple mechanical scenario and not flexible enough for other scenarios.
We design two cases in mechanics scenarios: calculating the acceleration and the relative acceleration in the center-of-ng mass frame. SymDL is designed for handling formula discovery in the simple case, flow. Comparatively, our method moves forward and challenging setting without specifying a message-passing flow representing the formula skeleton. To ascertain the correctness of the learned formulas, besides the baseline, we further introduced variant of our method without searching the pressage-passing flow, i.e., directly using the full message- passing block in the first stage. As listed in Figure 4 our model achieved the same performance with two baselines in the simple mechanics case. In the complex case, our model outperformed the SOTA baseline and the variant of our model by large margin w.r.t. R² metric. Specifically, these two competitors both failed and got a rather low R², while our method rediscovered the correct sessagnn--assingg
with specified
more general
formula with R2 0.917, indicating the advantage of searching for the correct flow. For our model, the difference between two formulas for two cases is about the latter two terms corresponding to the additional message-passsingg flows V - V and V - u' - V', which SymDL cannot handle. The formulas learned by baselines wrong for lack of necessary dependencies, which fail to have physical meaning and differ largely from the ground truth.
Our problem differs from SymDL, requiring prior knowledge of the formula skeleton for designing the deep learning architecture, which is almost impossible to know in new real-world scenarios. For the rest three cases that the SOTA baseline cannot handle, as listed in Figure 4 the performance gain over the variant using full message-passing flow indicates that optimizing Pareto-optimal score is essential in obtaining correct formulas, which is less subject redundant message- pasiggeg-aainggggg flows that hinder the subsequent SR process, including unnecessary inputs and redundant computation steps. The detailed searching process of the electricity case is analyzed in Section [.1] For the complex case of thermology, it can be observed that the learned formula successfully captures the effect of externally conducted heat compared with the simple case, while other baselines fail to have physical meanings due to unnecessary inputs and redundant computation steps. Although the change is slight (whether external heat conduction exists), the skeleton and the entire formula are quite different, and so is the entire formula.
To be more precise, the learned message-passing flows by our model are shown Figure 5 Besides, the time cost of each independent part is shown in Table 2 where can observe that searching message- passagg-passing flows only takes a small part of the whole procedure, and our model's runming times similar to SymDL and usually shorter than FulIGN. Furthermore, we conduct experiments to demonstrate the design philosophy of our method, which is reported in Appendix [A.9
Qualitative Results for Understanding the Searching Process. We show the searching process of message-passing flow in Figure 6) where the upper row shows the learning curve in terms of error (RMSE), complexity, and score, which is weighted summation of error and complexity. From Figure 6 we observe that if the message-passing flow is a sub-structure of the ground-truth message-passing flow, the performance will drop significantly. On the other hand, message- -aassing flows with redundant layers/inuus/embedding sizes will have similar performance, echoing the rationale of our pruning strategy. The core idea is to search for the most compact message-par -assing that is expressive enough, and the four-step searching process is as follows, (i) the model tried the number of blocks as 1~3 3, finding similar errors with a rise in complexity, so it opted for 1; (ii) from the searched message-passing flow at the previous stage, it tried to delete every layer associated with 6. It turned out that deleting the edge layer would cause a huge error increase, so the edge layer was preserved, after which it tried to delete the node layer and found that the score decreased (error does not change a lot and the complexity decrease), so it decided to delete node layer; (iii like the previous stage, it tried to delete each input and found that only deleting the V - 6¹ connection
7

---

Published as a conference paper at ICLR 2023
Search #Blocks
Layers
Useful Input
Compress Representation
Enor/ /
Complexity
100
20
Scarching Time (mins)
Searched
Mlesage-Passinng
Flow
E
#Block-1
Removv Node Update Layer
Remove V-p" Connection
Reduce Emhedding Size
Figure 6: The searching process of message- saaagge-assinggggg flows in four steps in the circuit scenario.
would not cause an error increase, so this connection was deleted. (iv) finally, it tried to compress each representation and found that when the embedding size was 1, the score minimized, so 1 was chosen as the embedding size. After the whole process, the message-passing flows, including embedding (intermediate variable), functions, and topology, have explicit physical meanings, paving the way for symbolic regression.
3.2 EXPERIMENTS ON REAL-WORLD SCENARIOS OF PEDESTRIAN DYNAMICS
To better show how our model discovers unknown graph-structured physical mechanisms in the real world, we conduct carefully-designed experiments of our model formula discovery for pedestrian dynamics.
Problem Formulation. We aim find a formula that can approximately describe the relationship between acceleration a and the velocity v, pedestrian position æ, and destination position &dest- The graph G describe the interaction relationship, which constructed as follows when the two pedestrians' distance is less than R, they are connected, and otherwise, they are unconnected. Formally, the problem can be described as finding a formula F. that fits F(C, x, U, W.test).
Datasets. We conduct experiments on two real-world datasets of crowd trajectories: several experi- mental datasets from studies about pedestrian dynamics (Boltes & Seyfried 2013 I
including the following scenarios, (i) Unidirectional flow in corridor: group of people goes through corridor the same direction, as shown in Figure Saa); (ii) Bidirectional flow in corridor: a group of people goes through corridor in opposite directions, as shown in Figure b).
Comparing Models. For pedestrian dynamics, a well-known manually-desio designed model is the social force model (Helbing & Molnar 1995), ). in which pedestrians are with two forces: an attractive force from the destination and a repulsive force from the surrounding pedestrians and obstacles (refer to Appendix A.4 for details).
Learned Formulas. The learned formulas and the corresponding physical meanings are reported in Figure which demonstrate that model can learn different skeletons and formulas that are more precise than the social force model with explicit physical meanings. The performance comparison also reported in Figure 7, where we can observe that our model has about 10% improvement compared to the social force model.
RELATED WORKS
Symbolic Regression (SR). Distilling elegant symbolic expressions from vast experimental data has always been the mainstream method used for finding new formulas and verifying hypotheses throughout the history of physics. SR is classic topic (Schmidt Lipson! 2009 Petersen et al.. 20200 Biggio et al. 2021; Guimerà et al. 2020) that tries to emulate the process to learn an explicit symbolic model that can describe a projection from input to the output y/ as accurately as possible while maintaining its compactness. Traditional methods of discovering formulas from data are primarily based on genetic programming(GP) (Schmidt Lipson 2009
Koza 1994 Worm
Chiu 2013).
https://ped. fz- jz-juel: ch. chh.de/dattaaase/8kkk
8

---

Published as a conference paper at ICLR 2023
Secnario
Unidirrctional flow corridor
Hidirectiona flow in corridor
Searched
Messag-pasing
Flow
Searched
o Ri.
(((x,)))) ("((z)
Skeleton
JeNDo
JeHio
Learned Formula
4t ke Ar sign(v.) (ncemal componcnt)
wht o-Mrsign(y. (nornal component)
(Difference with
f (langential component)
6e: 6 (langential component)
Social Force)
6²: Oheenn mean(v))
*" Eyeng V (xx,,x) + (v APp)/t
q² Eeenn 9 + (v A2 A/p")/'
Physical Meanings
There repuisive force
front faster
Therr repulsive than the back Tend to dodge to the right
the front faster
the back
Tend dodge to right
Herding: Tend to follow the averuge speed of the crowd
Performance (R2)
0.71 (Social Force Model) 0.78 (Ours)
0.68 (Social Force Model) - 0.81 (Ours)
Improvements
Figure 7: Learned formulas for pedestrian dynamics in two scenarios.
Hitherto, there have been promising results yielded by GP-based SR methods such as Burlacu et al.
(2020), Virgolin et al. (2019), and the famous commercial SR method Eureqa (Dubčáková [201]] 1).
etc. More recently, methods based on DL (Zheng al. 2021 Qian et al., [021] Martius & Lampert
2016; Kusner et al. 2017| Udrescu & Tegmark 2020 Udrescu et al. 2020; Daniele et al. 2022) for
symbolic regression introduced with better expressive ability than GP. Furthermore, Cranmer
jet al. (2020) first proposed to learn graph-structured physical mechanisms (especially kinematics)
given formula skeletons. Beyond that, we propose searching for formula skeletons automatically, where existing SR methods can be exploited to look for basic components of the whole formula.
Graph Neural Network (GNN). GNN (Kipf & Welling, 2017 Veličković al. 2018 Gilmer
ct al. 2017) can be viewed in message-passing manner (Battaglia et al. 201 1 Veličkovid [022]
Bronstein al. 2017), and most of them can be summarized as message-passingg passsng among three levels: edge/node/graph level, while the message- pesssagepasssing flows and messag//agreegatiinn functions can be customized very differently based on the specific characteristics of applications. It has been widely used for physical systems by capturing the interacting mechanisms, such as simulating mechanical system (Sanchez-Gonzalez et al. 2020 Huang et al. 202 Sanchez- et al. [2018), designing circuits (Zhang et al. 2019: Ren et al.] 2020). simulating heat conduction (Chamberlain et al.
2021 Xhonmeux et al. 2020 simulating pedestrian dynamics (Shi et al. 2023 Zhang et al. al. 2022).
Furthermore, there are some works (You et al. 2020 Yoon et al.) 20200 Cai et al.] | Gu et al.)
2021) that adopt automated machine learning techniques for searching the best GNN architecture for specific prediction task. Unlike them, focus the SR problems on graphs and inspired by symbolic regression (Udrescu & Tegmark, (2020), we propose to search the Pareto-optimal message- passing flows, which is both accurate and simple and can benefit learning of symbolic models.
Pareto-optimal Search. The previous Pareto-optimal solutions proposed in Neural Architecture Search (NAS) area (Lomurno et al. 2021 Lu al.) 2020; Dong et al. 2018) focus on finding the model architecture with both high prediction accuracy and low inference latency, which does not meet requirements for solving graph SR problem. Instead, our proposed method is based on a novel insight in the SR scenario: the performance would be similar when the message- -aassgggg flow (skeleton) a super-structure of the ground-truie truth one. In contrast, the performance degrades a lot if it is sub-structure of the ground-truth one.
5 CONCLUSION
In this paper, we generalize the problem in Cranmer et al. (2020) by learning the formula skeleton rather than manually designing, which is crucial for learning formulas in a new physical area without much prior knowledge. We propose new SR method that first transforms the discovery of the formula skeleton to the search of the Pareto-optimal message-passing flow with accuracy and compactness and then symbolizes its message functions to obtain the underneath formula. We conduct experiments on five datasets from three different physical domains, including mechanics, electricity, and thermology, demonstrating that our method can consistently learn plausible formulas describing the graph-structured physical mechanism. Furthermore, to show that our model is practical for learning unknown formulas in the real world, we conduct experiments on two real-world datasets about pedestrian dynamics, which learn different formulas with explicit physical meanings for different scenarios more precisely than mainstream empirical formulas.

---

Published as a conference paper at ICLR 2023
ACKNOWLEDGEMENT
This work was supported in part by the National Key Research and Development Program of China under 2020YFA071 14033, the National Nature Science Foundation of China under 61971267, U1936217, 61972223, 62171260. Q. Yao was in part supported by NSFC (No. 92270106) and CCF-Baidu Open Fund.
REFERENCES
Peter W Battaglia, Jessica B Hamrick, Victor Bapst. Alvaro Sanchez-Gonzalez Vinicius Zambaldi. Mateusz Malinowski. Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, al. Relational inductive biases, deep learning, and graph networks. arXiv preprint arXiv: 1806::100.012000 082610161, 2018.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. arXiv preprint arXiv:2166.06-22 2021.
Maik Boltes and Armin Seyfried. Collecting pedestrian trajectories. Neurocomputing, 100: 100:127-133, 2013.
Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst. Geometric deep learning: going beyond euclidean data. IEEE Signal Processing Magazine, 34(4):18 42, 2017.
Bogdan Burlacu, Gabriel Kronberger, and Michael Kommenda. Operon C++: An Efficient Genetic Programming Framework for Symbolic Regression, pp. 1562-15700 Association for Computing Machinery, New York, NY, USA, 2020. ISBN 978 3781450371278. URL https::/doii org/10
1145/33 337999 3398099
Shaofei Cai, Liang Li, Jinean Deng, Beichen Zhang, Zheng-Jun Zha, Su, and Qingming Huang. Re- thinking graph neural architecture search from message-passing, In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 6657-6666, 2021.
Ben Chamberlain, James Rowbottom, Maria I Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi. Grand: Graph neural diffusion. In International Conference on Machine Learning, pp. 1407-1418. PMLR, 2021.
Miles Cranmer. Pysr: Fast & parallelized symbolic regression in python/julia, September 2020. URL http:/dii.orrg//// wtt.:////.....111..//
httt:://di..oo org.
4041459
Miles Cranmer, Alvaro Sanchez Gonzalez, Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. Advances in Neural Information Processing Systems, 33, 2020.
Alessandro Daniele, Tommaso Campari, Sagar Malhotra, and Luciano Serafini. Deep symbolic learning: Discovering symbols and rules from perceptions, 2022. URL
abs 0tt5/208..115
Jin-Dong Dong, An-Chieh Cheng, Da-Cheng Juan, Wei Wei, and Min Sun. Dpp-net: Device- progressive search for pareto-optimal neural architectures. In Proceedings of the European Conference on Computer Vision (ECCV), pp. 517-531, 2018.
Renáta Dubčáková. Eureqa: software review, 2011.
Félix-Antoine Fortin, François-Michte De Rainville, Marc-André Gardner Gardner, Mare Parizeau, and Christian Gagné. Deap: Evolutionary algorithms made easy. The Journal of Machine Learning Research, 1311:2277-22755 2012.
Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. Neural message passing for quantum chemistry. In International conference on machine learning, pp. 1263-1272. PMLR, 2017.
10

---

Published as a conference paper at ICLR 2023
Yu-Chao Gu, Li-Juan Wang, Yun Liu, Yi Yang. Yu-Huan Wu, Shao- Lu, and Ming- Ming-Ming Cheng. Dots: Decoupling operation and topology in differentiable architecture search. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pp. 12311-12320, 2021.
Roger Guimerå, Ignasi Reichardt, Antoni Aguular-Mogar Francesco A Massucci, Manuel Miranda, Jordi Pallarès, and Marta Sales-Pardo. A bayesian machine scientist to aid in the solution of challenging scientific problems. Science advances, 6(5):eaav6971, 2020.
Dirk Helbing and Peter Molnar. Social force model for pedestrian dynamics. Physical review E, 51 (5):4282, 1995.
Wenbing Huang, Jiaqi Han, Yu Rong, Tingyang Xu, Fuchun Sun, and Junzhou Huang. Equiv- ariant graph mechanics networks with constraints. In International Conference on Learning Representations, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
Thomas Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In International Conference on Learning Representations, 2017.
John Koza. Genetic programming as means for programming computers by natural selection. Statistics and computing, 4(2):87-112, 1994.
Matt J. Kusner, Brooks Paige, and José Miguel Hernández- Lobato. Grammar variational autoencoder. In Doina Precup and Yee Whye Teh (eds.), Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 1955-1954.
PMLR, 06-11 Aug 2017. URL
ml
html
JW Leech. The feynman lectures physics. Physics Bulletin, 17(10):367, 1966.
Pablo Lemos, Niall Jeffrey, Miles Cranmer, Shirley Ho. and Peter Battaglia. Rediscovering orbital mechanics with machine learning. arXiv preprint arXiv:2202.02300 2022.
Li Li, Minjie Fan, Rishabh Singh, and Patrick Riley. Neural-guided symbolic regression with asymptotic constraints. arXiv preprint 1eXiv:/901.07700 077714 2019.
Eugenio Lomurno, Stefano Samele, Matteo Matteucci, and Danilo Ardagna. Pareto-optimal pro- gressive neural architecture search. In Proceedings of the Genetic and Evolutionary Computation Conference Companion, pp. 1726-1734, 2021.
Zhichao Lu, Ian Whalen, Yashesh Dhebar, Kalyanmoy Deb, Erik D Goodman, Wolfgang Banzhaf, and Vishnu Naresh Boddeti. Multiobjective evolutionary design of deep convolutional neural networks for image classification. IEEE Transactions on Evolutionary Computation, 25(2):277-29 2020.
Georg Martius and Christoph H. Lampert. Extrapolation and learning equations, 2016. URL
org/ abtp//lss/1110
02995
Brenden K Petersen, Mikel Landajuela Larma, Terrell N Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathemat- ical expressions from data via risk-seeking policy gradients. In International Conference Learning Representations, 2020.
Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela van der Schaar. D-code: Discovering closed-form odes from observed trajectories. In International Conference on Learning Representations, 2021.
Haoxing Ren, George Kokai, Walker J Turner, and Ting-Shens Ku. Paragraph: Layout parasitics and device parameter prediction using graph neural networks. In 2020 57th ACM/IEÉE Design Automation Conference (DAC), pp. 1-6. IEEF, 2020.
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and PMLR, 2018.
control. In International Conference on Machine Learning, pp.

---

Published as a conference paper at ICLR 2023
Alvaro
Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Riedmiller, Raia Hadsell, and Peter Battaglia. Graph networks as learnable physics engines for inference and control. In International Conference on Machine Learning, pp.
PMLR. 2018.
Alvaro Sanzhez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International Conference on Machine Learning. pp. 8459-8468. PMLR, 2020.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data, 324(5923)881-85,, 2009.
science,
Hongzhi Shi, Quanming Yao, and Yong Li. Learning to simulate crowd trajectories with graph networks. In Proceedings of the ACM Web Conference, 2023.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: A physics-inspired method for symbolic regression. Science Advances, 2666) 2020.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Neural Information Processing Systems, 33, 2020.
Petar Veličković. Message passing all the way up. arXiv preprint arXiv:220..1109) 2022.
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In International Conference on Learning Representations, 2018.
Marco Virgolin, Tanja Alderliesten, and Peter A. N. Bosman. Linear scaling with and within semantic backpropagation-baase genetic programming for symbolic regression. In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO '19, pp. 1084-1092, New York, NY, USA, 2019, Association for Computing Machinery. ISBN 0781450361 18, doi: 10.1 10.1145/321707.3331755
/ wtg.ddi..or//// 1145/ /1114/////170.... 332177 3321758
Tony Worm and Kenneth Chiu. Prioritized grammar enumeration: symbolic regression by dy- namic programming. In Proceedings of the 15th annual conference on Genetic and evolutionary computation, pp. 1021-1028, 2013.
Louis-Pascal Xhonneux, Meng Qu, and Jian Tang. Continuous graph neural networks. In International Conference on Machine Learning, pp. 10432-10441. PMLR, 2020.
Minji Yoon, Théophile Gervet, Bryan Hooi, and Christos Faloutsos. Autonomous graph mining algorithm search with best speed/accuracy trade-off. In 2020 IEEE International Conference on Data Mining (ICDM), pp. 751-760. IEEE, 2020.
Jiaxuan You, Zhitao Ying, and Jure Leskovec. Design space for graph neural networks. Advances in Neural Information Processing Systems, 33:17009-17021, 2020.
Guo Zhang, Hao He, and Dina Katabi. Circuit-gma gnn: Graph neural networks for distributed circuit design. International conference on machine learning, pp. 7364-7373. PMLR, 2019.
Guozhen Zhang, Zihan Yu, Depeng Jin, and Yong Li. Physics-infused machine learning for crowd simulation. In Proceedings of the 28th ACM SIGKDD Conference Knowledge Discovery and Data Mining, pp. 2239, 2022.
Wenqing Zheng, Tianlong Chen, Ting- Hu, and Zhangyang Wang. Symbolic learning to optimize: Towards interpretability and scalability. In International Conference on Learning Representations, 2021.
12

---

Published as a conference paper at ICLR 2023
A EXPERIMENTS
A.1 DATA GENERATION
Besides the scenarios of mechanics and electricity, we further illustrate the scenario of thermology as follows.
Example 4 (Thermology: Heat Conduction) The objective of this problem is to compute the en- tropy production rate. The edge update function 6e corresponds to Fourier's s Law of Heat Conduction. and the node update function & corresponds to the Clausius entropy expression, followed by an aggregation p nu that sums up individual entropy production rates.
In different scenarios, we devised different inputs and computed the theoretical outcome according to the known mechanisms.
In the scenario of Mechanics, we randomly (standard normal distribution) set the (x, y) coordinates of the particles in 2-D cases. We assign the masses of the particles randomly according to the lognormal distribution. The original lengths of the springs are all set to I. the complex case, external forces with all dimensions following the standard normal distribution randomly. We then compute the accelerations each particle with Hooke's Law, the independent action principle of force, and Newton's 's Second Law of Motion.
exerted. Graph topology is picked
In the scenario of Electricity, we randomly give topology on the graph and set the electric potential of nodes following the standard normal distribution. Resistances of resistors edges chosen uniformly randomly from 0.01 to 1.01 to avoid extremely large power outputs. We then compute the power of each edge (resistor) according Joule's Law and add them up reach the overall power of the resistor circuit.
In the scenario of Thermology. the graph topology is given as 'grid', echoing the core idea of Finite Element Analysis. We randomly set the temperature of each node between 0 and and the thermal conductivity between 1 to globally. We then compute the discrete laplacian on the grid and the heat flow according to Fourier's 's Law of Heat Conduction. With each node's heat flow and temperature, we compute their entropy production rate separately and add them up to reach the overall entropy production rate.
The basic information of our used datasets is listed in Table 3
Table 3: Basic Information about Datasets.
Domain
Mechanics
Electricity
Thermology
Scenario
Spring connected particles
Spring-connected
Resistor circuit
Thermal system
Topology
Unoriented graph
Unoriented graph
Grid
Node Input
Mass, position, (force)
Electric potential
Temperature, (external heat)
Edge Input
Spring constant
Resistance
Thermal conductivity
Global Input
None
None
None
Output
Acceleration (relative to mass center)
Overall power
Overall entropy production
# Graphs
100
100
50
# Nodes
5
9
X- Y pairs
100000
100000
100000
A.2 REPRESENTATIVE SNAPSHOTS OF PEDESTRIAN DATASETS
To better understand the pedestrian scenarios, we show two representative snapshots of two pedestrian datasets in Figure & unidirectional flow in corridor and bidirectional flow a corridor.
A.3 BASELINE DETAILS
The message-passing flows of baselines, SymDL, and FulIGN are shown in Figure and their applicability in different scenarios are demonstrated in Table
13

---

Published as a conference paper at ICLR 2023
(a) Corridor, Unidirectional Flow.
(b) Corridor, Bidirectional Flow.
Figure 8: Representative snapshots of two pedestrian datasets.
V
SRindBias
FullGN Mechanics
FulGGN InElecotdly
Pullok in Thermology
Figure 9: message-passing flows of SymDL and FulIGN.
A.4 DETAILS OF SOCIAL FORCE MODEL
In the social force model, the baseline model for pedestrian scenarios, the dynamics of pedestrians are driven by two factors: (a) a pedestrian attracted by his/her destination with force F FD: = (Uace: vT) e: [x,
where Vai is the value of desired velocity, is the current velocity, is the relaxation time, e; is the unit vector toward the destination; (b) the nearby ones repulse a pedestrian with force F.y = A, exp (-ray B,) ey. where F, is the repulsive force, Tıj is the distance between pedestrian i and j. The joint force is F; FD: + EjEN; Fij. where AG means the set of pedestrians whose distance to pedestrian is less than 5 meters. The social force model is widely used as the foundation of much commercial software such viswalk and and anylogid In this paper, we assume that the mass of pedestrian is 1, and thus hand, the social force model is manually designed, which may have discrepancies with real-world pedestrian dynamics. On the other hand, different scenarios usually have very different pedestrian interaction mechanisms, which one single model cannot precisely model. So it is meaningful to learn data-dríven formulas to describe different pedestrian interaction mechanisms.
[x,-xx]| x
Fi/m = F.. However, on the one
A.S IMPLEMENTATION
We implement our model Python using Pytorch library and optimize all the models by Adam optimizer (Kingma & Ba 2015). We use parallel symbolic regression in Python (PySR Cranmer 2020 to extract formulas from each message functions ø.
A.6 PARAMETER SETTINGS
For the DL part, we set the learning rate 10-4, tolerance in early stopping 10, #layers and embedding size in MLP as 4 and 16, the max number of epochs as 20000 and the weight À as 0.1. The choice of parameter À is analyzed in Appendix A.8 For the traditional SR part, our candidate symbols include both binary operator {+, x,/3 and unary operator (sign, exp), and we set batch size 1024.
14

---

Published as a conference paper at ICLR 2023
Table 4: The applicability in different scenarios.
Method
Mechanics (simple) (complex)
Electricity Thermology (simple) (complex)
SymDL
X
Ours
Qurs Random
Score
10-
1-
0
10
30
40
50
60
80
90
100
110
Searching Time (mins)
Figure 10: The searched best score v.s. searching time for two different searching methods in the circuit scenario.
A.7 COMPARING HIERARCHICAL PRUNING WITH RANDOM SEARCH
Furthermore, to demonstrate the effectiveness of our search method, we compare it with the random searching strategy and plot their searching processes in Figure 10] From that, we can observe that our method is much more efficient than the random search algorithm, which suggests that the searching problem is difficult and that our method effectively reduces the original colossal search space. Specifically, even if the random search algorithm takes ten times longer than ours, the score of the best-searched skeleton is still 4.5 times worse than ours, and the searched skeleton is wrong.
A.8 PARAMETER SENSITIVITY
One of the most critical parameters in model is the weight À that balances the complexity and errors. We normalize the input-output pairs to make the outputs have a standard deviation of 1. Specifically, for each dimension of features, we divide the features by their standard deviation, maintaining the fitting errors with similar magnitudes. Since the complexities of different skeletons are also with similar magnitudes, the best value of À is similar for different datasets. As shown in T we test different values of À on three diverse scenarios, including the circuit scenario and two mechanical scenarios. For each value of X, use ten different seeds to train the model to test whether our model can learn the correct message message-passing flows and the correct formulas. We choose the best formula among ten formulas, so all these values of À can allow us find the correct formula in different scenarios. Among them, we choose À = 0.1 for achieving the highest success rate among all three scenarios.
A.9 DEMONSTRATION OF DESIGN PHILOSOPHY
To demonstrate the design philosophy of our method, we plot several learning curves of different message-passing flows in the simple mechanical scenario.
We first plot the learning curves of three different message-passinne passing flows in Figure I1 to the one lacking the necessary message-passinng connection, the ground-truth one, and the one with redundant message- pasging connection, respectively. We find that the ground-truth flow and the one with redundant message-passingg connection have only little variations (the RMSE is around 0.1 and the MAPE is around 1.0) in performances after the loss function converges. However, the performance of message- sassage-passing flow lacking necessary connections decreases significantly (the RMSE is about 1.4, and the MAPE is 4.0).
corresponding
We further test the impact of the embedding size on the performance, where the learning curves are shown in Figure 12 From that, can see that when the embedding size is less than a certain number
15

---

Published as a conference paper at ICLR 2023
Table 5: Success Rate in the Different Scenarios.
Value of
0.025
[0.05
|0005
0.1
0.2
0.4
learning paassage-passinggg flows
3/10
9/10
10/10
10/10
8/10
Cireuit
learning formulas
3/10
9//00
10/10
10/10
8/10
learning message-passing flows
2/10
8/10
9/10
6/10
1/10
Mechanics (Simple)
learning formulas
2/10
7/10
8/10
6/10
1/10
learning caassage-passingg flows
1/10
6/10
7/10
5/10
2/10
Mechanics (Complex)
learning formulas
1/10
5/10
6/10
5/10
1/10
Message-Passinns Flow Lacking Necessary Connection
Ground-truth Doossage-Passing Flow
Dressage-Passinggg Flow with Redundant Connection
3.00
332
RMSE
01
325
400
100m
000
Epoch
Epoch
Epoch
MAPE
www.
Fpoch
Epoch
Epoch
Figure 11: Learning curves message-passing flows.
the performance decreases significantly (RMSE more than 1.4), while the performance is similar (RMSE is around 0.2) when the embedding size is not less than 2.
Last, we whether the softmax function can learn the ground-truth aggregator. As we can see in Figure [3] the learning curve with softmax and with the ground-truth aggregator is quite similar, and we verify that the learned aggregator the same as the ground truth.
B
METHOD DETAILS
B.1
USAGE OF THE PROPOSED METHOD
The result obtained by our method is fairly stable in terms of both searching formula skeleton and symbolizing learned neural networks. First, our designed searching method can guarantee that only better graph structures imessage-passing flows) get selected during the learning process. Based on an insightful observation that the loss would increase lot when the graph structure is subset instead of superset, propose to search four components of graph strueture (block, layer, connection, dimension) by starting with a full structure and then pruning to obtain a compact one. Second, the stability of SR results has also been demonstrated by applications in Cranmer et al. (2020).
When applied to a new physical domain, we can directly use the same DL-related (such the learning rate and embedding size of MLP) old domains and slightly tune the search- related
(weight X) near the previous optimal value on old domains. The stability of the proper À value is demonstrated in Table 5 Due the randomness of DL. training and the
16

---

Published as a conference paper at ICLR 2023
Embedding Size
Embedding Size
7.00
1.0
1.0
RMSE
18-
mude
RMSE
17-
W
200
211
30
Epoch
Epoch
Embedding Size
Embedding Size
RMSF 1.4
RMSE
0/5
Epoch
Epoch
Figure 12: Learning curves w.r.t. embedding sizes.
Softmxx
Ang
30
10-
10-
w
200
300
130
200
309
500
689
300
Figure 13: Learning curves with without softmax.
genetic algorithm, we train model ten times with different seeds to get ten formulas and select the Pareto-optimal formula according to the score (see Appendix(B.6 [.6) as the final formula. Overall, the above process in new domain does not require heavy human labor.
B.2 TRAINING ALGORITHM
The training algorithm of our model is summarized in Algorithm
B.3 DETAILS OF SYMBOLIZING THE AGCREGGATIOO//EESSAGR FUNCTIONS
Symbolize Aggregation Function p. We choose several commonly used aggregators as candidates. including sum, mean and max, while other candidate aggregators such min can be achieved by min(x)= -max( -x), and root mean square be represented by (mean(z²)) Harmonic mean and I-norm can also be decomposed into two ø functions and mean aggregators. These operations can be learned by Φ.
17

---

Published as a conference paper at ICLR 2023
Algorithm The process of obtaining the graph-structured symbolic model.
Require: including training data {(X, Y, G)
1: (start stage I) search message-passing flow follows Section 2.3
2: step I: search the Pareto-optimal block number;
3: step 2: search the required layers;
4: step 3: search the necessary input variables;
5: step 4: search the Pareto-optimis embedding sizes;
6: train the model with Pareto-optimal message-passing networks;
7: (start stage 2) Symbolize Aggregation Function / and Φ
8: replace each p by the aggregator with largest weight;
9: train the model and record the input-output pairs of each o:
10: replace each ø by a formula obtained by classic SR, with constants left to be fitted on the data;
fit the correct constants from the data by gradient descent and get the final graph-structured
symbolic model.
Symbolize Message Function ø. After training the GNN model, we use the tool of symbolic regression to extract symbolic model with correspondence to message and update functions. Specifically, for message functions in the neural network, we record their input-output pairs and implement classic symbolic tools symbolize them.
Retrain and Fine-tune Constants. Finally, to eliminate accumulated errors, we set all constants in the entire formulas parameters to be trained and fine-tune them get the final graph-structured mechanism represented by symbolic models. Specifically, after extracting the message-passing flows, replace the MLP in message/update functions with the corresponding formulas. And set the constants in the formulas as parameters in the deep model to optimize (all the operations in symbolic models are differentiable).
B.4 DETAILS OF ERRORS
In this paper, we use RMSE the enror measure. RMSE can be calculated as
E (yA 0z)²
RMSE
where 9J is the i-th predicted value and y/s is the i-th ground-truu -ruth value and i 1, n
B.5 DETAILS OF COMPLEXITY
The design of complexity measurement for message-passing flows is flexible. In this paper, calculate the complexity as follows, (i) for each layer, the complexity can he calculated as the product of the embedding size of this layer and the number of inputs in this layer; (ii) the whole complexity of the message-passing flow can be calculated as the summation of the complexity of each layer. To illustrate, the complexity
four message-passee passing flows in the Figure can be calculated as
2 x 3 + 2 x 3 + 2 x = 14, x 3 + 2 x 1 8, 2 x 3 + 1 x 1 = and 2 x x 1 = 3. As we can see, the complexity consistently decreases in the search process.
B.6 DETAILS OF SCORE
The score of a graph-structured formula is + Ac, where / is RMSE, is the complexity of the graph-structured formulas, which defined as the complexity of cessage-passsigg flow multiplying the average complexity of component formulas.
C DETAILED RELATED WORKS ON SR
Martius & Lampert 2016) proposed model named EQL and extracted the symbolic formulas with
neural network with symbolic models as building blocks. Kusner et al. (2017) managed to eschew the problem of discrete optimization by converting discrete data into parse tree. AI Feynman (Udrescul
18

---

Published as a conference paper at ICLR 2023
& Tegmark 2020
Udrescu al. 2020) split the function into sub-functions. and performed regression on each module separately. The partition of functions was achieved with a trained neural network. SymDL (Cranmer et al. 2020) exploited the inductive biases of the correspondence of physical mechanisms (kinematics, especially) to GNN structure and established the method PySR to tackle the problem. They first revealed the link between GNNs and physical mechanisms. Unlike these works, our model can learn graph-structure physical mechanisms without requiring the information of formula skeletons, which can hardly be obtained in new physical scenarios for discovery.
19


Reviews:
[{"review_id": "s3a07vTSkN", "reviewer": "Reviewer_oCCc", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper improves SymDL by searching the Pareto-optimal message-passing flows to learn an additional formula skeleton. Such improvement empowers the method with more generality (requiring no prior knowledge), compactness, while maintaining correctness. Additionally, this paper also extends the applicability to Electricity and Thermology. Overall, the task formation is very interesting and the idea is simple but effective.", "summary_of_the_review": "As mentioned in my reviews regarding strengths and weaknesses, I think this paper has tackled a more practical and challenging problem than existing methods. The proposed method is also with novelty and soundness. Moreover, besides finding correct formulas in several known mechanisms, the authors further use it to discover new formulas for predicting real-world pedestrian dynamics. It seems that the obtained formulas are not only more precise (compared with a manually designed one) but also with good interpretability. Hence, I am leaning on the positive side.", "strength_and_weaknesses": "Strengths:\n1. The studied problem is important and interesting. AI for science is a highly promising application field that has not been well explored. This paper has made an essential footstep toward finding scientific formulas in a more general and automatic way.\n2. Compared with existing works, this paper has tackled a more practical and challenging problem of learning the formula skeleton rather than manually designing, which is crucial for learning formulas in new physical domains with less prior knowledge.\n3. The proposed method is novel, especially the idea of transforming the discovery of the formula skeleton to the search for the Pareto-optimal message-passing flow with accuracy and compactness.\n4. The two main designs are technically sound. First, the transformation from learning skeleton to searching message-passing flow is correct as the latter corresponds to explicit meanings in the symbolic calculation for graph-structured mechanisms. Second, the design philosophy behind the proposed pruning-based search procedure is reasonable, based on an insightful observation that when the searched graph structure is a subset of ground truth instead of a superset, the optimized score/loss will increase significantly.\n5. The authors conduct extensive experiments. The proposed method can find correct formulas in five mechanisms with different difficulties, while compared methods only succeed in a simple case, i.e., mechanics as in Cranmer et al. 2020. Moreover, it can discover a new analytic formula that can predict real-world pedestrian dynamics more precisely.\n6. Overall, The paper is well-written and easy to follow.\n\nWeaknesses:\n1. For easy-to-use concerns, the paper should have a specific subsection that provides a practical guide on using the proposed method when discovering formulas in new physical domains. Currently, the related information is scattered in Sec. 3.1, Appendix A.3 and A.4, not concentrated enough.\n2. To better demonstrate the superiority of the designed pruning-based search procedure, authors can consider adding further studies in the appendix. For example, how about comparing with a random search method that does not decompose searching steps and does not leverage a pruning strategy? Also, Figure 6 can be improved by plotting four steps together, with the x-axis being absolute training time. This can make the explanation text much easier to follow.\n3. The clarity of this paper can be further improved as follows,\n    - Use a specific figure to demonstrate the design philosophy of the pruning-based search procedure.\n    - Give more details on the normalization process, as in “We normalize the input-output pairs to make the outputs have variance 1” (Appendix A.3).\n    - Make related work part shorter and save space for results analysis, especially the real-world example in A.4.\n4. A few typos. e.g., on page 8, “the time cost of each independent part is shown in Table 6” should be Table 3.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "aF46WVPUhTB", "reviewer": "Reviewer_H8ss", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper presents an approach that generalizes symbolic regression\nto graph-structured physical mechanisms. As opposed to classical\nSymbolic Regression, this work assumes that X and Y in y=F(x) can be\nboth represented as graphs. The method is based on a two-level\noptimization procedure where first the formula skeleton is modeled\nwith a message-passing flow and parameters are learned. In a second\nstep, symbolization is used to each Deep Learning component.\n", "summary_of_the_review": "This paper presents an approach that generalizes symbolic regression\nto graph-structured physical mechanisms. As opposed to classical\nSymbolic Regression, this work assumes that X and Y in y=F(x) can be\nboth represented as graphs. The method is based on a two-level\noptimization procedure where first the formula skeleton is modeled\nwith a message-passing flow and parameters are learned. In a second\nstep, symbolization is used to each Deep Learning component.\n\nQuestions, comments\n-------------------\n\ns2.3, p4: for each layer: what do you mean by layer (in a general\ngraph)? Why do you need layers?\n\nFor the simple case in mechanics scenario The detailed --> ??\n\nStep 1: predicted acceleration: shouldn't the model be generic?\n\nObtaining the final formula: you haven't mentioned genetic algorithms\nin the methodology. What are the parameters? (population size,\npopulation format, fitness function etc)\n\nIs it fair to compare your model with SymDL, since their purposes and\nknowledge representation are different?\n\nThere are some logic-based systems that can also learn formulas and\nfunctions. Why aren't they included in the related work and why aren't\nthey compared against your work? (e.g. https://arxiv.org/abs/2208.11561).\n\n\nOther comments, typos etc\n-------------------------\n\nfollowing bi-level the optimization --> following the bi-level optimization\nin this layers --> in this layer (in these layers?) (2x)\n\ncannot by handled --> cannot be handled\n\ncan describe a the projection --> ??\n \n", "strength_and_weaknesses": "+ the idea of formulating the problem as a two step optimization problem, although not new, seems good in the context of this work\n- writing is a bit difficult to follow with some typos and sentences difficult to understand.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "KyKwjyQpTV", "reviewer": "Reviewer_dEBy", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "\nThe paper studies the problem of symbolic regression given the graph structure. It generalizes a previous approach (Cranmer et al. 2020) by additionally learning the formula skeleton. The method is seperated into two steps: searching the Pareto-optimal message passing flows and component-wise symbolic regression. Compared to (Cranmer et al. 2020), the second stage is not new, and the main novelty is to learn and prune the formula skeleton rather than manually designing one. Experiments on different scenarios demonstrate that the proposed method can learn better formulations and it has better prediction performance.", "summary_of_the_review": "The paper is solving an interesting problem and it provides a decent solution. However, the significance and practicality are not convincing. The clarity also needs some improvement.\n", "strength_and_weaknesses": "Strengths:\n1. It provides a good problem setting. Symbolic regression on graph-structured data is an interesting topic. \n2. The proposed method is generally intuitive and reasonable. It is easy to understand that not all flows are useful and pruning the flows to get a compact one is a good way to obtain a concise formula.\n3. Nice visulaization and good experimental results.\n\nWeaknesses:\n1. The methodology of this paper is basically following (Cranmer et al. 2020). Its generalization seems essentially a combination of NAS and Symbolic Regression. It is useful, but slightly overstated. For most of the physical systems, we know the mechanism and the design of the \"flow\" seems definite and there is no need to search it. The authors showed one case in Appendix A.4, but it is still not so attractive, because the search space of the flows is not large. The choices of reasonable flows are actually limited. \n2. The definition of complexity and its weight may largely impact the results. The model is unstable. We may need to run multiple times to get the best formula.\n3. Presentation can be improved. For example, A.4 is important to demonstrate the advantage of the method (as described in the conclusion), I believe it should be put in main body of the paper instead of Appendix. The Figure 3 is described in a wierd place. In fact, in Sec 2.3 it only described the Fig 3(a), but not the whole framework. There are also some other unclear points as below.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster