Graph-based Symbolic Regression with Invariance and Constraint Encoding
Authors: Ziyu Xiang
Kenna Ashen
Xiaofeng Qian
Xiaoning Qian
Venue: NeurIPS 2025 poster ()
Keywords: Symbolic Regression
Expression Graph
Hybrid Neural-guided Monte Carlo Tree Search
Expression Equivalence
Permutation Invariance
Physics Constraints
Primary area: machine_learning_for_sciences
TL;DR: Graph-based symbolic regression that captures expression equivalences on graph representations and incorporates constrained search utilizing hybrid neural-guided Monte Carlo tree search.

Abstract: Symbolic regression (SR) seeks interpretable analytical expressions that uncover the governing relationships within data, providing mechanistic insight beyond 'black-box' models. However, existing SR methods often suffer from two key limitations: (1) *redundant representations* that fail to capture mathematical equivalences and higher-order operand relations, breaking permutation invariance and hindering efficient learning; and (2) *sparse rewards* caused by incomplete incorporation of constraints that can only be evaluated on full expressions, such as constant fitting or physical-law verification. To address these challenges, we propose a unified framework, **Graph-based Symbolic Regression (GSR)**,  which compresses the search space through the permutation-invariant representations, expression graphs (EGs), that intrinsically encode expression equivalences via a term-rewriting system (TRS) and a directed acyclic graph (DAG) structure. GSR mitigates reward sparsity by employing a hybrid neural-guided Monte Carlo tree search (hnMCTS) on EGs, where constraint-informed neural guidance enables the direct incorporation of expression-level constraint priors, and an adaptive $\epsilon$-UCB policy balances exploration and exploitation. Theoretical analyses establish the uniqueness of our proposed EG representation and the convergence of the hnMCTS algorithm. Experiments on synthetic and real-world scientific datasets demonstrate the efficiency and accuracy of GSR in discovering underlying expressions and adhering to physical laws, offering practical solutions for scientific discovery.

Graph-based Symbolic Regression with Invariance and Constraint Encoding
Ziyu Xiang'
Kenna Ashen Xiaofeng Qian' Xiaoning Qian 1.2
Texas A&M University, College Station, TX 77840, USA Brookhaven National Laboratory, Upton, NY 11973, USA
{zxiang, kashen13, feng, xqian)@tamu edu; xqiani @alan1@bnl..m gov
Abstract
Symbolic regression (SR) seeks interpretable analytical expressions that uncover the governing relationships within data, providing mechanistic insight beyond 'black-box' models. However, existing SR methods often suffer from two key limitations: (1) redundant representations that fail to capture mathematical equiv- alences and higher-order operand relations, breaking permutation invariance and hindering efficient learning; and (2) sparse rewards caused by incomplete incor- poration of constraints that only be evaluated full expressions, such constant fitting or physical-law verification. To address these challenges, we pro- pose unified framework, Graph-based Symbolic Regression (GSR), which compresses the search space through the permutation-invarinn representations, expression graphs (EGs), that intrinsically encode expression equivalences via term-rewriting system (TRS) and directed acyclic graph (DAG) structure. GSR mitigates reward sparsity by employing hybrid neural-guided Monte Carlo tree search (hnMCTS) on EGs, where constraint-ifformed neural guidance enables the direct incorporation of expression-level constraint priors, and an adaptive e-UCB policy balances exploration and exploitation. Theoretical analyses establish the uniqueness of our proposed EG representation and the convergence the hnMCTS algorithm. Experiments synthetic and real-world scientific datasets demonstrate the efficiency and accuracy of GSR in discovering underlying expressions and adhering to physical laws, offering practical solutions for scientific discovery.
Introduction
Symbolic regression (SR) 2 is an approach to unveil the inherent dependencies governing the sys- tem under study in symbolic form. SR explores the space of closed-form mathematical expressions to find the most ideal analytical expressions that capture the underlying variable relationships from observed data. Compared black-box models (e.g., neural networks) that lack transparency or direct mechanistic understanding, SR strikes a balance between prediction accuracy and interpretability. highlighting its potential in advancing Al for scientific discovery B
However, solving SR problems is inherently NP-hard [4), due to the combinatorial nature of possible mathematical expressions (5). To improve search efficiency, recent SR methods have leveraged modern supervised deep learning (DL. and reinforcement learning (RI.) techniques 1 6 7 g)
Despite these advancements, little attention has been given to capturing the inherent invariances of expressions and effectively incorporating constraints in SR, leading to suboptimal learning strategies.
Prevailing works in SR 9 8 10 adopt semantic expression tree (ET) representations, which fail to account for the equivalence of mathematically identical expressions, breaking permutation invariance with redundant states that hinder fast convergence. Additionally, the limited structural information in ET may not fully leverage the effectiveness of DL-based encoding, limiting its expressive power.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Furthermore, existing search-based methods [T] [2] primarily enforce constraints evaluated at the operation level by incorporating priors during the sampling stage to avoid invalid expressions. However, constraints that depend on complete expressions, such as those involving constant fitting and physical law violation checking, are typically used as post hoc reward penalties instead of direct prior constraint incorporation in search. Such current practice leads to a sparse search space, making SR less effective for complex, constrained scientific research applications.
To address these challenges, we propose Graph-based Symbolic Regression(GSRR) which encodes permutation invariance and constraint priors using a Relational Graph Convolutional Network (R- GCN) on expression graphs (EGs). EGs inherently capture equivalences of expressions via term rewriting system (TRS) and its directed acyclic graph (DAG) structure, reducing redundant representations. To fully exploit EGs, we introduce a hybrid neural- neural-guided Monte Carlo tree search (hnMCTS), which balances exploration and exploitation by combining vanilla MCTS with neural- guided MCTS through an e-UCB policy, where R-GCN leverages node and edge features of EGs to encode expression-level constraint priors, allowing them be incorporated directly during the sampling process for the optimal search efficiency. Evaluating widely recognized SR benchmark datasets and real-world applications demonstrates that GSR effectively compresses the search space, mitigates reward sparsity challenges, and reduces the reliance costly constant fitting through its structured encoding. These results highlight the efficiency, accuracy, and robustness GSR in solving complex, baysscc-informed SR problems for scientific discovery.
2 Preliminary and Related Work
2.1 Symbolic Regression (SR)
Given dataset D {a yı), (222, V2).... ..., (ZN. BN)), SR aims to find mathematical expression f to map the input variable vector to the corresponding output y, with the minimum error ((:-)) over all data points in D:
min Σ T(f().). v.);
where / belongs to the mathematical expression space constructed using a predefined grammar (e.g.. arithmetic operations, functions, constants, and variables) as dictionary set $ = {10, 01,..., ). In our method, & denotes the operations or variables x,, We also include operator const for inserting constants and trans for transplantation strategy in & by default, detailed in Appendix C.
In this work, we develop GSR based on hnMCTS by leveraging vanilla MCTS 8 for lightweight exploration and neural-gune -urieed MCTS [14] for enhanced exploitation. The proposed framework not only optimizes the expression representation, enabling more structured and efficient search.
trade-off but also seamlessly integrates with graph-based
2.2 Expression Representation
SR can be formulated as Markov Decision Process (MDP), where new operator or variable is iteratively sampled based on the current expression state. The expression tree (ET) is the most commonly used state representation [2] as shown in Fig. [ In ET, operators are modeled as nodes, and edges represent hierarchical relationships between operators and operands, with trees growing from outer functions to inner ones. This structure reflects the sequential construction of expressions through iterative sampling steps.
However, ET representation suffers from a major limitation: it fails account for expression cquivalences and permutation invariance in operator generation sequences. As shown in Fig. [ two mathematically equivalent expressions with the same set of nodes can have different generation sequences, resulting in distinct parent and sibling relationships and tree structures. This redundancy inflates the state space, reducing learning efficiency, particularly for RL-based approaches (e.g.. PG and MCTS) and sequential encoding methods (e.g., RNNs), which struggle to converge due to unnecessary state diversity. Furthermore, ET encodes only basic operator types and hierarchical relationships, lacking higher-order features such as arity and operation order, making it hard to distinguish argument positions(e.g. base or exponent for power). These limitations weaken the predictive power of ba-based encoders.
2

---

0,7/0.7
4/
3,4,5
6/2
A din44) coc(R)
sin(A) (8 cns(
.101
cos
Operator
Operttor
E(sin(A) -B))
Sialeg
50h
Siilag
Figure I: An illustrative example of the permutation invariance problem in SR. (I) fi and f2: two equivalent expressions represented by ET's generated in different sequences as in the tables, yielding two distinct representations due to the different parent and sibling relationships highlighted by the colored boxes. This breaks the permutation invariance for sequential encoding methods. (II) fa: another equivalent expression represented both /1 sequence (blue numbers) and /2 sequence (red numbers), yielding the same representation, preserving the permutation invariance for the Graph Neural Network (GNN) encoder.
our proposed expression graph (EG), generated by
[6] rely either semantic ET-based representations, with little focus on improving expression representation. The authors of Al-Feynman [1 addressed this by pre- pre-training neural networks capture modularities and symmetries, decomposing SR problems into smaller sub- problems. However, this approach suffers from efficiency bottlenecks due to the need for extensive pre-training and brute-force search, which becomes computationally infeasible for complex systems. Though there are some works related to graph-based representation of expressions, such as which establishes the NP-hard nature of symbolic regression by linking symbol graph to degree- constrained Steiner arborescence problem, or Expression Directed Acyclic Graphs (DAGs) [T8] which is proposed to capture distributive equivalence and improve SR scalability, they fail to model the general equivalences, capture high-order dependencies, and preserve permutation invariance, as key insights our paper. Consequently, an effective expression representation that reduces the state's redundancy and facilitates richer structural encoding remains open problem in SR.
Most existing SR works
2.3 Constraint Incorporation
Constraints are ubiquitous in SR problems, which confine exploration to valid regions and guide SR for meaningful and robust solutions. Failure incorporate constraints properly result in sparse rewards and invalid solutions, particularly in real-world applications governed by fundamental physical laws. We classify constraints into two categories:
Operation-level Constraints: Constraints that can be evaluated based on a few operations during sampling. such as basic mathematical rules, void operations, maximum complexity (i.e., the number of operations and variables), and hand-crafted priors for a priori known physics constraints. These constraints are typically incorporated priors in the sampling stage by zeroing out probabilities of actions that would generate invalid expressions, preventing infeasible candidates from being explored.
Expression-level Constraints: Constraints that require a complete expression for evaluation, such as hidden mathematical constraints (e.g., 7 becomes invalid if "?" later sampled f(B) C' and S(B) C = 0), constants fitting constraints, and physics-based penalty functions (e.g., enforcing - 0) ao as discussed in Section 5.2). These constraints are typically incorporated as post hoc reward penalties, such as assigning zero rewards to invalid expressions after they are fully generated.
Existing SR approaches [TT 22 19] handle constraints by integrating operation-level constraints as hand-crafted priors during the operation sampling stage and expression-level constraints via reward penalty functions. However, the absence of expression-level constraint priors during sampling, especially in random roll-outs in MCTS and random crossover in GP, leads to inefficient exploration and sparse rewards, making it difficult for SR models to converge. While unit consistency can be modeled as operation-level constraint [19], the constant fitting process required for unit consistency remains an expression-level constraint, necessitating computationally expensive evaluations for every sample. This results in low efficiency and limits the scalability of SR methods. Hence, more efficient way incorporate expression-level constraint priors into the sampling stage remains an open challenge in SR. Additional related work on SR methods is reviewed in Appendix

---

3 Methods
To resolve the inefficient search limitations in symbolic regression (SR) as discussed Section 2 we propose Graph-based Symbolic Regression (GSR) with hybrid neural-guided Monte Carlo tree search (hnMCTS) to address both redundant search space and incomplete constraint incorporation, with an improved convergence and exploitation and exploration trade-off.
3.1 Expression Graph Representation
Expression Graph. The redundant search space originates from the failure to capture the expression equivalencies and preserve the permutation invariance. Consequently, we propose the expression graphs (EGs), denoted as G, which are unique expression representations aggregating their equivalent expressions, regardless of the generating sequence they went through during the sampling.
Expression
Expression Graph
7h Feature X
Input:
(x,]
(cry)
T2
Opprration Arily !Type
Node Embecding
Oiee Hiet
İ.
Ep
RL
Cij's Feature di
Linear
h
Operand Relation
R-GON
CU
W2
Semantic Expression
exp(A) x (0) C
Ejene
Euep
TRS: ((e))(((488)) 1///) Xeep(() expl B) 8] Tx 11/.).1/C))
[[A. 1/B,C)
i.
(A x A (8 C)
Expression Graph
MLP
(4-8YCC
(4"C)"B
A C))
R-GON Module
Tohmss
Σ,
Layer
Output
Expression Tree
Expression Graph
Output:
(c)
Figure 2: (a) Comparisons between the semantic expression representation, the ET representation (blue block), and the EG representation (orange block, where bullet points are the TRS rewrite rules). (b) Examples of (i) node and edge features of an EG, and (ii) EGs to unify the binary operations and align the commutative operands at the same level for permutation invariant representations. example of R-GCN encoder architecture for our EG representation, where node j is the neighbor of node i, (i) represents all neighbors of with the operand relation a, @ represents the element-wise multiplication, Wo and W. are self-loop weight and edge-dependent weight respectively, Ci denotes the coefficient for W.. and denotes the activation function.
(c) An
As shown in Fig. 2 (a), the construction of an EG involves two procedures: (1) implementing the term rewriting system (TRS) with rewrite rules for equivalences, detailed in Appendix C.2 which captures the expression equivalences through recursively normalizing the current sampled candidate expression into its canonical form; and (2) converting the canonical form of the semantic expression structure, which intrinsically encodes the permutation
into a directed acyclic graph (DAG)- invariance through unification of the binary operations (e.g. (+.-) Σ and (x ID and aligning commutative operands at the same level, as shown in Fig. I and (b). Note that for each EG, the node set = {n,} represents sampled operators, and the edge set relationships between operators and operands. Compared with ET, EG with the DAG structure extended its capability for richer structure information (e.g., additional node features of arity and edge features of operand relation a) to differentiate operands' roles (e.g., base or power for A, + or forr >); and for improved scalability handling complicated and highly dimensional datasets by node-sharing of the same variable node. Theorem 3.T establishes the uniqueness of the EG representation, and Appendix C.2 an extended discussion of the EG representation.
{e,y) encodes
Theorem 3.1. Let T(T) be the set of expressions over 4 and R be the TRS defined in Appendix C.2 Terms & to E T(D) are equivalent (t1 ~R t2) if they rewrite to common term under ++ (reflexive-

---

mrnsitive-symmeet closure of R). EG T()) - G is a function normalized via bottom-up R application. The EG representation is unique if for all tn ~n tz. EG(t)) = EG((t) (isomorphism))
Proof Sketch. Sketch. For uniqueness, we show that R is a terminating and confluent TRS, implying unique normal forms (irreducible terms under -*)). The EG then encodes this normal form as a minimal DAG with shared substructures, preserving uniqueness. Appendix C.2 details the complete proof.
Graph Neural Network Encoder. Considering the DAG structure of EG and discrete edge fea- tures capturing the operands' relations, we propose Relational Graph Convolutional Network (R-GCN) [20]. a novel graph neural network (GNN) architecture [21 [2]. as the foundation encoder in GSR. Its edge type-dependent message passing not only preserves the permutation invariance for EG during sampling, more importantly, the unique message passing depending on different operand relations also creates more accurate and expressive representation than traditional encoders with ET. Figure c shows an example of the R-GCN architecture, and Appendix C.3 provides the explanations for the detailed message passing updates. Given an EG G as input, R-GCN encodes the
expression state as the EG representation, formally expressed as:
R-GCN(G(V.E)) = ( To). = (P,r),
(2)
where P | is prior probability matrix, with each row corresponding the prior probability distribution for each E V. and each column in the distribution representing the prior probability for each operator & b to be added to the node. RR the predicted reward value of the state G. Telp and |- are trainable models predicting P and r given G, with parameters 9 of R-GCN. By self-learning with hnMCTS, EG representations encode expression-level constraint priors into P and providing structure-aware and constraint-informed neural guidance during the sampling.
3.2 Hybrid Neural-guided Monte Carlo Tree Search
To facilitate the EG-encoded constraint prior with MCTS (8. we focus on neural-guided MCTS 14 23 which enables self-learning of R-GCN on EGs from MCTS. However, the pure neural-guided MCTS suffers the initial training bottlenecks due to the biased estimates of the under-trained R-GCN module. Consequently, we propose new hybrid neural-guided MCTS (hnMCTS), which adopts an adaptive e-UCB policy to dynamically combine the vanilla MCTS with neural-guided MCTS, to provide warm-start training, faster convergence, and the balanced exploitation and exploration.
MDP settings. We model the expression generation process as a finite-horizon MDP with the sampling trajectory {so, ag. $1, a1...., St, ....... OH-1, $H }. where H is the maximum com- plexity. At step t, we define the state 8, to be the current EG, G., with the sampling node; action a, to be the sampled operator variable á & added to 8t. When is complete (attaining the maximum complexity or finishing the expression in closed form), we obtain the finalized expression fr- Then get reward R(T) through evaluating with the metries in Appendix C.7
Figure 3 illustrates the workflow of hnMCTS. To sample T, at each step t in T, we perform one batch of hnMCTS for each state St to update MCTS policy Thr that we can sam- ple ag ~ (s)) During each iteration of the batch, hnMCTS samples a trajectory т = {s2, at s/-1.11. o..1..... s444:
H1-11 _# , HH-)) from the root state & as the search, where the superscript denotes the search step in 7, and subscript denotes the current complexity.
Steps in hnMCTS. As shown in Fig. 3 hnMCTS will switch between the vanilla and the neural- guided MCTS at the beginning of each T, by the e-UCB policy, with the random variable E [0, and the given threshold c, E (0, 1]. Then, it will take four steps based on the selected policy.
Starting from the root state st, hnMCTS selects the next step iteratively before arriving at an expandable node or terminal node following the maximum Upper Confidence Bounds (UCB) 24 or the maximum
Selection:
medictorr-bseed UCB (PUCB) policy based on the selection of the e-UC policy. Define b as the next possible action, Ct and Cp as exploration rate controlling hyperparameters, Q(s, a) as the expected action value, ((s, a) as prior probability, and N N(s, a) as the total number of
5

---

One Epoch Trajectory
Evaluate
R-GON
$2
St
Backpro-
a1 Ri
02 m
constant fitting
pagation
UCB <
Expand
Simulate nxe e.. rollout
Evaluate constant fitting
...
Q ('(s,),
Select E-UCB policy
Iterate
Backpropagste
max PUCC .... F 2,6
Expand
Simulste
R. R-GON Modill
Predict
[()))),()),
hnMCTS Module
neurat guidianos
One Iteration Trajectory
Non-unary
Unery operator
node
Uneamoled
Sampung
Figure 3: An exemplar workflow of hnMCTS with R-GCN on EG for neural-guided MCTS.
times for taking action a from state 8, we have the e-UCB policy:
c-UCB policy
arg max, (s, a) = ((s, CU (PUCB( ( Q(x + CpP(s C
In Σ, ((sb)) [ss.a)
if < Ce
Σ, N(a.)) if 2 C.
arg max
[+N(sa)
Expansion: If hnMCTS traverses to a visited node with unvisited children, we call it an expandable node. We will select one of its unvisited children for the following simulations.
Simulation: After traversing to an unvisited node, based on the selection of the e-UCB policy, hnMCTS will either have the uniform prior and calculate R(s) through random rollout as the vanilla MCTS; or obtain (P(s, [(ss)) (Te (s), To lr (s)) through the R-GCN encoding on EGs as the neural-guided MCTS, where ((s,-) -) is the prior distribution P and R(s) is the expected reward r of state & according to (2. which incorporates expression- level constraint priors directly into sampling process in PUCB 33) by the self-learn mechanism according to (4).
Backpropagation: After obtaining R(s) of the expandable node, hnMCTS will increase N(s, a) by one and update Q(s, (1/N(s,a) a)) Ex R(s') ) for all those (s, a) pairs that have been traversed through. Here, s' represents the next state of by taking action a.
After one batch bnMCTS, we can update Ths = N(s,, 3/ Es Nss b) to step next T. When is completed, subsequently update the parameters 0 of R-GCN by minimizing the loss function:
Σ((R() Tol- TM (8t) log Tolr(s:)).
Appendix Couuliiee the proposed hnMCTS algorithm. With the dynamically controlled c-UCB policy, hnMCTS provides: (1) expression-level constraints incorporated proactively via R-GCN priors (P,r) within PUCB 33. densifying rewards without evaluations on full expressions; (2) smooth warm-start for training, initially relying on random roll-outs to explore broadly, and progressively shifting to ceural-guided MCTS R-GCN training stabilizes; (3) the reduced reliance on costly constant fitting of evaluations via R-GCN predictions, significantly reducing computational overhead; (4) improved convergence established in Theorem 3.2 and balanced exploration-exploiiane trade-off by controlling c, for the c-UCB policy as explained in Appendix C.5.2
Theorem 3.2. Given bounded rewards R(s, a) E (0, 1), finite action space A. Ct > 0, Cp 0, unbiased simulations, and O < Ce S 1, hnMCTS with the e-UCB policy ensures limy Nr (8,a)
0o for all a E A. Nr(s, a) = O(In T) for a + a" (the optimal action), and Nr(s, a") T.
Proof Sketch. The c-UCB policy with c. S alternates between the vanilla (UCB) and the neural-guides (PUCB) MCTS, ensuring infinite exploration of all actions due to UCB's 's exploration term, even if the neural prior P(s,a* a²) 0. Suboptimal actions have logarithmic visit counts (E(Nr(ss.a) (In 7)) while the optimal action's visits grow linearly (E[NT(s, (s,[N(ss,") a²)) ~ T), driven by Q-values from unbiased simulations. This guarantees policy convergence (m (8, + with relaxed condition compared to neural-guided MCTS. Appendix C.5.I provides the complete proof.

---

4 Experiments
4.1 Synthetic Dataset Benchmarking
Benchmarking Datasets: we have evaluated our GSR on diverse benchmark datasets, including 1) the black-box dataset consisting of 120 tasks without solution expressions from PMLB [25]. 2) Feynman dataset [T]] consists of 119 physics-informed problems with ground-truth solutions. 3) Nguyen's 's SR benchmark dataset [26] and Nguyen constant dataset [7]. For the black-box dataset and the Feynman dataset, we follow the experimental settings of state-of-the-arg symbolic regression benchmarking SRBench [27]. which includes 21 baseline models for the black-box dataset benchmarking and 14 baseline models for the Feynman dataset benchmarking. Additionally, we compared with the recent wansformer-based MCTS method, Transformer- -based Planning for Symbolic Regression (TPSR) 10. in the black-box dataset; and a neural-guided GP method, Unified Framework for Deep Symbolic Regression (uDSR) 191. in the Feynman dataset; as well as transformer- GP method, Deep Generative Symbolic Regression (DGSR) [6], for both datasets. We summarize our additional experimental settings in Appendix D.1
Figure summarizes the benchmarking results on both the black-box and Feynman datasets. On the complex black-box problems, GSR consistently achieves higher R scores than most baselines (Fig. a)), indicating stronger predictive accuracy in settings where expression structures are not explcitly provided. On the Feynman dataset, GSR recovers a larger number of ground-truth expressions (Fig- (d). demonstrating improved ability to identify physically meaningful symbolic forms rather than merely fitting numerical outputs.
Benchmarking Performances:
In terms of model complexity (Fig. bb). GSR attains relatively compact expressions while main- taining competitive accuracy. Regarding computational efficiency, GSR exhibits training times comparable to the pretrained transformer-based TPSR method [T0] (Fig. 4c),, despite not rely- ing on any pretraining or large auxiliary datasets, suggesting that the canonical expression graph representation effectively reduces redundant structures during search.
Overall, these results indicate that encoding invariance and expression- level constraints directly into the search process improves the trade-off between accuracy, interpretability, expression simplicity, and computational cost, particularly on more challenging benchmarks. Additional results on the Nguyen benchmark reported in Appendix D.2
R² Test
Model
Training
TPSA
Solution Rate (%)
Operon
GSR
SBP-GP
DGSR
FEAT
UDSR
DGSR
AlFeynman
GP POMMGMMAA
AFP FE
LOBM
ITEA
hardomhoriesi
GP SOMEA
Adadbooot
Operon
DSR
AFP
opleam
APp
SBP-GP
FFX
EPLEX
KemelRldge
BSR
DSR
FEAT
MRGP
HLF
FFX
Unear
ESR
AlFoynman
-0.25 0.00
0.25 0.50 0.75 1.00
10
10P
10a
10f
(a)
(b)
(c)
(d)
Figure 4: Performance comparisons of 24 algorithms regarding: (a) R²2 test score, (b) model size (solution complexity), and (c) training time on the black-box dataset; Performance comparisons of 17 algorithms regarding: (d) solution rate (the ratio of ground-trure truth equivalent solutions in mathematics to the total equations in the dataset) on the Feynman dataset. The dots represent the median value, and the shadow line represents the 95% confidence interval.
4.2 Ablation Study
Ablation Settings: To evaluate the individual contributions of invariance and expression-level constraint encoding in GSR, we conduct ablation studies six challenging symbolic regression
7

---

Model
GSR
MCTS-EG
hnMCTS-ET
MCTS
Invariance Encoding
Expression-lvvel Constraint Encoding
Solution Rate (%)
98.0
90.3
92.7
85,8
Number of Evaluations.
89,128
107.651
120,056
162.783
Training time (s)
13302.9
2089.2
1759.3
2782.6
Sparse Rate (%)
11.4
21.2
16.5
25.3
Table 1: Comparisons between GSR and three ablation models regarding the average solution rate, number of evaluations, the training time for the first solution, and the sparse rate (ratio of invalid generated expressions to the total expressions generated in an epoch) at 50% of the total training epochs on six tasks (Nguyen- 12 and five Nguyen benchmark with constant tasks).
tasks drawn from the Nguyen benchmark with constants 7 and the Nguyen- 12 task from the original
Nguyen benchmark [26], under consistent experimental settings described in Appendix D.2 In
addition GSR, we compare the following ablation models: 1) MCTS-EG employing the EG representation for invariance encoding but using standard MCTS without constraint encoding; 2)
hnMCTS-ET, using hnMCTS for constraint encoding but relying on representation without invariance encoding; 3) MCTS, a baseline using vanilla MCTS with ET, lacking both invariance and constraint encoding.
Performance Analysis: Table reports the results of the ablation study. Relative to vanilla MCTS, hnMCTS-ET, which incorporates expression-level constraint encoding, achieves higher solution rates, improved training efficiency, and reduced reward sparsity. These gains indicate that constraint -uudde sampling. together with the expressive capacity of the R-GCN encoder, provides more informative guidance during search, particularly for tasks involving constants and complex constraint structures.
In contrast, MCTS-EG requires substantially fewer function evaluations than vanilla MCTS, with a 33.9% reduction in the effective search space, demonstrating that equivalence-avee awaaae expression graph representations reduce redundant exploration through state-space compression. However, its overall training time remains higher than hnMCTS-ET, as constant fitting must still be performed explicitly during search, limiting the computational savings from reduced evaluations.
By integrating both expression based invariance encoding and constraint- -uidde neural search, GSR achieves further improvements in solution rate, evaluation efficiency, and reward sparsity mitigation. This result highlights the complementary roles invariance-awar aware state compression and expression-level constraint guidance, enabled by edge-aware message passing of the R-GCN over structurally enriched expression graphs. Additional analyses of reward sparsity and supporting visualizations are provided in Appendix D.4
5 Application
To demonstrate GSR's scalability and constraint- aware capability in highly complex, irregular, and physics-constrained real-world settings, in this section, we apply to a materials science application focused on predicting interatomic potential energy for atomistic simulations. This problem is central to understanding material behavior and accelerating the design and discovery of novel materials.
While traditional first-principles methods such as density functional theory (DFT) 28 29], offer high accuracy, their use is often limited by the substantial computational cost and extensive memory requirements. To address this, surrogate ML models increasingly explored, with symbolic regression showing strong potential for balancing efficiency and interpretability. Appendix [E. details the problem backgrounds, and Appendix E.2 describes the practical usage of GSR in this domain.
5.1 Problem Setting
We adopt dataset from first firttpprieiples molecular dynamics simulation of 32 copper atoms from [16].
which is described in detail in Appendix E.4 We aim to discover the interatomic potential energy
function f()) that effectively maps the atoms' pairwise distances to total formation energy E.
In this study, we consider four methods to learn the interatomic potential energy function: 1)
GNN-based blx method, CGCNN (30]: 2) genetic programming- SR method from [6:
3) our GSR in this paper, and 4) the vanilla MCTS method. GP1 and GP2 in Table 2 are two fitted
8

---

expressions by genetic programming reported in [16]. Appendix E.4 summarizes our experimental settings for baseline models.
5.2 Physics Constraints
To ensure the resulting f(-) is physically meaningful, we consider the following additional physics constraints besides the universal constraints listed in Section 2.3
(1) Scalar Output The output of (() has to be single scalar to match E for valid evaluation.
(2) Unit Consistency A scalar const should be multiplied before any polynomial terms and in any number in the exponential terms to ensure the unit calculation aligns with the physical meaning.
(3) Electrostatic Repulsion : The energy approaches infinity at an infinitesimal distance (i.e. f(r oo) due to Coulomb's law.
We include explanations for all constraints in Appendix F..3 Note that CGCNN only satisfies the scalar output (I) constraint; GPI and GP2 satisfy both the scalar output (1) and unit consistency (2) constraints; and MCTS and GSR satisfy all the constraints.
5.3 Performance and Transferbility
Test Results: Table 2 presents our comparison results for the formation energy prediction, where the Mean Absolute Error (MAE) is computed based on linear regression results of the model expressions regarding the formation energy E, divided by the number of atoms (32) in the crystal. Figure 5 depicts curves of the interatomic potential energy functions f()) according to the learned analytical expressions in Table 2 It illustrates the electrostatic repulsion that the energy goes to infinity
f(r) + oo short interatomic distance (r - 0), with more details described in Appendix E.6
Comparing training and testing MAEs in Table GSR achieves the best testing MAE while preserving all physical constraints (the 1/r"" terms reflect electrostatic repulsion at an infinitesimal distance). MCTS, though satisfying all the constraints as GSR, its largest MAF and highest complexity highlight the impact of sparse rewards to vanilla MCTS during random rollouts. In contrast, GSR guides MCTS with encoded physics constraints that improve the likelihood of generating valid expressions during simulations. Although CGCNN has the lowest training MAE, its larger testing MAE indicates the potential overfitting risk, showing less robustness than SR methods such as GSR and GP2 in liti-limited scenarios. GP1's 's low complexity comes at the cost of significantly poor performance,
MAE (meV/ato)
Model
Constraints
Expression (r) (eV/atom)
Complexity
Train Test
Transfer
x
2.08
CGCNN
44.89
x
3.09
(2)
11 0.211))(()) +0.97( 0.33"*(r))--
3.68
GPI
43.32
(3)x
3.53
(2
7.33 3.98 s(r) (27.32
2.57
GP2
(3)x
(11.13
0.08, (oroa)rrrrrrr))))) ((r)) (())
28
2.70
41.63
(2
2113.701 + 27.18/ 13.70re 0.08.2
3.91
MCTS
(3V
+ le 12.87r
3.63
36.31
(2
G6.04 10 10-" -11)
2.41
GSR
(3)
Σ 121.41(r +
24
2.63
34.15
Table 2: Analytical interatomic potential energy functions frr fitted to DFT total energy using GP1, GP2, MCTS, and GSR compared with the CGCNN model. Here, ((r is in the unit of eV eV/atom, and s()) is the smoothing function introduced in Eq (7) of [6]. "Constraints" denotes the satisfied constraints defined in Section 5.2 "Complexity" is computed by the of nodes in the expression graph or the expression tree. "Train/Test" represents the training and testing MAE of the copper dataset. "Transfer" represents the MAE of expressions directly tested on the newly generated dataset as described in Appendix E.7 Typically, DFT achieves an MAE of 10 ~ 20 me V/atom for energy predictions for reference.

---

1o3
GSR
1e1
MCTS
1e3
GP1
1e3 GP2
10
20
20
20
10
10
f(()
10
Figure 5: Interatomic potential energy functions f(r) learned by GSR, MCTS, GPI and GP2, as a function of interatomic distance r ranging between to 5Å. The z-axis is the interatomic distance and the y-axis is the interatomic potential energy. Only MCTS and GSR produce physically meaningful results, adhering to the strong electrostatic repulsion at short interatomic distances, while both GP1 and GP2 fail. A finite-valued model short interatomic distances, like GP1 and GP2, allows atoms to overlap. which can either cause downstream dynamics simulations to crash or lead to meaningless results violating physics laws.
while GP2, despite the comparable predictive power to GSR, has higher complexity that makes the underlying dependencies less straightforward. Besides, none of CGCNN, GP1, or GP2 adheres to the electrostatic repulsion constraint, resulting in less physically meaningful solutions short interatomic distances. Overall, GSR demonstrates strong expressive power with lower complexity than GP methods and reduced overfitting compared neural networks, while consistently maintaining physical meaning. Plots for Table results are included in Appendix B.5
To further assess the transferability of the interatomic potential energy model, we create a new dataset of 100 samples by performing DFT-based molecular dynamics simulations where volumetric compression of approximately 50% is applied to the unit cell of the original copper dataset. As a result, the new dataset consists samples with shorter interatomic distances (bond lengths) between neighboring atoms, described in Appendix E.7 This allows to assess the transferability of our physics-constrained analytical model compared to other conventional machine learning models. The "Transfer" MAE in Table 2 is obtained through zero-shot learning of the solution expressions directly tested on the newly generated dataset, where GSR and MCTS show significantly lower MAEs. As electrostatic repulsion has a direct impact at shorter interatomic distances, GSR and MCTS that integrate the physical constraints provide better generalization to shorter bond lengths on the new dataset. These results again verify that the incorporated physics constraints can prevent overfitting and foster the discovery of generalizable knowledge.
Transfer Results:
6 Conclusion and Limitation
In this work, identified two critical challenges in current Symbolic Regression (SR) approaches: (I) redundant expression representations that fail to capture expression equivalences and preserve the permutation invariance. and (2) the absence of expression-level constraint priors incorporated during the sampling stage, leading to low search efficiency, reward sparsity, and limited applicability in real-world settings. To address these challenges, (GSR) a method built on Relational Graph Convolutional Network (R-GCN) operating over expression graphs (EGs); and hybrid neural-guided Monte Carlo Tree Search (hnMCTS) framework with an adaptive c-UCB policy. EGs provide invariance-aware representations and higher-order structural features, enabling R-GCN encode constraint priors more effectively. hnMCTS with R-GCN encoding facilitates a balanced exploration-explooitaioo trade-off, reduces reward sparsity, and enhances search efficiency.
proposed Graph-based Symbolic Regression
Comprehensive benchmarks and ablation studies on synthetic datasets confirm that GSR achieves superior performance with compressed search space. Furthermore, real-world application in materials science, involving complex physics- based constraints, demonstrates GSR's 's capability to produce interpretable and physically consistent solutions, underscoring its practical value.
Limitations remain for the GSR, where the term rewriting system (TRS) is a handcrafted rule set that captures most of the basic algebraic equivalences but not higher-order symmetries or domain-specine invariances. Looking forward, extending GSR with automated learned TRS rules and pretrained symmetry-awand graph encoders could further enhance generalization and cross-domain applicability.
10

---

Acknowledgment
Z. X. and X.N. Q. acknowledge the support from U.S. National Science Foundation (NSF) grants DMREF-2119103, SHF-2215573, and IIS-2212419. First-principles calculations by K.A. were supported by the Center for Reconfigurable Electronic Materials Inspired by Nonlinear Dynamics (reMIND), an Energy Frontier Research Center funded by the Department of Energy under award DE-SC0023353. X.F.Q. acknowledges the support from NSF under grant CMMI--22690B Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.
References
[1] Dimitrios Angelis, Filippos Sofos, and Theodoros E Karakasidis. Artificial intelligence in physical sciences: Symbolic regression trends and perspectives. Archives of Computational Methods in Engineering, 3066):3845-3865, 2023.
[2] Nour Makke and Sanjay Chawla. Interpretable scientific discovery with symbolic regression: review. Artificial Intelligence Review, 57(1):2, 2024.
[3] Yiqun Wang. Nicholas Wagner, and James M Rondinelli. Symbolic regression in materials science. MRS Communications, 9(3):793-80S 2019.
[4] Jinglu Song, Qiang Lu, Bozhou Tian, Jingwen Zhang, Jake Luo, and Zhiguang Wang. Prove symbolic regression is NP-hard by symbol graph. CoRR, abs/2404. 13820, 2024.
[5] Marco Virgolin and Solon P Pissis. Symbolic regression is NP-hard. Transactions on Machine Learning Research, 2022.
[6] Samuel Holt, Zhaozhi Qian, and Mihaela van der Schaar. Deep generative symbolic regression. In The Eleventh International Conference Learning Representations, 2023.
[7] Brenden K Petersen, Mikel Landajuela Larma, Terrell N Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-secking policy gradients. In International Conference Learning Representations, 2020.
[8] Fangzheng Sun, Yang Liu, Jian-Xun Wang, and Hao Sun. Symbolic physics learner: Discovering governing equations via Monte Carlo Tree Search. In The Eleventh International Conference on Learning Representations, 2023.
[9] Mikel Landajuela, Chak Shing Lee, Jiachen Yang, Ruben Glatt, Claudio P Santiago, Ignacio Aravena, Terrell Mundhenk, Garrett Mulcahy, and Brenden K Petersen. unified framework for deep symbolic regression. Advances in Neural Information Processing Systems, 35:33985- 33998, 2022.
[10] Parshin Shojaee, Kazem Meidani, Amir Barati Farimani, and Chandan Reddy. Transformer- based planning for symbolic regression. Advances in Neural Information Processing Systems, 36:45907-45919, 2023.
[11] Silviu-Marian Udrescu and Max Tegmark. AI Feynman: A physics-inspired method for symbolic regression. Science Advances, 6 6(16).ea4 2020.
2020.
[12] Liron Simon Keren, Alex Liberzon, and Teddy Lazebnik. A computational framework for physics-informed symbolic regression with straightforward integration of domain knowledge. Scientific Reports, 13(1):1249, 2023.
[13) Pierre-Alexandre Kamienny, Guillaume Lample, Sylvain Lamprier, and Marco Virgolin. Deep In International Conference on
generative symbolic regression with Conte-Carlo-Tree-SS Machine Learning, pages 15655-15668. PMLR, 2023.
[14] David Silver, Julian Schrittwisser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, al. Mastering the game of GO without human knowledge. nature, 550(7666):354-355, 2017.

---

[15] Terrell N. Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio P. Santiago, Daniel faissol, and Brenden K. Petersen. Symbolic regression via deep reinforcement learning enhanced genetic programming seeding. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
[16] Alberto Hernandez, Adarsh Balasubramanian, Fenglin Yuan, Simon AM Mason, and Tim Mueller. Fast, accurate, and transferable many-body interatomic potentials by symbolic regres- sion. npj Computational Materials, 5(1):1 12, 2019.
[17] Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Al Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Neural Information Processing Systems, 15:4860-48711 2020.
[18] Paul Kahlmeyer, Joachim Giesen, Michael Habeck, and Henrik Voigt. Scaling up unbiased searchbased symbolic regression. In Proceedings of the Thirty-Thar International Joint Confer- ence on Artificial Intelligence, IJCAI-24, pages 4274 2024.
[19] Wassim Tenachi, Rodrigo Ibata, and Foivos Diakogiannis. Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws. The Astrophysical Journal, 959(2):99, 2023.
[20] Michael Schlichtkrull, Thomas N Kipf. Peter Bloem, Rianne Van Den Berg, Ivan Titov, and Max Welling. Modeling relational data with graph convolutional networks. In The semantic web: 15th international conference, ESWC 2018, Heraklion, Crete, Greece, June 3-7, 2018, proceedings 15, pages 593-607. Springer, 2018
[21] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip. A comprehensive survey on graph neural networks. IEEE transactions on neural networks and learning systems, 32(1 ):4 24, 2020.
[22] Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In International Conference on Learning Representations, 2017.
[23] Surag Nair. A simple Alpha(Go) Zero tutorial, 2017.
[24] Peter Auer. Using confidence bounds for exploitation-exploraiio trade-offs. Journal of Machine Learning Research, 3(Nov):397-422, 2002.
[25] Randal S Olson, William La Cava, Patryk Orzechowski, Ryan J Urbanowicz, and Jason H Moore. PMLB: a large benchmark suite for machine learning evaluation and comparison. BioData mining, 10:1-13, 2017.
[26] Nguyen Quang Uy, Nguyen Xuan Hoai, Michael O'Neill, Robert McKay, and Edgar Galván- López. Semantically-basee crossover in genetic programming: application to real-valued symbolic regression. Genetic Programming and Evolvable Machines, 12:91-119, 2011.
[27] William La Cava, Bogdan Burlacu, Marco Virgolin, Michael Kommenda, Patryk Orzechowski, Fabrício Olivetti de França, Ying Jin, and Jason Moore. Contemporary symbolic regression methods and their relative performance. Advances in neural information processing systems, 2021(DB1):1. 2021.
[28] Pierre Hohenberg and Walter Kohn. 136(3B):B864, 1964.
Inhomogeneous electron gas. Physical review,
[29] W. Kohn and L. Sham. Self-consistent equations including exchange and correlation effects. Phys. Rev., 130-A1133-A1144 1965.
[30] Tian Xie and Jeffrey C Grossman. Crystal graph convolutional neural networks for accurate and interpretable prediction of material properties. Physical Review Letters, 120 (0011):1443001 2018.
[31] G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave wave basis set. Phys. Rev. B, 54:11169-1111 1996.
12

---

[32] John R Koza. Genetic programming as means for programming computers by natural selection. Statistics and computing. 4:87-112, 1994.
[33] Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. Science, 324(5923):81-88 2009.
[34] Runhai Ouyang, Stefano Curtarolo, Emre Ahmetcik, Matthias Scheffler, and Luca M Ghir- inghelli. SISSO: A compressed sensinng method for identifying the best fow-dimensional descriptor in immensity of offered candidates. Physical Review Materials, 288):083802, 2018.
[35] Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3922-3337 2016.
[36] Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. SymbolieGPT: A generative transformer model for symbolic regression. arXiv preprint arXiv:2106. 14131, 2021.
[37] Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambatista Paras- candolo. Neural symbolic regression that scales. In International Conference on Machine Learning, pages 936 945. Pmlr, 2021.
[38] Maciej Swiechowski, Konrad Godlewski, Bartosz Sawicki, and Jacek Mañdziuk. Monte Carlo Tree Search: A review of recent modifications and applications. Artificial Intelligence Review, 56(3))2497-2562 2023.
[39] Bruno Buchberger. Bruno Buchberger's PhD thesis 1965: An algorithm for finding the basis elements of the residue class ring of a zero dimensional polynomial ideal. Journal of symbolic computation, 415--1):45--11, 2006.
[40] Gérard Huet. Confluent reductions: Abstract properties and applications to term rewriting systems: Abstract properties and applications to term rewriting systems. Journal of the ACM (JACM), 27(4) 821, 1980.
[41] Peter Auer, Nicolo Cesa-Bianche and Paul Fischer. Finite-time analysis of the multiarmed bandit problem. Machine learning, 47:235-256, 2002.
[42] Yoshitomo Matsubara, Naoya Chiba, Ryo Igarashi, and Yoshitaka Ushiku. Rethinking symbolic regression datasets and benchmarks for scientific discovery. Journal of Data-centric Machine Learning Research, 2024.
[43] Bowen Deng, Peichen Zhong, KyuJung Jun, Janosh Riebesell, Kevin Han, Christopher J Bartel, and Gerbrand Ceder. CHGNet as a pretrained universal neural network potential for infgge-i atomistic modelling. Nature Machine Intelligence, 10)):1031-1041, 2023.
[44] Songling Liu, Yonghao Sun, Zhaoyuan Liu, and Weihua Wang. Tetrahedral stacking of Copper its incubation period of crystallization: study of large-scale molecular dynamics simulation. The Journal of Physical Chemistry C, (21(44)21116-2188) 2023.
[45] Makov and Mike Payne. Periodic boundary conditions in ab initio calculations. Physical Review B, 51(7):4014 1995.
[46] P. E. Blöchl. Projector augmented-wave wave method. Phys. Rev. B, 50: 10:1799 79533-179777 1994.
[47] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made simple. Phys. Rev. Lett., 77:386538868, 1996.
[48] Hendrik J. Monkhorst and James D. Pack. Special points for Brillouin-zone integrations. Phys. Rev. B, 18:5188-5192, 1976.
[49] Chad A McCoy, Marcus D Knudson, and Seth Root. Absolute measurement of the Hugoniot and sound velocity of liquid Copper at multimegabar pressures. Physical Review B, 96(17):174109, 2017.
[50] DE Fratanduono, RF Smith, SJ Ali, DG Braun, A Fernandez-Pañelll,, S Zhang, RG Kraus, Coppari, JM McNaney, MC Marshall, et Probing the solid phase of noble metal Copper at terapascal conditions. Physical Review Letters, 12151)0 2020.
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: As detailed Section 3 and our main claims made in the abstract and intro-
duction our GSR's effectiveness and efficiency accurately reflect the paper's 's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Section Conclusions and Future Work briefly discussed our future research di- rection on addressing challenges in integrating more expression symmetries and a pretrained GNN encoder, which hopefully has stronger trustworthy guarantees on learned SR models.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged
the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
14

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: We provide complete proof of Theorem 3.1 in Appendix C.2 and complete proof of Theorem 3.2 in Appendix C.5.I
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof' can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have made all efforts to clearly describe our experimental settings in Section and [5] as well as in Appendix D.I C.6 and E..4] Multiple SR experimental runs have been conducted in Section and and their results have been explicitly reported.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to model. In general. releasing code and data often one good way to accomplish this, but reproducibilityy can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open- open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibiliity In the case of closed-source models, it may be that access to the model is limited in
15

---

some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Both the code and data will be released in an open-access repository after the paper is publicly available.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. ec/
public/guidess (eegli//uuudes//ooeebbbssssssss onlliicy)) for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
//nips. //nips.c/ppie http:liiigguddde//deeeeeeee
semmiissionni ssionPolicy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: We have made all efforts to clearly describe our experimental settings in Section 4 and 5 as well in Appendix D.I
C.6 and
E.4
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We have conducted experiments on the standard benchmarking dataset several times and reported their 95% confidence interval with the median value in Section
Guidelines:
The answer NA means that the paper does not include experiments.
16

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The computing platform together with runtime is described in the Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics httph:ttp/:/
Answer: [Yes]
Justification: We confirm that the research conducted and presented in the paper conforms, in every respect, with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there special consid- eration due to laws or regulations in their jurisdiciio))
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We have provided a brief statement on potential societal impacts, especially related to materials science, in Section
17

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., denloymemt of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: We do not foresee any potential risk of releasing related data and code we focus on fundamental research in Al for Science and follow the commonly adopted practice in the field.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We have cited the references for all the original code, data, and models used in our paper.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
18

---

The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g.. website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We have clearly described all the details of our GSR and the dataset we have generated via DFT in Sections] 5 and Appendix
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mooel part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The presented work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The presented work does not involve the risks listed above.
19

---

Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The presented work does not involve important original core methods originating from LLMs.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non non-standard components.
Please refer to our LLM policy for what should or should not be described.
20

---

Appendix
A
Computing Resources
We have performed our experiments on a platform with one CPU, Intel Xeon 6248R (Cascade Lake), 3.0GHz, 24-core, and one GPU, NVIDIA A100 40GB GPU accelcrator. For the reported Synthetic Dataset Benchmarking in Section 4 in the main text, our graph-based symbolic regression (GSR) can sample 60 expressions per second on average. For the materials science application in Section 5 GSR on average samples 20 expressions per second.
First-principles DFT calculations of the test dataset were performed on a compute node consisting of two Intel Xeon 6248R (Cascade Lake) CPUs with a total of 48 cores and 384GB DDR4 memory. The initial supercell structure of 32-atom fee copper was relaxed using VASP 31 with the conjugate gradient method in four steps within total of 18 seconds. For each temperature, VASP-based molecular dynamics simulations in the NVT ensemble were run for 6,000 steps for total of 900 minutes.
B Additional Literature Review on Symbolic Regression Methods
Symbolic regression (SR) has a rich history with diverse range of methods. Early ones primarily rely on Genetic Programming (GP)
which searches for candidate expressions through evolutionary processes such as selection, mutation, and crossover. Similar as brute- brute-force methods like SISSO 34 and basis-dependent methods like SINDy [35) GP-based methods often suffer from
16 32
poor scalability and high sensitivity to hyperparameters 7
Modern SR methods integrate deep learning (DL) and reinforcement learning (RL) to improve heuristic search efficiency. DSR [7] employs a recurrent neural network (RNN) with risk-seeking policy gradients (PG), performing well on simple tasks but incurring high computational costs due to repeated RNN inference for every sample and limited exploration in complex problems. Neural-guided genetic programming methods ( [ address these issues by combining RL with GP. improving both efficiency and exploration. Deep generative models 6 I0 36 37 offer fast inference but lack adaptability to out-of-distribution datasets due to their static pre-training. In contrast, Monte Carlo Tree Search (MCTS)-basee methods
X
13
38 have better exploration-exploitaiio balance,
achieving state-of-the-art performance.
C Methodological Detail
C.1 Transplantation and Constant Optimization
Transplantation: Inspired by the cross-over mechanism in Genetic Programming (GP), we introduce the concept of function modularity to break down complex problems into smaller sub-problems Leveraging a divide-and-come cnd-uuque heuristic search strategy, we store promising expressions with lower complexity than the current sampling average in budget set. These stored expressions can then be directly inserted into newly generated expressions as a single operator, which we refer to as the "transplantation" operator, denoted by Irans. The Irans operator functions as a variable that GSR can choose during exploration. When selected, GSR randomly picks a subset of expressions from the budget and uses them to generate new expressions, selecting the best one to calculate the action value of the trans operator. By default, GSR enables the trans operator after half of the total epochs, with the budget set to store the top models that have complexity lower than half of the maximum complexity.
Constants Optimization: Constants play a crucial role in symbolic regression, not only for main- taining unit consistency in real-world scientific problems but also for optimizing regression accuracy. In our approach, constants are introduced into expressions via the const operator, functioning as variable operator. While enabling constant optimization enhances the accuracy of the regression, it significantly increases the computational cost of evaluating an expression. To mitigate this, all constants within an expression are optimized only once per evaluation using a non-linear regression method, specifically the BFGS optimizer, as described in [7].
21

---

C.2 EG Representation Detail
EG as State. To leverage the expression graph (EG) representation G as the input for GSR encoding. we can directly feed 9 into a Relational Graph Convolutional Network (R-GCN) due the R-GCN's mechanisms, which naturally preserve the
inherent aggregation and cdge-dependent message-passiieg permutation invariance encoded EG and exploit the discrete edge type in EGs. However, to use EG as state in Monte Carlo Tree Search (MCTS), we first convert & into a canonical form by sorting the nodes. Next, we modify the upper triangle of the sorted adjacency matrix by replacing elements with the corresponding edge features indexed by numbers, where "0" denotes disconnected edge. The diagonal elements are replaced by the node features. Finally, we output this modified sorted adjacency matrix in its Voigt form, ensuring the representation of states while preserving symmetries.
Rewrite Rule
Level
Equivalence
a b+ b)
- b- (0,
Commutative/
Operation
x b II
[a,
Associative equivalence.
a b+ [1 1/b)
x (b+ + c) + Σ a x b, c)
Operation
Distributive equivalence.
(a') + (l(hcc
Operation
Power equivalence.
exp(a) x exp(b) + exp(a
Exponential/
Operation
log(a) + log(b) + log x b)
Logarithmic equivalence.
a + a + a x +
Expression
Identity equivalence.
a a+0 +
(a) + (a) +
Expression
Trigonometric equivalence.
a (a - (a + b) (a + b)2 a² + 2ab + 622
Expression
Square equivalence.
Table 3: An exemplar the TRS. Operation-level equivalences are judged through local operations without specific operand requirements, while expression-leve equivalences require specific operands to satisfy the equivalence.
Term Rewriting System (TRS). Table presents presentss example of the term rewriting system (TRS) with the rewrite rules and the corresponding equivalences for our expression graph (EG) representation. Equivalences categorized into operation and expression-level ones, depending on whether they can be evaluated using local operations without specific operand requirements. Operation-level equivalences, such as the commutative property + b b a, are identified solely based on the operator (e.g., addition or multiplication) and are encoded in EG by unifying operators into a single node type (e.g.. Σ or ID with unordered edge connections, preserving permutation invariance. In contrast, expression-level equivalences, such the trigonometric identity sin? (a) + cos? (a) = require specific operands to be satisfied (e.g., sin? h (a) + cos? (b) # 1) and are encoded only after the relevant operands are sampled, allowing EG to transform the expression into a simplified node structure. EG leverages node and edge features directly encode operation-lever equivalences during operator sampling, while expression-levv equivalences are deferred until the complete expression is formed, enabling verification and transformation. Due to their specific operand dependency, expression-level equivalences occur less frequently than operation-level ones, resulting in operation- level equivalences dominating the reduction of redundant representations during the search process.
Proof of Theorem 3.1
Uniqueness requires R convergent (confluent + terminating), ensuring unique normal forms nf (t), and EG encodes the normal form uniquely.
Termination: Finite expressions yield finite steps; R satisfies Termination as nules are non- recursive (no loops recreating left-hand sides).
22

---

Confluence: If - and f - v, then there exists w such that u -² w and v - w. R is confluent, as rules are orthogonal (left-linear, non-overlapping). Arithmetic rules (commu- mirror polynomial rewriting to sum-of-pre -fn-of-products form, confluent with multisets for commutative operands [39]. Trigonometriceexppeenti
rules do overlap with arithmetic rules, preserving confluence via modularity [40]. By Newman's Lemma, termination and local confluence imply global confluence, so It has unique (t).
FG Encodes Uniquely: Builds minimal DAG from n²((() (Appendix C.2): flattens chains (multi-operand nodes with multisets), expands distributivity, simplifies identiiiestrigonometii. Fixed variables ensure uniqueness up to isomorphism; node sharing for minimality, multisets for order-independence.
Consequently, t1 ~R t2 = nf(ti) = nf((t) = EG( = EG(tz). It also holds for nested cases (e.g.. (a (b+ c)) d) via recursive normalization.
C.3 R-GCN Message Passing Detail
We adopt the Relational Graph Convolutional Network (R-GCN) 20 for GSR, and modify it with the residual design for stable training with the following message passing equation for each R-GCN unit shown in Fig. 2
h42 = h4 + o EE Wah + Woh.).
(5)
where node j is the neighbor of node i, (i) represents all neighbors of i with the corresponding relation a, Wo and W. are self-loop weight and relation weight respectively, Cıj denotes the normal- izer for Wa. and a is the softplus activation function. To output the prior distribution P. we use multilayer perceptron (MLP) with softmax activations the node messages, while separate MLP with the softplus activation is employed to produce the global readout on node messages for predicted rewards.
The relation weight Wa is edge-feature dependent, which means that for different edge types, R-GCN will go through distinct message passing updates, leading unique expression representation for different operands' relations and more accurate and expressive encoding for expressions.
C.4 Hybrid neural-guided MCTS (hnMCTS)
Algorithm summarizes the pseudo-code for our GSR. The core of the algorithm the repeated sampling of expression trajectories 7 until either the desired solution is found or the maximum number of iterations is reached. At each trajectory step f, batch of hnMCTS with index and maximum batch size B is executed, following the four steps for hnMCTS iterations in Algorithm 2
Algorithm 1 Graah-bbsse based Symbolic Regression
Input: operator dictionary 4, Batch size B, maximum complexity H, e-UCB policy threshold e, repeat
for = 1 to H do for 1: = 1 to B do Update Q(s, a), ((, from hnMCTS(s,, (s.. c, end for
Update (st) N(s,, -)/ E,
Sample a; ~ (sr))s (s:), E $ grow with to 81+1
if s141 complete then
break end if end for Calculate R(T), record (x))) and (s)
Calculate the loss / ((R(T) (s))())))² (og Tolp(s:))
Update F. according to until Optimal f is found Output: solution expression f
23

---

Algorithm 2 hnMCTS
Input: root state St, value function Q(s, a), visit count (s, operator dictionary 4, maximum complexity H. e-UCB policy threshold Ce
Sample [0. 1)
for to H - t - do
if C. then
Select "++ arg max, UCB(s a(.)
else
Select a'+i arg max, P a/+4)
end if
Grow s'+; with a'+.. step to i+1
$14141 and update N (1)+21 a(+1) + (844)))444) +
if $4i41 has unvisited children then break end if end for
i++
if s|4; is expandable then
if <<c. c, then
repeat
Grow 8445 with a random مي_ي. step to #ffit1 and update 4- i
until 8:4; is complete
Evaluate R(s)+4) with constants fitting and record P(st,-) = 1/1$1
else
R(st) = nel (s), P(s).)) ne (s))
end if
else
Evaluate R(s(4;) with constants fitting
end if
Backpropagate Q(s, a) (1/N a)) Es) R(s') for all (s, pairs traversed through Output: value function (s, a), visit count (ss
C,5 Proof of Convergence and Exploration-Exploittiio Trade-Off
This appendix provides detailed proofs for the convergence (Theorem 3.2 and exploration- exploitation trade-off of the hybrid neural-guided Monte Carlo Tree Search (hnMCTS) with an c-UCB policy, as outlined the main text. We first prove convergence for the edge cases (ce = 1, 0) and the general case (0 < 1), showing that hnMCTS achieves better convergence rate than pure neural- -uiddl- MCTS. We then analyze the exploration-expllitttin trade-off, highlighting how the c-UCB policy balances these objectives.
C.5.1 Convergence
We consider state with a finite action space A, where |A| K. Rewards are bounded, R(s, a) [0, 1). and the optimal action is a* arg max.EA (* a), with Q* (s, a) as the true action value. The hnMCTS policy is defined as:
Nr(8,a)
TM (s,
EsEA Nr(s,b)'
where Nr a) the visit count for action after T iterations. Convergence requires:
[1 if =" if a a*
lim (s, - -T (s, T-00
0
assuming a unique optimal action for simplicity. For multiple optimal actions, the policy distributes probability equally among them.

---

The c-UCB policy selects vanilla MCTS with probability and neural- MCTS with probability 1 Cc. as defined in Equation (1) of the main text. Naive MCTS uses the UCB selection policy:
In E, (s,
UCB(s, a) = Qr (s. + V
(T)
where Qr(s, a) Nr(0,0 E Re(8, a), and Cy 0. Neural-guided MCTS uses the PUCB selection policy:
Nr, PUCB(s,a) = Qr(s, CCP(ss,)
1+ NN((s,))
where P(s,0 a) is the neural network's ' prior probability, Ea P(s,a) 1, P(s,a) a) 2 0. and Cp > 0.
Case 1: Pure Naive MCTS (c. =
Theorem C.1. Given bounded rewards R(s, a) (0, 1), a finite action space A, and Cy > 0, MCTS with the arg max, UCB selection policy ensures limy- + (r a) = for all E A.
Proof. Assume there exists an action a E A such that (r a) S M < oo for all T. After some iteration To. Nr(s, = m S M is fixed. The UCB value is:
IT
UCB(s, Qr(s, a) + +Cv
m
Since Qr(s,a) E [0, 1), and T = E2 Ny (s. b). as 00, In S0:
UCB(ss, a) + c0.
For another action a' # with (r) + 0o:
In
UCB(s,d') = Qr(s,a') + Cup
(s.
As Nsr + x., UCB(ss d'), forcing selection of a, contradicting the assumption that Nr a) = m. Consequently, lim (r(ss oo.
Nr(a.a") + 0, so UCB(s, - Qr a') S 1. For large T, (s)
In
Theorem C.2. Given the assumptions of Theorem unbiased random rollouts, and unique optimal action a", we have (r a) = ((n for a a" and Nr a²) ~ T.
Proof. For a suboptimal action + a", define A, (s, (s) > 0. Action a is selected when:
IT
In
Qr(s, + Cu
Qr(s, + Cu
NT(s,a) As T + xo, Qr(s,a) + as) Qr(s,a") + (²(ssa)))) The exploration term must overcome A.:
Nr(s,a*)
IN
Cu
Nr(8, 2 A..
Squaring gives:
C2 In
Nr(8, <
A2
From UCB analysis
41 Nr(8, a) s
8ln T A2
constants, so (s) = O(n For K
suboptimal actions:
Σ Nr(s, O(KK 11(K)))
agta
Since E. ar(
Nr(3,a")=TT
= - T.
asa-
25

---

Remark 1, Theorems C.1 and C.2 imply (s,0a) + 1, (s) + 0 for o a", ensuring
convergence the optimal policy.
Case 2: Pure Neural-Guiee MCTS (c, =0)
Theorem C.3. Given the assumptions of Theorem C.I and Cp 0, MCTS with the arg max, PUCB policy ensures lim7-+00 (r(s, = oo for all with P(s, a) > 0.
Proof. Assume (s) S M < oo for some with P(s,a) > O. After To, Nr(3,0) = m S M. The PUCB value
T
aUCB(s,a)= (Grlla)))eG( CpP(s,a) 1 1+m
As T + x0, VT oo, so PUCB(s,a) - co. For a' # with N7 a') + oc:
T PUCB(s,d') Qr(s,d) (r(ssdd)+ CpP(s,,') 1+ Nr(ssa))
Since
T V (+,7() O, PUCB(s, + Qr(s, a') S 1. Then, PUCB(s, > PUCB(s,a')
for large T, forcing selection of contradicting Nr(s,a) Hence, limr Nr(s, co if ass,a)= 0.
Theorem C.4. Given the assumptions of Theorem C.1
unbiased reeural-guides simulations, unique optimal action a", and P(s, a²) > 0, we have (r(s.a) = O(n f for a + a' and Nr(s,, a²) T.
Proof. For a * a", action is selected when:
Qr(s,a)+ CPP((sa)) I+ NT
1+ NT(saa))
As T - + 0o, ()( + (s, a). - ('((s with A. = (*( - Q" 0.
+ Assuming NT(8,a") >> Nr(8, a), the condition approximates to:
T (T AA.
CpP(s,a)
Squaring:
Nr(s, < CpP(s,a)2T
A2
Adapting the convergence analysis of UCB 41] if P(s, a) > 0, (s = O(n If P(, a) = 0. (s, a) may be bounded. Aggregating suboptimal visits:
Σ azaa
NT(s, =((( (Ink))
Accordingly:
= -T- ~ T.
Remark Theorems
C33 and C.4 require P(s,a*) > If P(s, ((s,a") = 0. convergence may fail due
to under-exploration of
26

---

Case 3: Hybrid MCTS (0 < 1)
Theorem C.5. Given the assumptions of Theorem C.1 hnMCTS with the e-UCB policy ensures limr-xc N1 8, a) for a + a', and Nr(8, a²) T.
unbiased simulations, and 0 < c. < 1,
for all E A. Nr(s, O(I
Proof. The expected visit count is:
(8, C = NT c
+ N (s. ().
By Theorem C.1 NyCB
(8, a) for all Hence, E[NT (s, a)] + oc, even if P(s, 0. For (s, a) = 0 7) and NPLCB (8.a) O(hT)),
suboptimal actions, Theorems C.2 and C.4 give NyCB
so:
E = =
for + a'.
For the optimal action:
B[NT(8,a") = Σ B|Nr(s,a). ~
Consequently, (,, + 1, M (8, 0 for + a².
Remark 3. The condition P(s, a²)= ) > is relaxed in hnMCTS due to the vanilla MCTS component
(c. > 0), improving convergence over pure neural-guided MCTS.
C.5.2 Exploration-
Trade-Off
The e-UCB policy, selecting vanilla MCTS with probability & and neural-guided MCTS with probability Ce. balances exploration and exploitation:
Naive MCTS (c.): The UCB exploration term Cu In Eg Nr(s,b)
(r(s,a) promotes exploration when
(s, a) is small, and the Qr(s, a) term promotes exploitation when (r, a) is large.
Neural-Guided MCTS (1 - c.): The PUCB term CpP((s, a) / for low (s, a), moderated by P(s, a). while QT a) drives exploitation.
E. (x,6,) 1+Nr(ss,) encourages exploration
The parameter & controls the trade-off:: High c,: Emphasizes vanilla MCTS, favoring exploration via uniform sampling and random rollouts.
Low c.: Emphasizes neural-guided MCTS, enhancing exploitation via the neural prior P(s, a).
The regret of hnMCTS:
R S c.((KTTTT)) + +1-
where vanilla MCTS contributes O(VKT In T) and neural-guider MCTS contributes O O(VT) (as- suming accurate priors). Tuning c. optimizes the trade-off. with c, > 0 ensuring robust exploration correct misleading priors, enhancing convergence reliability over pure neural-guided MCTS.
C.6 Hyperparameter Setting for GSR
Table 4 describes the hyperparameter settings of GSR for our reported experiments in Sections 4 5) and Appendix [ In these experiments, c, not only promotes better convergence and exploration- exploitation trade-off deseribed in Appendix C.5.I but also facilitates smooth start-up in R-GCN training. Initially, GSR lacks sufficient information about the constraints, making random roll-outs more effective. Therefore. for the first 10% of the training epochs, we set G. exploration. As training progresses and GSR becomes more informed, c, is gradually reduced to 0.1, allowing R-GCN to guide more simulations and prioritize exploitation. Experimental results of the tuning are detailed in Appendix D.4 To promote better complexity control, we also adopt dynamic complexity upper bound strategy. Given the maximum complexity upper bound H*, the actual dynamic complexity upper bound H in Algorithms) and is initialized as the half of H* for the first 20% the total training epochs, then linearly increased to H² by the 60% of the training epoch mark, leading to more cxplorations in lower complexity region.
encourage
27

---

Notation
Value
Explanation
The switching threshold for the e-UCB policy in (3), initialized to be I for the first 10% training epochs, then linearly decreased to 0.1 by the 80% the training epoch mark.
Ge
1 1.1.11
Cy
V2
The exploration coefficient for UCB in (3)
Cp
The exploration coefficient for PUCB in
The normalizer for the relation weight Wa in (
Cij
+
[ denotes the total number of operations in the dictionary
Dynamic complexity upper bound of expressions in Algorithms 1 and 2 initialized to be half of the maximum complexity H for the first 20% training epochs, then linearly increased to H by the 60% of the training epoch mark.
H
0.5/1 + H*
The batch size for hnMCTS in Algorithm dependent on the dataset.
B
n
20
The budget size to store expressions for the transplantation.
Table 4: Hyperparameter settings for GSR
C.7 Evaluation Metrics for hnMCTS
During the training, we evaluate I through the reward R(T) = 1/(11 NMAE), where NMAE the normalized mean absolute error defined as:
1
NMAE
M fx())
(6)
By 72
D Additional Experimental Detail
D.1 Experimental Setting
Black-box and Feynman dataset: The Black-box and the Feynman datasets are widely accepted datasets adopted by the state-of-the-art SR benchmarking, SRBench (27). The Feynman dataset is synthetic dataset consisting of 119 physics-inspired solution equations derived from the Feynman Lectures on Physics, and the Black- box dataset is challenging dataset consisting of real-world observation and simulation data without explicit solution expressions [42]
For both the Black-box and Feynman datasets, each of the equations is run ten times with different random seeds. Each of the benchmarking methods is restricted to million evaluations, time budget of 10 hours, and a maximum complexity of 200 for the expressions. We follow the same dictionary
(available operators and variables) and other default experimental settings in SRBench, where the train-test split ratio is 3 1, and the default operation set includes & (+--, x, exp, log, sin, cos, ]. Here æ; denotes the variables in the dataset. Baselines included in SRBench adopt the same hyperparameter settings reported in (27) GSR adopts a hnMCTS batch size 100 and other hyperparameters described in Appendix C.6 For the added baselines, DGSR
[6] and uDSR [9] utilized the reported hyperparameter setting for the corresponding datasets. TPSR
[10 follows the reported À 0 model configurations. SPL [8] for Nguyen benchmarking utilized the maximum Module Transplantation as 20, episodes between module transplantations as 10,000,
the regularization factor 7) as 0.99, and the exploration rate as v/2
Our proposed GSR method is also evaluated in Nguyen's 's SR benchmark dataset 26 and Nguyen's 's SR benchmark with constants dataset 77, two widely adopted synthetic datasets for evaluating the performance and robustness of various SR algorithms.
Nguyen dataset:
The "Expression" column in Sable the ground-truth expressions used to generate synthetic data. In this set of experiments, adopt results of GP, NGGP and SPL from [8] and follow the same
28

---

Expression settings to implement GSR, which uses dictionary set Pu (x +. for the benchmark Ngugen-I to Ngugen-8, and 1 Do {y} for the benchmark Ngugen-9 to Ngugen-12. 12. It is worth noticing that for the benchmark Nguyen-8, vz can be recovered from exp ェサオ log (z)). For Nguyen- 10, x can be recovered from exp (y log (z)). For Nguyen-7 and Nguyen-10, 10, they can also be recovered from log (v* + x² + 1) and sin (ær + y). The maximum complexity H" for each expression is set 50. The batch size and the maximum number of epochs is capped at 1,000, with a total search space of up to one million expressions for MCTS and GSR. The training and testing datasets are divided equally, with 20 randomly generated data points for each.
sin, cos, ezp, log, x}
to 100 for MCTS simulations,
D.2 Nguyen's Benchmark Result
Table shows that GSR outperforms competing SR approaches in almost all tasks, particularly in all complex tasks, including Nguyen-12 12. This advantage is primarily due to GSR's effective management of polynomial terms, facilitated by the equivalence captured by EG, which compresses the search space more efficiently than the expression tree other representations utilized by alternative methods, highlighting considerable potential in the domain of Symbolic Regression.
Benchmark
GP
NGGP
SPL
GSR
N
1² +
99%
100%
100%
100%
Nguyen-2
x²
+
x + +x² tx
90%
100%
100%
100%
Nguyen-3
5 + x + r3 + +2 +x
34%
100%
100%
100%
Ng iguyen-4
6
x + x² + x3
x² +I
54%
100%
99%
100%
Nguyen-5
sin (x²) cos (x)
12%
80%
95%
100%
Nguyen-6
sin (I)) sin (r
11
100%
100%
100%
Nguyen-7
log (x + 1) log (x? +1)
17%
100%
100%
96%
Nguyen-8
V
100%
100%
100%
100%
gguyen-9
sin (I) sin (y²)
76%
100%
100%
100%
Nguyen-10
2sin (I)) cos (y)
86%
100%
100%
100%
Nguyen-11
x
13%
100%
100%
100%
Nguyen-12
x4 x + 1 y
0%
4%
28%
88%
Nguyen-1
3.39x² + 2. 12x + 1.78r
0%
100%
100%
100%
Nguyen-2
0.48x* 4 + 3.39x + 2.12x + 1.78r
0%
100%
94%
100%
Nguyen-5
sin (z²) cos (I) 0.75
1%
98%
95%
100%
Nguyen-8
VI.23z
56%
100%
100%
100%
Nguyen-9
sin (1.5x) + sin (0.5y²)
0%
90%
96%
100%
Average Solution Rate
38.2%
92.5%
94.5%
99.1%
Table Solution Rate of GSR and other baseline models benchmarked on Nguyen's 's SR benchmark dataset and Nguyen's SR benchmark with constants dataset (marked with upper index c). The solution rate is the ratio of ground-truu -trth equivalent solutions in mathematics to the total of parallel experiments for the same equation.
D.3 Additional Scalability and Efficiency Analysis
EGs are DAGs with node sharing, resulting fewer nodes than equivalent ET's (Fig. which helps reduce memory and computational overhead, Moreover, the hnMCTS design does not rely solely on neural-guided search, where the UCB component ensures efficiency during exploration. Together, these design choices contribute the method's practical efficiency as shown in the ablation study (Section 4.2 where GSR achieves strong performance while requiring the least training time among competitive baselines.
Regarding scalability, GSR performs robustly across datasets with varying input dinensionality. For example, in the black-box benchmark (Fig. 4c)), the average training time remains manageable across increasing dimensions, as shown Table
29

---

Dataset Dimension
3
8
13
18
36
Average Training Time (s)
12676
31483
20781
12783
9172
Table 6: Average training time of GSR in the black-box dataset across different dimensions
Sparse Rate Over Epochs
Sparse Rate Over Epochs
OSA
80
OSR(c decay)
MCTS-EG
=SR(c,=1)
hnMCTS-ET
MCTS
Sparse Rate (1)
60
OSR(c 0)
CSR(c = 0.5)
40
20 Percentage Mark of Training Epochs (96)
50
70
30
Percentage Mark of Training Epochs (%)
(a)
(b)
Figure 6: Sparse rate over epochs in the (a) ablation study and (b) c. tuning. Sparse rate is calculated based on the average ratio between the number of invalid expressions sampled per epoch and the total number of expressions sampled per epoch.
D.4 Additional Ablation Analysis
Figure a(a) illustrates the average reward sparsity across different ablation models in the ablation study in Section 4.2 During the initial 10% epochs, all models exhibit similar trends due to the random rollout induced by the warm-up e-UCB policy in hnMCTS. However, beyond this phase, models incorporating hnMCTS with constraints encoding demonstrate significant reduction in reward sparsity, and the models incorporating EG with invariance encoding expedite the decrease of the sparse rate. This highlights the effectiveness of hnMCTS in capturing constraints encoding in integrating priors of expression-level constraints, along with EG's complementary effects of invariance encoding for faster convergence, ultimately enhancing search efficiency.
Ablation Sparsity Analysis:
Sensitivity Analysis on Cet Figure k(b) illustrates the average reward sparsity across different hyperparameter tuning strategies for c, of GSR, based on the same experimental settings as Fig. b]a). When c, is fixed to be 1, hnMCTS acts like vanilla MCTS, so GSR(c, = 1) in Fig. b is equivalent to MCTS-EG in Fig. a). When is fixed to be O, hnMCTS will fully rely on R-GCN guidance.
MCTS When c. is fixed to be 0.5, hnMCTS will switch between vanilla MCTS and neural-guide
by equal chance. GSR(c. decay) represents the same model as GSR in Fig. 6(a), which adopts decaying c, as described in Appendix[C.o. The results of the sparse rate in Fig. b demonstrate that pure neural-guided MCTS (c. = will result in slow convergence at the beginning due the sparse rewards of SR and lack of training of R-GCN, leading to the biased Q(s, a) estimation. But the vanilla MCTS (c. 1) also has the highest sparse rate after the long training due to the lack of constraint prior encoding to promote promising candidates. Consequently, compared to another hybrid MCTS with =055.,, hnMCTS adopting the decaying c, achieves a better balance between the early exploration and latter exploitation, achieving a faster convergence and higher efficiency.
E
Materials Science Application
E.1
Problem Background
The physical world is composed of numerous ions and electrons, governed by quantum mechanics, often referred as many-body problems. While density functional theory (DFT) can simulate simple materials with few ions and electrons, extending DFT to larger systems is challenging. One of the goals relevant to the physical properties of materials is to find an analytical approximation that relates DFT-calculated potential energy to the interatomic distance (i.e., pairwise distances) between atoms, adhering to the strong electrostatic repulsion at short interatomic distances and weaker interaction at long interatomic distances. Using copper (Cu) crystals as an example, we aim to identify analytical relationships between interatomic distances and total potential energy. This functional form is crucial for materials science as it enables large-scale molecular dynamics simulations to study materials'
30

---

mechanical, thermal, and kinetic properties and explain fundamental physical mechanisms at the atomic level.
The following provides some explanations for domain-specifie terms:
Molecular dynamics (MD) simulates the structure and properties of materials under constant or varying environments (e.g., temperature or mechanical strain). MD simulations, the atomic forces are computed from the derivatives of interatomic potential energy functions with respect to the atomic distance placement, which then moves the atoms based on Newton's second law of motion. The updated atomic positions lead to new interatomic potential energy. The calculations are performed iteratively until a desired number of time steps. Statistics of total potential energy, kinetic energy, atomic forces, and stresses, etc.. provide systematic understanding and quantitative estimate of the physical properties materials. MD simulations are often performed for systems with hundreds millions of atoms, making them computationally expensive. This necessitates machine learning algorithms like those employed in our study. Analytical expressions with low complexity and high accuracy particularly valuable for large-scale MD simulations.
In MD simulations, the system state can be specified using ensembles like NVT (constant number of atoms, volume, and temperature), allowing other properties such as pressure and chemical potential to vary.
Atomistic simulations based on machine learning-derived potential energy functions use atomic species and interatomic forces to determine material properties. This method is efficient but requires accurate fitting from first-principles DFT or other quantum mechanics-basee calculations.
Periodic boundary conditions (PBCs) enable the modeling of infinite or effectively infinite systems by repeating "unit cell" in all directions. This unit cell represents the smallest section of the structure that can be repeated to create the correct crystal structure, allowing large systems to be modeled efficiently without losing structural or symmetrical integrity.
E.2 Pratical Usage of GSR's Solution
Framework Generality: GSR is not limited by system size or atom types. The 32-atom copper dataset is chosen to serve as a benchmark to compare GSR with prior work and demonstrate its advantages. The derived potential energy function can be applied to larger systems.
Dataset Relevance: The 32-atom copper dataset is representative of materials studies, consistent with the average system size (~ 31 atoms/structure) of the widely used MPtraj dataset 43)
Scalability: While DFT methods are limited to simulating small systems (1-1000 atoms) over very short timescales (a few picoseconds), our method enables large-scale simulations (millions of atoms) over extended timescales (microseconds)) which can then help determine several physical properties of materials that are not pratical for quantum chemistry methods such as DFT. The analytical function from GSR can be applied to study real material problems with sizes far beyond 32 atoms for long-time dynamics that are completely inaccessible to DFT or other quantum chemistry methods- which is the key purpose of developing a force field based on accurate analytical potential energy functions from the GSR method. As result, one can simulate the mechanical deformation process of single or polycrystalline copper consisting of millions of copper atoms using large-scale molecular dynamics simulations with the potential energy function developed here. For example, studying the novel stacking of copper in the incubation period of crystallization 44 requires the simulation of millions of copper atoms.
DFT-Level Accuracy: The dataset is derived from DFT calculations, enabling DFT-level accuracy but without explicit electron degree of freedom which is the major purpose of developing the potential energy function (i.e., machine learning force field) from quantum chemistry datasets such as DFT. The high accuracy, the analytical nature, and the proper limiting trend at short bond length make this method and the derived potential energy function especially useful in practice.
E.3 Physies Constraints and Explanations
We offer detailed explanations and examples for some physics constraints defined in Section 5.2
(1) Scalar Output (operation-level constraint): Because the input variable r is vector of variable length for different samples, Σ operator must be included in the expression before any
31

---

(2) Unit Consistency (operation-level constraint): It is defined to ensure that the unit calculation aligns with physical meaning. as const can introduce an extra unit as a coefficient. For specific example of the unit consistency constraint, the expression of (r² +r) is not physically meaningful, as r has a length unit Å. and we would have an inconsistent unit (A² +Å) from the expression. But a constant с can include unit so that we have a consistent calculation with ,² + C т. The introduced const can ensure that the final output unit aligns with the target E's unit as electron-volt (eV; eV 1.6 10-19 Joule) from the input r's unit Angstrom (A; 1A= 10-10m).
(3) Electrostatic Repulsion (expression-level constraint): It's subtle consequence of Pauli's ex- clusion principle of quantum mechanics: When the distance between two atoms approaches zero, their electron wavefunctions start overlap significantly, which is excluded by Pauli's principle. As a result, the atomic orbitals hybridize and form molecular orbitals with bonding and antibonding characters, and electron density is hence pushed away from nuclei, leaving the repulsive nucleus-nuclens electrostatic interaction being the dominant one and approaching infinite potential energy at very short distances.
E,4 Application Experimental Setting
Dataset Description: The 32-atom copper dataset consists of 150 snapshots generated by the Vienna Ab initio Simulation Package (VASP), firstpriiciples DFT package. Each sample includes total formation energy E' as the target and a crystal structure of copper as the variable. To map a copper crystal structure E, we first convert the 3D coordinates of the structure into pairwise distances between atoms within the unit cell in the crystal under the periodic boundary condition (PBC) [45]. then use the surrogate model as the interatomic potential energy function f(-) to obtain E = (r). During the conversion, we only consider pairwise distances within a cutoff range of Teutoff 5Å.
SR Settings: In the symbolic regression configurations, we utilize a dictionary set + {+,-,x, A, E.er., const, r), where "E" is used for the summation-based aggregation opera- tion, "const" denotes constants optimized through non-linear regression. The maximum complexity H is set at 35, and the batch size B is 1 1,000, with 5,000 episodes allocated for both MCTS and GSR methods. For the black-box CGCNN method, we adhere to the default hyperparameter settings with maximum of 5,000 epochs. For the train-test split, follow [16] with 50% for training and 50% for validation.
E,5 Application Result Figure
Figure compares the data fitting plots by different SR solution models in Table the performance discussions in Section
2 corresponding to
(0lime/armmm)
GPL MAE: 3.53 3. 3.68mmWWatom)
VCTS MAE: 3.65 3.9L)meVatoom
GSR MAE 2.65 2 2:410eeisboom)
ttadktted Eteergy (mmea//etmm)
400-
Troe cood Einengg Inn-tlatomt
Pretcted mestrrrm
800
220
330
Syec Energy
220
300
::00
1.00
300
150
100
200
300
400
100
300
400
120
300 Energy imeViatom!
200
400
10C
200
300
450
Energy (meVianom)
Total Energy (meViattom)
3meVharom)
(a)
(b)
(c)
(d)
Figure 7: Training (orange dots) and Testing (blue dots) MAES of the formation energy predictions by (a) CGCNN, (b) GPI, (c) MCTS, and (d) GSR on DFT-based molecular dynamics simulations of FCC copper are provided in the captions of the plots on the top, including both testing MAE (left)
and training MAE (right). The dashed lines mark the identity mapping.
E.6 Physical Meaning
Figure pressnnts the potential energy curves as a function of the interatomic distance (bond length) for the expressions generated by the GP1, GP2, MCTS, and GSR in Table 2 These curves demonstrate
32

---

whether the fitted models satisfy electrostatic repulsion at zero bond length. Figure 5 shows that GP1 and GP2 violate the constraint due the finite value at zero bond length (r=0), In contrast, GSR and MCTS yield infinite potential energy at r=0, satisfying the principle. Such a constraint is crucial not only for the underlying physies but also for practical applications in materials system simulations. The missing divergence at 0 may cause atoms to collapse to each other during molecular dynamics simulations and yield incorrect results.
Additionally, Fig. 5 reveals two critical issues with GP1 and GP2: (i) a discontinuity r=1A, and (ii) an infinite value when approaches 5Å. These issues can lead to incorrect energy and force predictions, forcing atoms to remain unrealistically close together.
E.7 Generation of the New Dataset for Testing Model Transferability
To test the transferability of the models, we generated test dataset from DFT calculations using the VASP package [3]] where the projector augmented wave method was applied to treat core electrons [46]. We used the Berdew-Burke--rnzzene
(PBE) form of exchange-correlation energy k- k-point sampling grid of 4 X 4 X 4 [48]. A hydrostatic compression applied to FCC copper with strain of -0.2 along all lattice vectors. We then performed first- peisnprinssples molecular dynamics simulations in the NVT ensemble for total 6,000 steps with time step of 3 fs at two different temperatures, i.e. 300 K and 1,400 K. For the NVT calculation at each temperature, we excluded the first 1,000 steps of the initial equilibration process, and extracted the atomic structures for every 100 steps from the remaining equilibrated 5,000 steps, which yields total 100 snapshots (or samples) to test model transferability.
functional within the generalized gradient approximation [47 and a Monkhorst-
The compressed dataset has the following physical meaning:
Understanding matter in extreme conditions such as high pressure and high temperature is an important and active subject of materials research. Copper is one of the systems of particular interest. As shown in [49], the density in experiments reaches ~~ g/em", double the density the standard condition of 8.95 g/em², that is, the volume contraction by 50%. Another example is done by [50) at the National Ignition Facility (NIF) at the U.S. Lawrence Livermore National Laboratory (LLNLL,, in which the solid copper was even compressed 28 g/em² terapascal conditions, corresponding to ~677 volume contraction.
Scientific Context:
Practical Necessity: Another key motivation to apply large compression is to provide a more ac- curate trend away from equilibrium towards - 0. This is particularly important as the machine learning interatomic potentials or machine learning force fields very often have wrong limiting behavior. When they were applied to simulating long-time dynamics at high temperature high pressure, there will be an increasing probability of "direct crossing or fusion" atoms which are pure artifacts due to the wrong limiting trend, consequently, the results can be completely non non-physical and wrong.
33


Reviews:
[{"review_id": "EqJcsjxVOS", "reviewer": "Reviewer_RJ7g", "summary": "Based on the invariance-aware design utilizing graph representation, this paper addresses the long-standing challenges of redundant representations (due to the inability to capture equivalences) and reward sparsity in symbolic regression. The proposed framework is validated across both ​​synthetic benchmarks​​ and ​​real-world scientific datasets​​, demonstrating superior performance in expression recovery and constraint satisfaction. Furthermore, it provides ​​convergence guarantees​​ for the hybrid neural-guided Monte Carlo Tree Search (hnMCTS) algorithm, as detailed in the appendix.", "questions": "Is the expression represented by a graph unique? If not, it seems that the limitation has not been completely solved. Can the relative improvement of uniqueness be quantified? Are there any future work directions that can be pointed out for the redundant representation problem that graph representation still cannot completely solve?", "limitations": "1. Theoretical Limitation: Unverified Uniqueness of Expression Graphs​​\nWhile the Expression Graph (EG) representation resolves redundancy for basic operator equivalences (e.g., commutative operations in Fig. 2b), ​​this paper does not establish EG's capacity to guarantee canonical forms for expressions involving complex nested operators​​. \n\n​​2. Empirical Limitation: Narrow Validation of Physical Constraints​​\nThe physics-informed validation ​​exclusively relies on a single-domain case study​​ (interatomic potential energy in materials science, Sec. 5).", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength\n1. This paper solves the ​​redundant representation problem​​ in symbolic regression by introducing Expression Graphs (EGs) that inherently capture mathematical equivalences.\n2. This paper solves the ​​sparse reward challenge​​ by embedding expression-level constraints (e.g., physical laws) directly into the search process via R-GCN-guided MCTS.\n3. This paper addresses the reproducibility gap through open-sourced code and detailed configurations (Appendix D.1/F.3).\n\nWeakness\n1. Uniqueness gap​​: Lacks proof for EG's representation uniqueness in complex cases (e.g., nested operators).\n2. Only material science case validates physical constraints—needs cross-domain tests (e.g., fluid dynamics).", "was_revised": false, "final_justification": "", "rebuttal": ""}, {"review_id": "ncXjwGhqgU", "reviewer": "Reviewer_YFGA", "summary": "The authors present Graph-based Symbolic Regression (GSR) to address two limitations in symbolic regression: permutation invariance in expression trees (ET), and sparse rewards due to inadequate constraint incorporation during sampling. GSR uses Expression Graphs processed by Relational Graph Convolutional Networks (R-GCN) and employs hybrid neural-guided Monte Carlo Tree Search (hnMCTS) to integrate constraints directly into the search process. The method is evaluated on standard benchmarks and a materials science application for predicting copper interatomic potentials.", "questions": "1. For the Feynman and Strogatz dataset, could you provide numeric accuracy (R2) scores, expression complexity, and training time comparisons (similar to the black-box dataset analysis)? This would help assess whether the benefits extend beyond recovery rate. Also, if the authors can add some other strong and recent baselines to the comparison, the paper will be improved. I suggest baselines such as PySR (Cranmer), uDSR (Landajuela et al.), TPSR (Shojaee et al.), and MDLFormer (Yu et al.).\n\n2. How does your approach to expression-level constraints differ from post-hoc evaluation in tree-based methods? For constraints like f(r→0) = ∞, it appears complete expressions must still be generated and evaluated - could you clarify how this constitutes \"direct constraint incorporation\"?\n\n3. The method looks to be computationally expensive. What additional complexities graph representation adds to the method? Could you provide wall-clock time comparisons across different dimensions (scalability)?", "limitations": "The authors mention some of the limitations throughout the paper; this is a minor suggestion but they could be more explicit in bringing the limitations in the final discussion/conclusion of the paper.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n* The paper approaches real limitations in symbolic regression. The expression graph representation appears to have the potential to address permutation invariance issues that affect traditional expression trees, and the R-GCN architecture seems well-suited to leverage edge features for encoding operator relationships. The hybrid hnMCTS approach provides a reasonable way to balance exploration and exploitation while incorporating learned guidance.\n* The experimental evaluation is comprehensive, covering multiple standard benchmarks (Feynman, Black-box, Nguyen) + an application problem to molecular dynamics inter-atomic potentials --  However, the paper could benefit from some additional comparisons and discussions (see questions)\n* The paper is extensive and well-written and the technical execution appears sound based on the algorithmic details provided in the appendix. The combination of graph neural networks with MCTS represents a reasonable engineering approach that shows consistent improvements across benchmarks.\n\nWeaknesses:\n* The novelty claims may be overstated. Mainly, the claimed advantages in constraint handling are unclear. Both operation-level and expression-level constraints appear implementable in tree-based methods using the same mechanisms (action masking and post-hoc evaluation). The distinction between traditional \"post-hoc penalties\" and the proposed \"direct constraint incorporation\" seems unclear since expression-level constraints inherently require complete expressions for evaluation.\n* The practical application scope seems limited. Training on 32 copper atoms to predict properties of the same system doesn't demonstrate the generalizability typically expected from scientific methods. While the physics constraints are interesting, they appear to be relatively basic requirements (unit consistency, electrostatic repulsion) rather than complex physical principles.", "was_revised": true, "final_justification": "I still remain slightly positive overall about this work and keep my score.", "rebuttal": ""}, {"review_id": "WrbWTadTdi", "reviewer": "Reviewer_bn25", "summary": "The paper claims two improvements to symbolic regression: first a new DAG representation which removes the well-known permutation problem of expression trees. Second, use of a graph neural network to propose Monte Carlo moves. The GNN proposals are used as part of a policy which can also draw conventional random proposals. The GNN takes only encodings of the expression graph, and does not look at the data.  Acceptances depend on the reward, ie the fit of the numerically-optimised tree to the data.   A RL algorithm receives and propagates rewards to states, where rewards are the fit from model to data.     Numerical parameters are fit with standard BFGS for each expression (which also handles constraints), and a form of chunking type memory is used in \"trans\" to store small useful tree fragments which then get their own special syntax to insert into new expressions.   The method beats state of the art methods including AI-Feynman, Operon and DGSR on a mixture of two standard datasets for scientific discovery, Feynman and blackbox, measured by solution rate;  and beats them all except AIFeynman (which appears to give a NaN score) measured by R^2;  and beats a smaller collection of models on a new dataset.", "questions": "The conversion from expression tree to expression graph looks similar to moves made by optimizing compilers, could you comment on any links here?  Is there a standard compiler theory for how to do this optimally?  How different is the proposed method from such standards?\n\nDoes the data play any role in generating the expression graph proposals, or are they proposed based only on the current expression graph?  If the latter, why not also encode and look at the data when making proposals?\n\nIs the learning done from scratch each time a new data set is presented for fitting, or it is a foundation style model where such proposals are learned for the general case then applied to each data set?\n\nWhy does AIFeynman have a NaN score for R^2?   Could this be tweaked in any way to obtain an actual number (e.g. removing a few problematic items?) to make it more comparable with the new method?\n\nWhy is the new physics dataset only tested against 4 less well known models in Appendix Table 2, rather than the full set of state of the art models used in the Feynman+blackbox experiment?   (maybe there is something about this dataset than means only models with some certain property can be run on it?)\n\nBeating all of the state of the art methods on the standard data sets seems impressive.   But it is sometimes possible to do this by developing new algorithms against the test data sets throughout their lives.  Even if numerical parameters are not specifically optimised to fit the test set, different choices in the structure of the algorithm itself can be tried out and used to evolve the system to fit the test set. And this is a very complex algorithm which has many such decisions to be made.   Can the authors confirm that this was not done during algorithm development, and that all development was done using different data sets from the test sets presented here?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 1, "significance": 3, "originality": 2, "strengths_and_weaknesses": "The main paper does not contain enough detail to reproduce the work, though this is provided in a very long appendix which is effectively the real paper.   (I'm not sure what is NeuroIPS policy on this, many papers now seem to be basically ignoring page limits via this approach, so look more like journal papers but split between the two documents. AFAIK appendices are supposed to be for optional, nice-to-have extensions rather than for the main algorithmic content?)  \n\nLikewise the lit review in the main paper is very short, but the appendix does contain the full one, where it is good to see the history of SR traced back to Koza's GP in the 1990s (though not Simon's original BACON from the 1970s, which should maybe get a mention?), through to recent relevant DNN-based systems like AI-Feynman.\n\nConvergence proofs are provided in the appendix.\n\nI'm struggling to understand the part of the algorithm for making the neural network expression proposals.   I can see the expression graphs being encoded with graph NNs.  I can see discussion of an MDP and rewards (i.e. the expression, when numerically optimised, being a good fit to the data, appendix eqn 6), and Q-values. There is GNN message passing in adx eqn 5.   From my reading, idea idea appears to be to encode the expression trees via a GNN, then use this encoding as state, in a reinforcement learning algorithm which selections actions that modify the graph (via encoding representations).   I'm struggling to understand why the proposal method doesn't look at the data -- only the current state of the graph -- and what gets learned by doing so?    I can see the  learning makes use of the data via the rewards, so it looks like the network is learning heuristics which say \"if the current graph has these kinds of features, then try modifying it in these ways\".", "was_revised": false, "final_justification": "", "rebuttal": ""}, {"review_id": "El7IvERT8A", "reviewer": "Reviewer_ytHy", "summary": "This work tackles the challenging problem of symbolic regression by introducing novel methods to learn more effective representations of mathematical equations and improve reinforcement learning (RL) signal design.\n\nThe key contributions of the paper are twofold:\n\n1. **Expression Graphs (EG) Representation and R-GCN-Based Encoding**: The authors propose a graph-based representation of equations—referred to as Expression Graphs (EG)—which captures structural and semantic relationships within mathematical expressions. To encode these graphs into meaningful representations, they employ R-GCNs, enabling the model to better understand and generalize from complex equation structures.\n\n2. **hnMCTS for Reward Densification in RL**: The second contribution introduces hnMCTS. This component is used to densify sparse reward signals in the RL framework, allowing the agent to evaluate and learn from partial equation expansions rather than relying solely on complete solutions.\n\nThe paper includes a series of experimental evaluations and ablation studies that support the effectiveness of the proposed approaches.", "questions": "Please refer to weakness 3,4, and 5.", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 1, "significance": 4, "originality": 2, "strengths_and_weaknesses": "**Strengths:**\n\n1. **[Significance]** The paper is well-motivated and tackles two rather significant issues in symbolic regression.  \n2. **[Novelty of EG]** The proposal of using EG to represent equations is novel.  \n3. **[Quality of experiments]** The author conducted experiments and relevant ablation studies to verify the various components of the algorithm.  \n\n**Weaknesses:**\n\n1. **[Clarity]** The mathematical formulation is somewhat unclear. For example, $\\mathbf{P}$ was initially defined as a \"prior probability matrix\" on line 130, but was redefined as a predicted output of R-GCN on line 132. The dimensionality of $\\mathbf{P}$ is also unclear. As another example, the actions $a_t$ are not defined in the text. It is thus unclear how the authors parameterized $a_t$. There are many other examples of ambiguous notations in the manuscript. The unclear formulation poses challenges to thoroughly understanding the method.  \n2. **[Novelty of hnMCTS]** The proposed hnMCTS seems to be a simple combination of two well-known algorithms (MCTS and neural-guided MCTS). While the motivation and effectiveness of hnMCTS are clear, the method is lacking in novelty.  \n3. **[Limitation of EG]** The proposed EG method can only account for a limited set of invariance relations in equations (i.e., commutative equivalences). In practice, there exist other modalities that lead to two equations being equivalent (e.g., distributive laws and other common algebraic equalities). It would be beneficial to analyze the amount of invariance EG can resolve versus the amount of residual invariance.  \n4. **[Limitation of sequential generation]** In the MDP formulation, the EG is generated *sequentially*. For the same EG, the generations will be different if the nodes are traversed in different orders. In other words, although EG captures the invariance in equations, the sequential traversals on EG might sacrifice efficiency again (thereby unnecessarily complicating the search space). The authors should clarify and disucss this aspect. \n5. **[Incorporating equation-level constraints]** The author mentioned that the new method allows incorporation of equation-level constraints without seeing a completed equation (line 92). However, the exact way to incorporate such constraints is not clearly stated.", "was_revised": true, "final_justification": "The rebuttal partially addresses my concerns. Point 1, 4, and 5 in Weakness needs further elaboration.", "rebuttal": ""}]

Author rebuttal:

Meta-review: This paper proposes a graph neural network–based symbolic regression approach that addresses two key challenges in the current symbolic regression literature: **(1) redundant representations** and **(2) sparse rewards**. The method introduces permutation invariance into the graph representation of expressions to eliminate redundancy, and it encodes constraints during the expression sampling process, thereby avoiding unnecessary rollouts of full expressions.

**Strengths:**

* Novel approach that tackles fundamental issues in symbolic regression.
* Clear methodological contributions: permutation-invariant representations and constrained expression sampling.
* Potential to significantly improve efficiency and robustness of symbolic regression.

**Concerns:**

* Some reviewers found the algorithmic process not fully clear.
* The lack of proof of uniqueness was noted, though this was addressed in the rebuttal.
* Evaluation is somewhat limited, focusing primarily on material science applications.

**Assessment:**

The reviewers agree that the work is novel and makes progress on fundamental problems in symbolic regression. The concerns raised were satisfactorily addressed in the rebuttal, which also included useful clarifications.

**Recommendation:**

I recommend **acceptance**. I strongly encourage the authors to incorporate the clarifications and additional results presented in the rebuttal into the final version, and to follow through on their commitments to strengthen the paper. Doing so will ensure the work has maximum clarity and impact.
Decision: Accept (poster)