GenSR: Symbolic regression based on equation generative space
Authors: Qian Li
Yuxiao Hu
Juncheng Liu
Yuntian Chen
Venue: ICLR 2026 Poster (2026)
Keywords: Symbolic Regression; Equation Generative Latent Space
Primary area: neurosymbolic & hybrid AI systems (physics-informed, logic & formal reasoning, etc.)

Abstract: Symbolic Regression (SR) tries to reveal the hidden equations behind observed data. 
However, most methods search within a discrete equation space, where the structural modifications of equations rarely align with their numerical behavior, leaving fitting error feedback too noisy to guide exploration.
To address this challenge, we propose GenSR, a generative latent space–based SR framework following the "map construction $\rightarrow$ coarse localization $\rightarrow$ fine search" paradigm. Specifically, GenSR first pretrains a dual-branch Conditional Variational Autoencoder (CVAE) to reparameterize symbolic equations into a generative latent space with symbolic continuity and local numerical smoothness. This space can be regarded as a well-structured "map" of the equation space, providing directional signals for search. At inference, the CVAE coarsely localizes the input data to promising regions in the latent space. Then, a modified CMA-ES refines the candidate region, leveraging smooth latent gradients. 
From a Bayesian perspective, GenSR reframes SR task as maximizing the conditional distribution $p({\rm Equ.}|{\rm Num.})$, with CVAE training achieving this objective through the Evidence Lower Bound (ELBO). This new perspective provides a theoretical guarantee for the effectiveness of GenSR. Extensive experiments show that GenSR jointly optimizes predictive accuracy, expression simplicity, and computational efficiency, while remaining robust under noise.

Published as a conference paper ICLR 2026
GENSR: SYMBOLIC REGRESSION BASED ON EQUATION GENERATIVE SPACE
Qian Li Yuxiao Hu Juncheng Liu', Yuntian Chen
Shhangaai Jiao Tong University, Shanghai, China
"asterrn Institute of Technology. Ningbo, China
Thhe Hong Kong Polytechnic University, Hong Kong, China
*Imperial College London, London, England
qian1012050 jtu. edu. cn, yuxias ao hu hwwconnect.
junchengl i 1u23e wetee3emmmmprieeeeee imperiiplaaie ac. uk, wwwen@eitech.eeeeeeee
ABSTRACT
Symbolic Regression (SR) tries to reveal the hidden equations behind observed data. How- ever, most methods search within a discrete equation space, where the structural modifi- cations of equations rarely align with their numerical behavior, leaving fitting error feed- back too noisy guide exploration. To address this challenge, we propose GenSR" generative latent space-based SR framework following the "map construction - coarse localization -> fine search" paradigm. Specifically, GenSR first pretrains a dual-branch Conditional Variational Autoencoder (CVAE) to reparameterize symbolic equations into generative latent space with symbolic continuity and local numerical smoothness. This space can be regarded as well-structured "map" of the equation space, providing di- rectional signals for search. At inference, the CVAE coarsely localizes the input data to promising regions in the latent space. Then, modified CMA-ES refines the candi- date region, leveraging smooth latent gradients. From a Bayesian perspective, GenSR reframes SR task maximizing the conditional distribution p(Equ. Num.)), with CVAE training achieving this objective through the Evidence Lower Bound (ELBO). This new perspective provides theoretical guarantee for the effectiveness of GenSR. Extensive ex- periments show that GenSR jointly optimizes predictive accuracy. expression simplicity, and computational efficiency, while remaining robust under noise.
INTRODUCTION
Symbolic Regression (SR) a core machine learning task aimed at discovering interpretable mathematical expressions that explain observed data. Unlike traditional regression, which assumes a fixed functional form (e.g., linear, polynomial), SR explores the equation space to identify functions that best fit the data. Formally, given N observation samples {x,, 3/1 ) the objective is to find mathematical expression f + y that captures the underlying data distribution. Owing to its flexibility and interpretability, SR widely used in scientific discovery and engineering (Xu et al. 2023 Chen et al. 2022), where understanding the underlying equations is as important as predictive accuracy. SR effectively bridges the gap between black- box models and interpretable scientific insight. Consequently, growing research efforts focus on developing SR algorithms that balance fitting accuracy, computational efficiency, and expression simplicity.
"Corresponding author
https: //gpthub...ttttt// wwwp://otaaka22///..R22666GG

---

Published as a conference paper at ICLR 2026
Many prior SR methods search for equations directly in a discrete symbolic space, using heuristic strategies such as genetic programming (GP) (Stephens et al. 2015 Virgolin et al. 2021). ). Monte Carlo Tree Search (MCTS) (Sun al.). or other combinatorial optimization methods (Yu al. 2025| Zhang et al.] 2025a). These methods iteratively adjust the symbolic structures of equations through various operations such as crossover, mutation, and tree expansion to search for better solutions. Unfortunately, such modifications do not reliably reduce fitting errors and may even lead the search into suboptimal regions. This is because edit distance and numerical behavior similarity are unrelated measures: equations with similar structures can exhibit drastically different numerical behaviors, and vice versa. Consequently, in the discrete equation space where similarity is measured by edit distance, the numerical error feedback neither corresponds to edit distance nor provides reliable search direction. The search relies heavily on random mutations, combina- tions, and backtracking, resulting in an unstable search trajectory and high time complexity. Therefore, the discrete equation space with edit distance inefficient and even unsuitable for SR.
In this paper, we reparameterize the discrete equation space into a continuous latent space, which can be regarded as a "map" of the equation world (where equations are represented as vectors and dimensions naturally define search directions). To enable efficient search in this space, two key points must be addressed: (1) The space should be generative. A generative latent space, constructed by generative models, supports continuous sampling such that latent vectors decode into syntactically valid equations. This property enables smooth interpolation and fine-gre graiigr search in the equation latent space. In contrast, pretraining- based methods (Meidani et al., Kamienny et al. 2022) typically learn discriminative embeddings for similarity or prediction, rather than modeling a generative distribution. Such discriminative spaces are fragmented and filled with non-decodable regions, hindering continuous search and producing invalid samples.
(2) The space should be well-organized. Although some methods attempt to construct a generative latent space for equations (Mežnar et al. 2023 Popoy et al. 2023). the Euclidean distances within it reflect only edit distances and are unrelated to the equations' numerical hehavior. Consequently, numerical errors cannot be mapped into meaningful Euclidean distances for efficient optimization, and such latent spaces still cannot support smooth, efficient search. Ideally, equations with similar numerical behaviors should be close in the latent space, allowing smooth numerical error signals to guide exploration. Yet expressions like cos² x, 1- -in? 2 x, and their Taylor expansions exhibit the same numerical behavior but entirely different symbolic forms. Placing such expressions in close proximity would force the decoder to map nearby vectors to drastically different expressions, destabilizing decoding. Thus, a key design challenge is balancing symbolic continuity and numerical smoothness. Inspired by human localization first identifying a coarse region before precise targeting- we construct generative latent space with global symbolic continuity and local numerical smoothness. The former ensures stable decoding and induces clear structural clustering for coarse localization, while the latter allows fitting errors to provide directional signals for fine-grained search once the relevant region is located.
Building this insight, we propose GenSR, novel framework that follows the "map construction - coarse localization + fine search" paradigm. GenSR first pretrains a dual-branch Conditional Variational Autoencoder (CVAE): The posterior branch learns distributions over symbolic structures given symbolic and numerical inputs, whereas the prior branch learns distributions over numerical features from numer- ical inputs only. Joint optimization of reconstruction loss and KL divergence yields a Generative latent space with global symbolic continuity and local numerical smoothness for SR, serving as a "map" of the equation world. At inference, input data are mapped to an initial localization distribution by prior branch, which highlights promising latent regions. A modified CMA-ES algorithm then contracts this distribution toward the optimum, exploiting smooth directional signals in the latent space for efficient convergence. From a probabiliitic view, GenSR reformulates SR as a Bayesian optimization problem, i.e., maximizing (Equ. [Num.), with CVAE training corresponding to ELBO optimization. This perspective provides the theoretical guarantee for GenSR. Extensive experiments show that GenSR consistently achieves optimal balance of predictive accuracy, expression simplicity, and runtime efficiency, while maintaining robustness under noisy conditions.
2

---

Published as a conference paper ICLR 2026
Our contributions come in four parts: (1) GenSR constructs continuous generative latent space unifying symbolic continuity and local numerical smoothness, addressing the inefficiency of discrete equation spaces; (2) GenSR follows a new search paradigm of "map construction + coarse localization -> fine search" to achieve effective SR; (3) GenSR provides a new Bayesian perspective for SR and uses ELBO optimization solve it; (4) Extensive experiments provide a thorough analysis of latent space properties and demonstrate that our method achieves a better balance among accuracy, equation complexity, and efficiency.
2 RELATED WORK
The field of SR has evolved through variety of approaches. The most common SR methods rely on combinatorial optimization over discrete equation space using heuristic search. Genetic Programming
(GP) techniques (Stephens et al. 2015 Virgolin et al. 2021 Burlacu et al. 2020, Cranmer 2023) evolve
populations of expressions through selection and mutation, while Monte Carlo Tree Search (MCTS) meth- ods (Sun et al.) systematically explore the expression tree space. More recent neural-guided approaches have integrated reinforcement learning (Petersen et a..). Transformer-based planning (Shojaee et al., 2023). retrieval- -aggeeeedd generation (Zhang et al.] 2025a). and minimum description length principles (Yu et al., 2025) to enhance search efficiency. Large language models have also been applied to generate equations through prompting techniques (Shojaee et al.; Grayeli et al.] 2024); however, due to inherent hallucination issues, LL-M-based methods still resemble enumeration in the discrete space and require many queries reach high-quality solutions. Ultimately, symbolic regression over discrete equation spaces remains funda- mentally inefficient, as structural similarity poorly reflects numerical behavior, leaving search unguided and error feedback uninformative.
Alternative strategies have sought avoid discrete search entirely. Ecquence-to-sequenee models (Kami- enny et al. 2022 Li al. [2222] Biggio et al., (022) generate equations directly but do not enforce alignment between symbolic forms and numerical behavior during training. limiting their generalization. The SNIP framework (Meidani et al.) made progress by learning a shared latent space through contrastive learning, creating representations useful for similarity assessment. However, because SNIP learns discriminative rather than generative space, many points in its latent space do not correspond to valid equations, making it inadequate for the generative search demands of SR.
Our work addresses these limitations by learning continuous, generative latent space using dual-branch CVAE. This approach transforms the discrete equation space into a structured manifold where nearby points correspond to expressions with similar symbolic forms and numerical properties. The resulting space pro- vides meaningful gradients that enable an efficient three-stage search strategy- map construction, coarse localization, and fine-gra amaiigr search- addressing the core weaknesses of previous methods.
3 GENSR: SYMBOLIC REGRESSION BASED ON GenERATIVE LATENT SPACE
As illustrated in Fig. [ the core idea of GenSR is to construct continuous generative latent space for equations (which can be regarded perform efficient search directly within this space. The dual-branch CVAE consists of a Transformer-based encoder-decoder network. two branch-specific networks, and a feature fusion module. The CVAE is pre- trained on approximately 5 million synthetic forms. At inference, GenSR leverages the prior branch to produce latent space and refines the solution in the latent space using a degenerate version of CMA-ES (Hansen 2016), enabling fast and precise search for the distribution of candidate equations.
"map" of the equation world) using dual-branch CVAE and to
pairs to ensure coverage of diverse functional initial localization distribution in the
3

---

Published as conference paper ICLR 2026
Training: Generative Space Construction
Reconstruction Loss
MGG of)
f(x) Prefix Order
Posterior Network
Position
Transformer Encoder
KL Loss
Feature Fusion
Transformer Decoder
Prior Network
N(poof)
Inference: Search the Generative Space
Schematic of the evolution
Prior Network
Feature Fusion
(d)
Tranformer Encoder
Decoder
Transformer
Refinement
Constants
f(x)
Ym
NGn
Tokenize &t Embbed
Evaluation: Fitting (roo.
Complexity
Update
Degraded CMA-ES
Frazen
Sampling
Figure I: The overview of GenSR. During training, the dashed lines denote the prior branch, while the solid lines indicate the posterior branch. During inference, only the prior branch is used.
3.1 PREPARATIONS
To construct the generative latent space, we first prepare synthetic data and preprocess input numerical samples and equations, laying the foundation for dual-branch CVAE pre-training.
Preparation of synthetic data. To construct the generative latent space, we pretrain conditional VAE (Fig. [T on large synthetie dataset. Following Kamienny et al. (2022); Meidani et al., construct dataset in which each equation f' is paired with numerical examples (æ, y) (= where 1) denotes the dimensionality the input variables. Note that we control the expression length to push complex equations
toward the periphery of the latent space, mitigating overfitting to complex solutions during search.
Pre-processing of numeric samples. Following Kamienny et al. (2022), we tokenize each numerical sample using base- 10 floating point representation, rounding to four significant digits and decomposing into sign, mantissa (0-9999), and exponent (E- 100-E100), e.g., 0.7895 is tokenized as [+, 7895, E-4]. The m samples for each equation f are represented as {(x, }) E I=*x3 Icarnable embedder with dimensionality reduction maps these token embeddings to E Rmxd which serves as the numerical input to the Transformer encoder (Kamienny al. 2022).
Pre-processing of equation. Each expression f is represented as binary tree linearized in prefix order (Lample & Charton). Operators, variables, and integers mapped to dedicated tokens, while con- stants are tokenized in the same way as numeric samples. Special tokens (BOS)) and [(EOS)] are inserted mark the start and end of each sequence. Particularly, sequences padded fixed length match the sequence length numeric samples. Each token is embedded into a dddimennsional space with positional encoding, yielding F Roxd as the equation input the Transformer encoder.
4

---

Published as a conference paper ICLR 2026
3.2 PRE-TRAIN THE GENERATIVE LATENT SPACE FOR EQUATIONS
We adopt the tions, which can be regarded as a "map" of the equation space. This map provides initial localization of high-probability regions and enables smooth directional search toward optimal solutions, achieving efficient equation discovery. The branch CVAE consists of a posterior branch and a prior branch. Follow- ing Kingma& & Welling 2013); Sohn et al. (2015). we assume the conditional distribution of latent variables is a multivariate independent Gaussian distribution. Next, we describe each branch and its training details.
CVAE shown in Fig. T to learn well-organized generative latent space for equa-
Posterior Branch (solid lines in Fig. [ Comprises 8-layer Transformer encoder, posterior network, a feature fusion MIP,, and an 8-layer Transformer decoder. takes both X (numerical samples) and F (symbolic equation) as inputs, learning the posterior distribution g( X, F) [N o network. Using the reparameterization trick, we sample n latent vectors Z = { pi + σi O e) where E (0, I) and @ denotes element-wise multiplication. The feature fusion MLP aggregates salient information across sampled latent vectors, and the Transformer decoder reconstructs the original symbolic equation, ensuring the latent space captures the continuity of symbolic structures.
If via the posterior
Prior Branch (dashed lines in Fig. I The prior branch shares all modules (Transformer encoder/decoder, feature fusion MIP) with the posterior branch, except that the posterior network is replaced by a prior net- work. It takes only numerical samples X as input and learns the prior distribution p( X) = ( T) which captures the numerical behavior of equations.
Training Process: During training, both branches are jointly optimized on paired data (F, X) using a combined objective: (1) Reconstruction loss Lise enforces continuity of symbolic structures by accurately decoding equations from the latent space: (2) KL divergence DxL aligns the posterior distribution g(= X F) with the prior distribution p(z X), ensuring local numerical smoothness. This alignment of distributions regularizes the prior network to support reliable inference when the true equation is unknown. We apply KL annealing during training to mitigate posterior collapse.
3.3 ANALYSIS OF GENERATIVE LATENT SPACE PROPERTIES
We further explain why GenSR in Sec. 3.2 yields generative latent space that is globally continuous in symbolic structure and locally smooth in numerical features.
Continuity of symbolic structure. Let (z) = ( o (1)21) and (z) = NN o (2)2 denote the two output posterior distributions for input sample pairs (F(1), x()) and (2)). respec- tively. We define the high-confidence region for each distribution at confidence level
(3)
R.d
E Rd
Σ
xa
i= 1,2,
12
where x².0 is the a-quantile of the chi-squared distribution with degrees of freedom. For two distributions, their n-sampled sets g = (2,"0 n²0) + ( @e) where ~ NO, (N.O))))) J. = 1, 2, typically fall within the respective high-probability regions R.!") and R!!2) Here, unlike the traditional VAE framework that sam- ples only single 2, we employ repeated probabilistic sampling and a feature fusion module to ensure that the reconstruction loss influences the high-probability regions of the distribution, thereby strengthening the
symbolic continuity of the latent space. If there exists z" E R!!") n !?? it is shaped by the reconstruction objectives of both f()) and f(2), yielding a decoded structure that interpolates between them. This encour- ages smooth symbolic transition in the latent space. Trained on diverse equations, such local interpolation
generalizes globally, clustering structurally similar equations in adjacent regions and inducing a latent space with continuous symbolic structure.
5

---

Published as a conference paper ICLR 2026
Smoothness of numerical feature. We leverage the prior branch capture characteristics of numerical distributions and align them with the symbolic latent space through KL divergence loss, imposing numer- ical constraints on distribution. This design resembles contrastive learning in SNIP (Meidani al.), but with key distinction: we align entire distributions rather than individual points, which ensures local numer- ical smoothness within the generative latent space. Specifically, the KL loss (x(]x( , F) || p(=[X)) regularizes the posterior distribution to stay close to the prior distribution. For the independent Gaussian distributions assumed in our CVAE (E1 = of Σ₂ = =oiI), this divergence simplifies to:
d
(113)-22 - (2.j) )
of
o. 1 In
DxL ( 1l
Σ j=1
2
2
2.j
2
Minimizing this term in our bual-branch CVAE constrains the distance between the posterior and prior means while aligning their variances, ensuring that the high-probability regions-ratter than isolated points-and are matched to preserve local numerical smoothness.
Consequently, GenSR establishes well-structured generative latent space for SR. We will provide more validation of the properties of the space in Sec. 5.3
3.4 SEARCH IN THE GENERATIVE LATENT SPACE
During inference, since the ground-true equation is unavailable, only the prior branch used. The inference process is illustrated in Fig.D
For given set of numerical samples {(x, }) we first preprocess them as described in Sec. 3.1 The prior branch then encodes them into an initial Gaussian distribution N ((Ho. I) )) in the latent space. This distribution provides coarse localization, highlighting regions likely to contain equations that fit the data well. For example, high-probabillyy regions align with plausible symbolic families (e.g., trigonometric, exponential) and input dimensionality, offering a structured starting point for fine-grained search.
Jim
Since the prior and posterior latent spaces are only approximate, directly decoding the vector μ4 in an end-to-end manner may not yield the most accurate result. To refine the initial distribution A ((Ho, o toward the optimal solution, use modified Covariance Matrix Adaptation Evolution Strategy (CMA- ES) (Hansen 2016). leveraging the numerical smoothness of the latent space for efficient convergence. The algorithm follows:
1. Initialization: Set the initial search distribution N N(uo I) (with no and do the output of the prior network).
2. Iteration (for each generation i):
(a) Sample multiple latent vectors { from (i, of
(b) Decode {z,} into candidate equations (z) using the well-trained decoder.
(c) Refine the constants of {z) via BFGS optimization.
(d) Evaluate the fitness of candidate equations: Fitness = R - w complexity.
(e) Select the top-p candidates and update the distribution parameters + p441 and σ; + 0+
using the CMA-ES update rules. This yields the updated distribution
The prior branch provides well-organized generative latent space for CMA-ES. thereby enabling faster and more stable convergence. However, CMA-ES becomes costly as the latent dimension grows. To improve efficiency, we introduce two key modifications to standard CMA-ES:
1. Diagonal Covariance Assumption: Inspired by Ros & Hansen (008)), we restrict the covariance matrix Σ to diagonal form の² updating only the variance of each latent dimension indepen- dently. This aligns with the latent variable assumption in VAEs (Kingma & Welling (013) (mod- eled as independent Gaussian components) and reduces the number of parameters to optimize.
6

---

Published as a conference paper ICLR 2026
2. Top-k Variance Update: setting others to zero). This focuses search on the most uncertain and relevant directions, acceler- ating convergence without compromising solution quality.
Update only the top latent dimensions with the largest variances (while
These modifications enable rapid and precise contraction of the initial distribution, facilitating rapid local- ization of the target equation. Fig. | provides a schematic of the CMA-ES distribution contraction process.
4 GENSR: SYMBOLIC REGRESSION FROM BAYESIAN PERSPECTIVE
GenSR can be naturally interpreted as probabilistic inference, reformulating SR task as a Bayesian problem: given numerical samples X mR the goal is to infer posterior distribution p( [ over equation space that best explains the data, i.e., maximize p( F ). Unlike traditional SR tasks seeking a single optimum, this Bayesian formulation reasons over the entire posterior distribution of candidate equations, enabling principled uncertainty quantification and probabilistic comparison across candidates.
To approximate the intractable distribution p( where Rd is latent vector. Following Sohn et al.
we introduce variational distribution q(> X, F).
22015) derive the Evidence Lower Bound (ELBO) of the marginal log-likk ligeelildoood log [ X) as follows. The full derivation is provided in Appendix
log logr(FIX = ((2 X, F) log p( F F
([F,(X) dz q(z X log qx F)
gg(= x X.F) [log p( F X =]]
DKL (q(z|X, p(z]x))
(1)
To optimize this objective, GenSR employs posterior encoder ( F), prior encoder pa.e (2 x) and decoder po (F) X, 2) to parameterize q(||X, F), p(=[X) and PF| 2). respectively. Then, training our dual-branch CVAE with Lyee and DKL is equivalent to maximizing the ELBO maximizing [=)()X_F)) [log p( (F) X, =)] corresponds to max- imizing the probability of reconstructing F using z sampled from the output of a (x) i.e., minimiz- ing the reconstruction loss Lree; and the second term in inequality
Maximizing p( X) amounts to maximizing the ELBO
I
Inequality[Ioooreppponn corresponds to the KL-divergence loss Sec. 3.2 which aligns the variational distribution (= X, F)) with the prior distribution (=|
During inference, we sample from the prior distribution Pa (= to reconstruct the equation F, which approximates the intractable expectation ( [log pe [F) F| =]] with [=,(= [log Pe
F X =)]. To
mitigate the inherent approximation error, we refine the prior distribution explained in Sec. 3.4
Notably, GenSR differs from Holt et al., which approximates y(Token Num.) at the token level without constructing latent distribution over symbolic equations, and also from Mežnar et al. (2023); Popov al. (2023), which only models (.Eq..). GenSR explicitly formulates the SR task as maximizing the conditional distribution p(Equ./ Num. from Bayesian inference perspective. To the best of our knowledge, GenSR is the first framework that realize SR based on estimating and optimizing the ELBO of p(Equ. NNum.). Under GenSR framework, our dual-branch CVAE pretraining model corresponds directly to ELBO, while the CMA-ES refinement process serves as an approximate optimization of the conditional variational distribu- tion. This not only provides GenSR with solid theoretical guarantee, but also introduces new Bayesian perspective and technical route for SR.
5 EXPERIMENT
5.1 EXPERIMENTAL SETTINGS
To validate our method, we conduct evaluations on the SRBench benchmark (La Cava et al. 2021 Cavalab 2022). which consists of 119 Feynman equations, 14 ODE-Strogatz challenges, and 57 black-box regres-
7

---

Published as conference paper ICLR 2026
sion tasks without known underlying functions. We compare against 18 baseline algorithms, covering both pretraining-base based methods and diverse heuristic-based methods. The evaluation metrics include accuracy R², time complexity, and equation complexity. To better balance and visualize the trade-offs among these three metrics, we employ the Pareto front for representation. Comprehensive information on datasets, met- rics, and baseline methods can be found Appendix G.1
G.2 anddG.3 For clarity, this section presents only the key experimental results; the full results are deferred to Appendix G.5 Additional analyses, including
the ablation study and other task based on our generative space, are provided in Appendix G.4
5.2 COMPARISON OF GENSR WITH OTHER METHODS
FG4T
Aipatiin Compleatt
Atr
Dum
00
62000
Adan Ronk
Hea Zan
Figure Pareto front results on the Feynman dataset. The z-axis shows the mean test R² rank, while the y- axis shows equation complexity rank (left) and time complexity rank (right). Solid lines indicate the optimal Pareto front, and dashed lines show lower-ranked fronts from bottom-left to top-right.
AP7
PBA7
APP-P
SHP
SPP-GP
Opman
GP-GOMEA
14:00
RAG-5R
Ourz
0.80
0.55
0.30
0.35
1.00
10'
16
10'
20
40
100
g' Score
Tine Compleattvis)
Equation Complexity
Figure 3: Comparison on the Strogatz dataset under different noise levels. Subplots (left to right) report R² score, time complexity (s), and equation complexity. Noise levels represented by blue circles (0.000), orange squares (0.001), green triangles (0.01), and red diamonds (0.1), with error bars indicating standard deviations. Only methods whose mean R² across noise settings exceeds 0.9 included.
Pareto performance analysis. As shown in Fig. 2 our method achieves balanced optimization across accuracy (R²),, equation complexity, and time complexity, and consistently lies on the rank-1 Pareto front, outperforming competing baselines overall. The superiority stems from the generative continuous latent space learned via CVAE training, which preserves global structural continuity while ensuring local numer- ical smoothness. Built upon this representation, the GenSR efficiently enables hierarchical optimization
8

---

Published as conference paper ICLR 2026
paradigm of "map construction - coarse localization + fine search." In contrast to discrete symbolic spaces that rely on heuristic mutations and backtracking, our well-organized generative latent space provides sta- ble directional signals, making the search both more efficient and less prone to overfitting. Consequently, GenSR achieves robust improvements in accuracy, expression brevity, and computational efficiency.
Performance under noisy targets. As shown in Fig. 3 our method achieves the highest R² across varying noise levels, with smaller variance than competing baselines, indicating stable and robust fitting under perturbations. Moreover, it requires significantly less runtime and yields equations of lower com- plexity. These results highlight that our generative latent space preserves structural organization and local smoothness even in the presence of noise, thereby providing search with stable optimization directions and avoiding the random oscillations and inefficient exploration typical of discrete symbolic spaces. Benefiting from this design, GenSR maintains favorable balance of accuracy, efficiency, and compactness under noisy conditions, demonstrating strong robustness.
5.3 VISUALIZATION COMPARISON OF LATENT SPACE
We compare GenSR with E2ESR (Kamienny et al. 2022) and SNIP (Meidani et al.). E2ESR employs Transformer to predict target equations, while SNIP leverages contrastive pretraining to align symbolic forms with numeric data. Unlike GenSR, Etheir spaces are all discriminative rather than generative.
Structural separability comparison of latent spaces. As shown in Fig. 4 different methods induce distinct distributional characteristics of equation latent vectors. The latent space of E2ESR fails to effectively separate function types and input dimensions: although can partially distinguish two input variants of the log operator, the trig and exp components remain highly entangled, resulting in mixed representation. SNIP achieves slightly better separation at the operator level. yet lacks intra-class compactness; for instance, the exp-5D) cluster splits into two distinct regions. accompanied by noticeable overlaps across categories. contrast, our GenSR exhibits clearly structural disentangled latent space, where both different function types and input dimensionality are well separated. This space provides GenSR with favorable conditions for initial equation localization while ensuring decoding stability.
8 op 3D
0 40-40
Sp3D
(a) EZESR
SNIP
(c) GenSR
Figure 4: 2D t-SNE visualization of latent variables from E2ESR, SNIP, and GenSR. The legend distin-
guishes six categories, corresponding to equations from three representative function families, each evalu- ated under 2D and 5D input dimensionality, illustrating the clustering behavior of the leamnd latent spaces.
9

---

Published as conference paper ICLR 2026
1.00
Avvengee Normllzzd
0.50
0.00
(4) exp SD
(b) exp 2D
(4) trig
(4) trig 2D
(e) log 5D
(1) log 2D
Figure 5: 2D t-SNE visualization of GenSR latent variables for equations from three function families (exponential, trigonometric,, logarithmic) under 2D and 5D input settings, shown in subplots (a)-(f). Colors indicate the average of normalized values, as displayed in the accompanying color bar.
Numerical continuity within function Classes. As for each observed sample {æ;, ji normalized mean of (ys } to represent the numerical feature. Fig. visualizes the local numerical feature distribution of GenSR, where each subfigure shows the numerical feature distribution of one cluster in Fig, (c). For each structural cluster, GenSR exhibits smooth distribution of numerical features, which enables the search algorithm refine convergence along smooth numerical directions after localizing high probability areas. Details and additional latent space experiments are provided in the Appendix
the
CONCLUSION AND FUTURE WORK
We introduced GenSR, a generative latent space-based SR framework that maps the discrete symbolic equa- tion space into continuous, well-structured space. GenSR first jointly enforces symbolic continuity and local numerical smoothness through dual-branch CVAE, and then enables efficient symbolic regression through coarse localization and fine-grained search. This design bridges combinatorial symbolic search with continuous optimization, yielding balanced trade-off among accuracy, complexity, and efficiency. More- over, GenSR organizes distinct equation families into separate latent regions, naturally supporting tasks such as equation classification, which we further demonstrate in the appendix. These results highlight the broader potential generative latent spaces foundation for equation discovery and, more generally, for combina- torial scientific search. Future works include incorporating richer generative priors and physical constraints, and leveraging more advanced generative models to construct more powerful and generalizable equation spaces, thereby advancing interpretable and domain-aware machine learning for scientific discovery.
7 ACKNOWLEDGEMENT
This work was supported by the National Natural Science Foundation of China (12572266), and the High Performance Computing Centers at Eastern Institute of Technology, Ningho, and Ningbo Institute of Digital Twin.
10

---

Published as conference paper ICLR 2026
REPRODUCIBILITY STATEMENT
For the reproducibility of our results, we have provided a detailed description of our methods and experi- mental setups in Sec. 3.1 3.2 and Appendix F We also confirmed the robustness of our results through the
experiment (Appendix G.5 In addition, to further facilitate the reproduction, will release our codes and the checkpoints for the trained models.
REFERENCES
Ignacio Arnaldo, Krzysztof Krawiec, and Una-May O'Reilly. Multiple regression genetic programming. In Proceedings of the 2014 Annual Conference on Genetic and Evolutionary Computation, pp. 879-886, 2014.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neu- ral symbolic regression that scales. In International Conference on Machine Learning, pp. 936-945. Pmir, 2021.
Bogdan Burlaeu, Gabriel Kronberger, and Michael Kommenda. Operon c++ an efficient genetic pro- gramming framework for symbolic regression. In Proceedings of the 2020 genetic and evolutionary computation conference companion, pp. 1562-1570, 2020.
Cavalab. Srbench: A living benchmark framework for symbolic regression, 2022. URL https: Accessed:
77ccvva
2024-09-26.
Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree boosting system. In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pp. 785-794, 2016.
Yuntian Chen, Yingtao Luo, Qiang Liu, Hao Xu, and Dongxiao Zhang. Symbolic genetic algorithm for discovering open-form partial differential equations (sga-pde). Physical Review Research, 4(2):023174, 2022.
Miles Cranmer. Interpretable machine learning for science with pysr and symbolicregression. jl. arXiv preprint stXiv:2305.01582 2023.
Fabricio Olivetti de Franca and Guilherme Seidyo Imai Aldeia. Interaction Interaction-transfrmmaaa transformation evolutionary algorithm for symbolic regression. Evolutionary computation, 29(3):367-390, 2021.
Qihao Duan, Bingding Huang, Zhenqiao Song, Irina Lehmann, Lei Gu, Roland Eils, and Benjamin Wild. Janusdna: A powerful bi-directional hybrid dna foundation model. arXiv preprint 1rXiv:2505.17257,
2025.
Chen Feng, Ziquan Liu, Zhuo Zhi, Ilija Bogunovic, Carsten Gerner-Beuerke and Miguel Rodrigues. Prosac: Provably safe certification for machine learning models under adversarial attacks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, pp. 2933-299 2941, 2025.
Chen Feng, Minghe Shen, Ananth Balashankar, Carsten Gerner-Beuerle, and Miguel RD Rodrigues. Noisy but valid: Robust statistical evaluation of LLMMs with imperfect judges. In The Fourteenth International Conference on Learning Representations, 2026. URL https: http::::///ppereviiiee net forummid-
hEhxreaLdd
Yoav Freund and Robert Schapire. decision-theoretic generalization of on-line learning and an applica- tion boosting. Journal of computer and system sciences, 55(1):119-139, 1997.
11

---

Published as a conference paper ICLR 2026
Arya Grayeli, Atharva Schgal, Omar Costilla Reyes, Miles Cranmer, and Swarat Chaudhuri. Symbolic regression with learned concept library. Advances in Neural Information Processing Systems, 44709, 2024.
Systsmsm, 37:44678
Nikolaus Hansen. The ema evolution strategy: A tutorial. arXiv preprint arXiv:1604.00072, 2016.
Samuel Holt, Zhaozhi Qian, and Mihaela van der Schaar. Deep generative symbolic regression. In The Eleventh International Conference on Learning Representations.
Yuxiao Hu, Qian Li, Xiaodan Shi, Jinyue Yan, and Yuntian Chen. Domain knowledge-cnhaneed multi- spatial multi-temporal pm2. 5 forecasting with integrated monitoring and reanalysis data. Environment international, 192:108997, 2024.
Yuxiao Hu, Qian Li, Dongxiao Zhang, Jinyue Yan, and Yuntian Chen. Context-alignment Activating and enhancing Ilm capabilities in time series. arXiv preprint #rXiv:2501.33747, 2025.
Tao Huang, Rui Wang, Xiaofei Liu, Yi Qin, Li Duan, and Liping Jing. Detecting misbehaviors of large vision-language models by evidential uncertainty quantification. arXiv preprint arXiv:2602.05535 2026.
Yanhao Jia, Ji Xie, S Jivaganesh, Hao Li, Xu Wu, and Mengmi Zhang. Seeing sound, hearing sight: Un- covering modality bias and conflict of ai models in sound localization. arXiv preprint 1221, 2025.
Ying Jin, Weilin Fu, Jian Kang, Jiadong Guo, and Jian Guo. Bayesian symbolic regression. AAAI, 2020.
Alexa-Ale Kamienny, Stéphane d' Ascoli, Guillaume Lample, and François Charton. End-to-end symbolic regression with transformers. Advances in Neural Information Processing Systems, 35:10269- 10281, 2022.
Zhaolu Kang, Junhao Gong, Jiaxu Yan, Wanke Xia, Yian Wang, Ziwen Wang, Huaxuan Ding, Zhuo Cheng, Wenhao Cao, Zhiyuan Feng, et al. Hssbench: Benchmarking humanities and social sciences ability for multimodal large language models. arXiv preprint mrXiv:2506.03922 2025.
Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. Lightghm: A highly efficient gradient boosting decision tree. Advances in neural information processing systems, 30, 2017.
Diederik Kingma and Max Welling, Auto-encoding variational bayes. arXiv preprint arXiv:1312.611 2013.
William La Cava, Lee Spector, and Kourosh Danai. Epsilon-lexicase selection for regression. In Proceedings of the Genetic and Evolutionary Computation Conference 2016, pp. 741-748, 2016.
William La Cava, Tilak Raj Singh, James Taggart, Srinivas Suri, and Jason Moore. Learning concise representations for regression by evolving networks of trees. ICLR, 2019.
William La Cava, Bogdan Burlacu, Marco Virgolin, Michael Kommenda, Patryk Orzechowski, Fabrício Olivetti de França, Ying Jin, and Jason H Moore. Contemporary symbolic regression meth- ods and their relative performance. Advances in neural information processing systems, 2DBI( 2021.
Guillaume Lample and François Charton. Deep learning for symbolic mathematics. International
Conference Learning Representations.
12

---

Published as a conference paper ICLR 2026
Qian Li, Yuxiao Hu, and Yang Cao. Adaptive perona-malik model based on dynamical threshold for image multi-noise removal with details preservation. Computers & Mathematics with Applications, 137:28 -4,, 2023a.
Qian Li, Yuxiao Hu, Ye Liu, Dongxiao Zhang, Xin Jin, and Yuntian Chen. Discrete point-wise attack not enough: Generalized manifold adversarial attack for face recognition. In Proceedings of the IEEE/CVE conference on computer vision and pattern recognition, pp. 20575-20584, 2023b.
Qian Li, Yuxiao Hu, Yinpeng Dong, Dongxiao Zhang, and Yuntian Chen. Focus hiders: Exploring hidden threats for enhancing adversarial training. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 24442-24451, 2024.
Wenqiang Li, Weijun Li, Linjun Sun, Min Wu, Lina Yu, Jingyi Liu, Yanjie Li, and Songsong Tian. Transformer-basee model for symbolic regression via joint supervised learning. In The Eleventh International Conference on Learning Representations, 2022.
Kazem Meidani, Parshin Shojace, Chandan K Reddy, and Amir Barati Farimani. Snip: Bridging mathemat- ical symbolic and numeric realms with unified pre-training. In The Twelfth International Conference on Learning Representations.
Sehastian Mežnar, Sašo Džeroski, and Ljupčo Todorovski. Efficient generator of mathematical expressions for symbolic regression. Machine Learning, 1 2023.
Randal S Olson, William La Cava, Patryk Orzechowski, Ryan J Urbanowicz, and Jason H Moore. Pmlb: 2017.
large benchmark suite for machine learning evaluation and comparison. BioData mining,
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 riek-ng policy gradients. In International Conference on Learning Representations.
Sergei Popoy, Mikhail Lazarev, Vladislav Belavin, Denis Derkach, and Andrey Ustyuzhanin. Symbolic expression generation via variational auto-encoder. PeerJ Computer Science, 9:e 1241, 2023.
Steven J Rigatti. Random forest. Journal of Insurance Medicine, 47(1):31-39, 2017.
Raymond Ros and Nikolaus Hansen. A simple modification in ema-es achieving linear time and space complexity. In International conference on parallel problem solving from nature, pp. 296-305. Springer, 2008.
Michael D Schmidt and Hod Lipson. Age-fitness pareto optimization. In Proceedings of the 12th annual conference on Genetic and evolutionary computation, pp. 543 544,, 2010.
Parshin Shojaee, Kazem Meidani, Shashank Gupta, Amir Barati Farimani, and Chandan K Reddy. Llm- sr: Scientific equation discovery via programming with large language models. In The Thirteenth International Conference on Learning Representations.
Parshin Shojace, Kazem Meidani, Amir Barati Farimani, and Chandan Reddy. Transformer-based planning for symbolic regression. Advances Neural Information Processing Systems, 36:45907-45919, 2023.
Kihyuk Sohn, Honglak Lee, and Xinchen Yan. Learning structured output representation using deep condi- tional generative models. Advances neural information processing systems, 28, 2015.
Trevor Stephens et al. gplearn: Genetic programming in python, with a scikit-learn inspired api. Journal of Engineering for Gas Turbines and Power, 30, 2015.
13

---

Published as conference paper ICLR 2026
Steven H Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. CRC press, 2018.
Fangzheng Sun, Yang Liu, Jun Wang, and Hao Sun. Symbolic physics learner: Discovering gov- crning equations via monte carlo tree search. The Eleventh International Conference on Learning Representations.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Neural Information Processing Systems, 33:4860-4871, 2020.
Marco Virgolin, Tanja Alderliesten. and Peter AN Bosman. Linear scaling with and within semantic backpropagation--ban baasd genetic programming for symbolic regression. In Proceedings of the genetic and evolutionary computation conference, pp. 1084-1092, 2019.
Marco Virgolin, Tanja Alderliesten, Cees Witteveen, and Peter AN Bosman. Improving model-based genetic programming for symbolic regression of small expressions. Evolutionary computation, 29(2):211-2)
2021.
Hao Xu, Jinglong Lin, Dongxiao Zhang, and Fanyang Mo. Retention time prediction for chromatographic enantioseparation by quantile geometry-enhanced graph neural network. Nature Communications, 14(1): 3095, 2023.
Yilong Xu, Yang Liu, and Hao Sun. Rsrm: Reinforcement symbolic regression machine. ICLR, 2024.
Xiaokang Yang, Yuntian Chen, et al. Riesz neural operator for solving partial differential equations. In The Fourteenth International Conference on Learning Representations.
Zihan Yu, Jingtao Ding, Yong Li, and Depeng Jin. Symbolic regression via mdlformer-guided search: from minimizing prediction error to minimizing description length. In The Thirteenth International Conference on Learning Representations. 2025.
Hengzhe Zhang, Qi Chen, Wolfgang Banzhaf, Mengjie Zhang, et al. Rag-sr: Retrieval-augmented generation for neural symbolic regression. The Thirteenth International Conference on Learning Representations, 2025a.
Hengzhe Zhang, Qi Chen, Bing Xue, Wolfgang Banzhaf, and Mengjie Zhang. Lim-meta-se In-context learning for evolving selection operators in symbolic regression. arXiv preprint arXiv:2505.18602 2025b.
14

---

Published as a conference paper ICLR 2026
A THE USE OF LARGE LANGUAGE MODELS (LLMS)
Large Language Models (LLMs) were used solely as auxiliary tools for language polishing and minor edit- ing. They did contribute the research ideation, experimental design, or writing to an extent that would qualify them as contributors. All substantive research ideas, methodologies, and results presented in this paper are entirely the work of the authors.
B DERIVATION OF p(F)
FF X) MMXIMIZATION
We refer to Sohn et al. 2015)) to approximate the intractable distribution X) using a variational dis- tribution g( X, F) that can be any distribution (which is assumed multivariate independent Gaussian distribution in our work). The detailed derivation is as follows.
log logp(F (
g(2 X, F)) Fog F xX)d
p(F,z|X dz
a( X, F) log
p(= p(=|x, X, F)
F(F, x\) g( F)
a(> x(=|X, x, F) log
dz
(2||X, F) p(z F, x)
dz / [(2]|x F) log
[]=||x F)
a(2 Xx Fog log
dz
((2|X, F)
p(= F, x)
((F,z|X) ((2|X,F)
g(=1 X, Fog
log
dz DeL (,(=|XX F, X))
Since DKL (1(2 X, ())|p(=/FFF)))))))))))))))) X)) 2 0, thus, we have the lower bound log (
log > J. ((=]X,F)
( x
((=|X, F) dz
p(F|z, X x)x(((X)
((1|x, (EXX,)) FFllyy
dz
d(2|X, F)
((*|X,F))) log (FFlg((FFzz) X X)dz + q(211x F) log X()
dz
(x=x|XFF))
q(= F) log p 2, x) DRL ()(=\x,
(x(= blog F X, 2)] DEL (1(2|X, F)]lp(z(X)).
To optimize this lower bound, GenSR employ posterior encoder go (2 X, F), prior encoder (a, and decoder Pe(F [F x) to parameterize q( F), p(z and [ z). respectively. Here, o, y. and z denote the model parameters respectively, and 0 comprises the parameters of the Transformer decoder and the feature fusion module. Therefore, training the CVAE with reconstruction loss and KL divergence is equivalent to maximizing the lower bound, i.e., maximizing the conditional distribution [F] X).
the Transformer encoder, the posterior network, and the prior network,
15

---

Published as a conference paper ICLR 2026
C T-SNE DATASET CONSTRUCTION AND SUPPLEMENTARY EXPERIMENTS
C.I T-SNE DATASET CONSTRUCTION
Fig. F and Fig. 5 we constructed a dedicated dataset to analyze the separability of function categories and the continuity of numerical features in the latent space. The dataset covers three representative classes of functions: exponential functions, trigonometric functions, and logarithmic functions. The exponential class includes the operator exp(-), the trigonometric class includes sin(.), cos(-), and tan(-), while the logarithmic class uses log(-).
To systematically examine the impact of input dimensionality latent space distribution, we generated two cases for each function class: 2-dimensional input (D = 2) and 5-dimensional input (D = 5). This setting
allows us verify whether the latent space maintains clear structural separability and numerical continuity under both low- and higher-dimensional inputs.
The expression generation process is follows: first, we specify the dimensionality of the variable x, and then combine the hasic operators +, x with the corresponding function operators from each cate- gory. To ensure comparability, we constrain the expression length such that the number of tokens (including coefficients, operators, and variables) is less than 25.
In terms of scale, for each function class and dimensionality combination, generated 5000 equations, and for each equation, 200 data points were sampled to construct the corresponding dataset. This design provides sufficient statistical coverage and ensures that the t-SNE visualization results are representative.
In summary, this dataset systematically spans function categories, input dimensionalities, expression struc- tures, and sample sizes, serving solid basis for validating the global structural separability and local numerical feature continuity of the latent space.
C.2 INTERPOLATION RESULTS IN LATENT SPACE
Table I: Linear interpolation between two expressions (1+2)) in latent space across methods. For each pair, we decode from =rratio) = (1 ratio)=1 + ratio 22 with ratio {0,, 0.25, 0.5, 0.75, 1}. Expressions marked in red are syntactically invalid.
ratio
Ours
E2ESR
SNIP
Pair 1: Expression 1: cexp(x) sin(x + Expression 2: log(z + cos(x:)
0
cexp(r)sin(r cexp(x) sin(æ c) csin(#) cexp(c- c) sin ) + log(c x) cos(c) exp(c) + log(æ + c) cos(x)
cexp(x) sin(x c) Invalid Invalid
cexp(x) sin(x - c) exp(# s)s((()) csin(x) exp(x) (z) + ccos(x + c) Invalid
0.25
0.5
0.75
log(x)) sin(x) cos(z)
1
log(æ + cos(2)
log(z + c) cos(2:
log(z + e) cos(x)
Pair 2: Expression I: cos(cx) + Exprcssion 2: sin(x - c)
0
x cos(cr) +c ccos(cz) + c cos(2) ccos(c- u) + c- cos(x c0)(c.2) sin(x c)
cos(c + Invalid
x cos(ca)) + cos(c) + csin(x) Invalid Invalid sin(x-
0.25
0.5
Invalid
0.75
Invalid
sin(r
x sin(x-
We compare the continuous interpolation capabilities of generative (our GenSR) and discriminative (such SNIP and E2E) latent spaces through interpolation experiments in the latent space. For each method,
16

---

Published as a conference paper ICLR 2026
first feed the numerical samples of Expression into its numeric encoding branch. In GenSR, this branch corresponds to the prior branch that only uses numeric inputs, while in E2ESR and SNIP it is implemented as their respective numeric encoders. The output of this branch taken as the latent representation 21- Expression 2 is handled in the same manner to obtain 22.
We then construct intermediate latent points by linear interpolation.
z(ralio) = (1 - -ratio) 2 + ralio 22
with ratio E {0, 0.25, 0.5, 0.75, 13. Each interpolated latent vector z(ratio) is subsequently passed through the corresponding decoder to generate the intermediate expression. This unified procedure enables fair comparison of different methods in terms of structural smoothness, and decoding robustness under latent- space interpolation, and allows to assess whether their latent representations support continuous interpo- lation between equations.
Table 1 reports the linear interpolation between two expressions (1~2) in the latent space. Our method consistently produces syntactically valid and meaningful expressions at each interpolation step, thanks to the design of continuous generative latent space that preserves the decodability of equations under inter- polation. In contrast, E2ESR and SNIP frequently yield syntactically invalid expressions (marked in red), indicating that their discriminative latent spaces lack sufficient syntactic regularization and often degenerate into structures that cannot be converted into valid prefix forms. Moreover, the transitions generated by our method are smooth: expressions gradually transform from 1 to demonstrating the structural continuity of the latent space. By comparison, E2ESR and SNIP exhibit drastic structural changes even under small variations in the latent variables, suggesting less stable latent representations. For clarity, all constants are uniformly denoted by to eliminate the influence of specific numeric values.
D RECONSTRUCTION PERFORMANCE OF POSTERIOR VS. PRIOR BRANCH
Table 2: Reconstruction error on synthetic expression datasets. We report Levenshtein (edit) distance (mean + std; lower is better) between an expression and its reconstruction. Two operator sets are considered: (i) LogExp = (+, X, +, log, exp}: (ii) Trig = LogExp U {sin, cos, tan). Each dataset has 2k samples; the maximum expression length is shown in the name.
Dataset (ops, length, vars)
Posterior branch
Prior branch
LogExp-15 (vars < 3)
0.094
(#0.023)
0.106 (#0.031)
LogExp-20 (vars < 5)
0.095 (+0.027)
0.112 (+0.029)
LogExp-25 (vars < 5)
0.107 (+0.018)
0.120
#0.002
Trig-15 15 (vars < 3)
0.115 (+0.032)
0.137
(+.0019)
Trig-20 (vars S 5)
0.113 (+:0.026)
0.142
(+.0.43)
Trig-25 (vars < 5)
0.146
(+-0.041)
0.171
(+0.038)
Synthetic dataset construction. To systematically evaluate the reconstruction capability of our varia- tional autoencoder, we construct six synthetic datasets under two operator sets: (i) LogExp, which includes basic arithmetic operators {+, -, x, } along with log and exp, and (ii) Trig, which extends LogExp by adding trigonometric functions {sin, cos, tan). For each operator set, we generate three datasets with max- imum expression lengths of 15, 20, and 25 tokens, constraining the number of variables to 3, 5, and 5, respectively. This range covers expression lengths slightly longer than those commonly seen in real-world benchmarks (the Feynman dataset has average length of 13.56), allowing performance under increasing structural complexity. Each dataset contains 2k randomly synthesized sam- ples.
evaluate reconstruction
17

---

Published as conference paper ICLR 2026
We quantify reconstruction fidelity using the Levenshtein (edit) distance, which counts the minimum number of insertions, deletions, and substitutions required to transform the reconstructed expression into the ground- truth one. This metric directly measures symbolic accuracy rather than numeric approximation quality, making it a natural choice for evaluating exact expression reconstruction.
reports the reconstruction performance using latent codes from both the poste- rior and prior branches. We summarize three consistent observations:
Results analysis.
Table 2
Posterior branch consistently outperforms the prior branch. This is expected, as the posterior branch is trained by directly maximizing the expectation (x [log (F X, x)] using the
ground-truth latent variables supervision, resulting lower edit distances and more accurate reconstructions.
Prior branch performs slightly worse but remains competitive. Unlike the posterior branch, the prior branch is not trained to maximize the above expectation directly. Instead, it optimized by aligning its output distribution with the posterior through the KL divergence loss, which leads to slightly higher edit distances. Nevertheless, the reconstruction quality remains reasomably good. This is crucial, since during inference can only sample from the prior distribution. The results indicate that the KL loss effectively aligns the prior with the posterior, enabling accurate recon- struction even without access ground-truth latent variables.
Expression complexity has a controlled impact reconstruction. As expression length in- creases and trigonometric operators are introduced, the reconstruction error grows gradually but remains well-behaved. The posterior branch maintains near-exact reconstruction even for length- 25 expressions, validating that the learned latent space preserves sufficient structural information.
Overall, these results demonstrate that the learned latent space is both expressive and well-regularized: it enables faithful equation reconstruction from posterior samples, while the prior remains sufficiently aligned support reliable generation at inference time.
ABLATION STUDIES
To better understand the design choices in GenSR, conduct ablation studies both the latent space configuration and the CMA-ES search procedure. These experiments reveal how different hyperparameters affect accuracy, equation complexity, and runtime, and guide choice of robust default settings.
18

---

Published as conference paper ICLR 2026
E.1
LATENT DIMENSION OF GENERATIVE SPACE
Noise-Free
Noise: le-2
Noise: le-1
Equatiio Complexity 3
3
Equation Complexity
Equariio Complexity
0.7
%
0.8
X
2
-2
0.0
0.2
61
R2 Eipattion Complestty
R² Squatio Coupieity
R² Gpuation Complemny
0.2
128 116
112
768
136
254
312
266
0.2
216 Latent dimension
128
312
266
Latent dimension
Latent dimension
Figure 6: Effect of latent dimension GenSR. The X-axis shows the latent dimension of the pretrained model (64, 128, 256, 512, 768), while the left and right y-axes report R² and equation complexity. Results are reported under three noise levels: Noise-Free, x 10 10-2, and 10-!
Fig. 6 presents the results of pretraining GenSR with different latent dimensions. The results in Fig. 6 To isolate the evaluation the latent space itself, we did not apply CMA-ES optimization and instead directly compared the outputs of the pre-trained model. In the Noise-Free setting (left subplot), R² improves steadily as the latent dimension increases, reaching its maximum at 512 before showing a slight decline at 768. Equation complexity grows moderately for smaller dimensions but exhibits marked surge beyond 512. indicating that excessively large latent spaces may encourage overly complex expressions. Under noisy conditions, the 768-dimensional model suffers from complexity, with performance degradation exacerbated as noise levels increase. Based on these observations, we adopt dimennnisli latent space for pretraining, as it offers the best trade-off and ensures robustness under noise.
pronounced drop in accuracy and further inflated
E.2 CMA-ES HYPERPARAMETERS
2200
210
lexitty
2.555
034
R²
Complecctyi
128
256
322
170
Figure 7: Effect on CMA-ES performance.
(1) Effect of updating principal dimensions. In the CMA algorithm, adopt diagonal approximation by retaining only the variance terms of the covariance matrix, and further update only the top-k principal dimensions with the largest variances to reduce computational cost. Fig. 7 illustrates the effect of different I values. We observe that R² improves as increases and reaches a near-saturation point at i = 256, with almost further gains beyond this dimension. In contrast, the search time grows approximately linearly
19

---

Published as conference paper ICLR 2026
with k, leading a significant increase computational cost for larger k. Balancing accuracy and efficiency, we therefore fix k 256 in all subsequent experiments, as it provides the best trade-off.
(2) Other hyperparameters. We further study three intrinsic hyperparameters of CMA-ES: population size initial step size t, and the fitness weight w in the objective Fitness R² complexity.
359
250
Time Completityys
Equation Complexity
0.590
12
Equation Complexity
12
-260
0.9
0.925
12
139g
0.000
0.94
035
160
R'
R2
Complesty
0.55
10
20
30
56
66
70
03
67
13
13
15
0.770
275
360
325
350
375
400
425
16
Figure Ablation results on CMA-ES hyperparameters. Left: effect of population size s on R and runtime. Middle: effect of initial step size R and equation complexity. Right: effect of weight in the fitness function.
Population size s. As increases. R² steadily improves from 0.95 to nearly 0.99 but saturates after = 50, indicating that larger populations enhance global exploration yet exhibit diminishing returns. Meanwhile, runtime (i.e., time complexity) grows almost linearly from tens of seconds to nearly 300 seconds, showing a significant computational cost. Thus, medium-sized populations strike a better trade-off between accuracy and efficiency.
Initial step size l. The value of R² rises and then falls with increasing I, achieving the best performance around t 1.1, while equation complexity grows monotonically. Small step sizes constrain the search within local regions, whereas excessively large step sizes cause unstable convergence. A moderate there- fore provides preferable balance between accuracy and complexity.
Fitness weight w. the fitness function Fitness = R²2 - w complexity, w governs the trade-off between accuracy and complexity. Smaller emphasizes maximizing R² but yields overly complex equations, while larger w favors simpler forms at the cost of accuracy. A moderate choice of w effectively balances these two aspects, avoiding both overfitting and underfitting.
Summary. Overall, our ablations show that optimal performance is obtained with a 5 512-dimensional latent space, k 256 for covariance updates, and moderate settings for and These choices consistently yield robust trade-offs between accuracy, complexity, and runtime.
MODEL ARCHITECTURE AND TRAINING CONFIGURATION
F.1 MODEL. ARCHITECTURE
As illustrated in Fig. GenSR consists of input embedders (for numerical values and equations), shared Transformer encoder, posterior and prior MLPs, a feature fusion linear layer, and a Transformer decoder.
20

---

Published as a conference paper ICLR 2026
The numerical input {(x, y) E I=x((++1) is first tokenized into {(,, y) E Ixxx3(++1) then mapped by the numerical embedder to X, E R=xx(D++1xdd, where d, 64 denotes the numerical embed- ding dimension. A reshape and linear projection yield the final X E Im Similarly, equation expressions are padded to the maximum sequence length m and embedded into F E IRmxd with 512.
For the posterior branch, X and F are concatenated and passed through an 8-layer Transformer encoder heads, 512 hidden units), followed by an MLP produce the posterior distribution NV ( at). For the prior branch, X alone (with masking applied) is fed into the same encoder to produce the prior distribution (H22 oï). Notably, the numerical input X does not include positional encodings, ensuring that only the functional relations, rather than ordering artifacts, are modeled. The KL divergence between posterior and prior distributions serves as a regularization term:
CKL = DKL DK(,, (l oī)).
During training, samples are drawn from the posterior (M,, d to form Zi E Rmxx, which is pro- cessed by the feature fusion linear layer and decoded by symmetric 8-layer Transformer decoder (8 heads, 512 hidden units). The decoder generates logits of symbolic expressions, which are compared with the ground-trum truth equations through a reconstruction loss:
172
LCE
lo po (fz z1),
where fi denotes the i-th token of the target equation and the decoder parameters.
The overall training objective combines these two terms:
L Lce + L
where À is balancing hyperparameter. We employ KL. annealing to gradually increase from 0 to its full value during training. This prevents the KI. term from dominating the optimization in the early stages and mitigates posterior collapse, while ensuring a well-structured latent space in later stages.
During the inference phase, we use the prior distribution (N,, o form latent variables Zg Rm which are then fed into the Transformer decoder to generate candidate equations. Each equation is evaluated with fitness score, which subsequently utilized by the CMA-ES optimization procedure iteratively refine the search and obtain the final predicted equation.
F.2 TRAINING CONFIGURATION
We adopt the Adam optimizer with batch size of 256 and train the model for 200 epochs, each consist- ing of 1000 steps. The learning rate set 1.0 x 10-3 and follows the Noam schedule with w 8000 warm-up steps. This schedule allows the learning rate to increase linearly during warm-up and then decay proportionally to the inverse square root of the step number. For the KL divergence weight X, we apply an an- nealing strategy: À linearly increased from 0 to 1.0 during the first 50% of training steps and kept constant thereafter. This prevents the KL term from dominating optimization at the beginning and mitigates poste- rior collapse, while ensuring structured latent space later stages. The model comprises 162.1 million trainable parameters. All experiments are conducted using two NVIDIA H800 GPUs (Hopper architecture, 80GB memory each), with each epoch requiring approximately hour of training time.
21

---

Published as a conference paper ICLR 2026
G DETAILED EXPERIMENT RESULTS
G.1 DETAILED DESCRIPTION OF EVALUATION DATASETS
For our experimental evaluation, we adopt the SRBench benchmark (La Cava al. 2021). widely used and challenging collection of datasets for symbolic regression. SRBench consists of three representative groups: the Feynman dataset, the Strogatz dataset, and the Black-box dataset. These datasets differ in input dimen- sionality, sample size, functional properties, and task difficulty, jointly forming a comprehensive evaluation framework for symbolic regression methods and providing standardized platform for comparison.
Feynman dataset. The Feynman dataset (Udrescu et al. 2020) derived from the Feynman Lectures on Physics and contains 119 equations spanning diverse domains such mechanics, electromagnetism, and quantum physics. These equations are structurally diverse and grounded in real physical laws, making them key benchmark for evaluating symbolic regression in scientific discovery. To ensure consistency across tasks, the dataset constrains the input dimensionality to D < 10. A unique advantage that the true underlying functions are fully available, eliminating the ambiguity present in black-box tasks. With a cumulative size of approximately 10 data points, this dataset allows us to assess both large-scale regression accuracy and generalization across diverse physical laws.
Strogatz dataset. The Strogatz dataset (Strogatz, 2018) originates from nonlinear dynamical systems, in- spired by the classical ODE-Strogatz collection. Each task corresponds system modeling problem, with input dimensionality fixed at D = 2. Unlike Feynman, which emphasizes physical interpretability, the Strogatz dataset focuses on dynamic behaviors and is particularly suited to test- ing the ability of methods to handle nonlinearity and temporal dependencies. Each sub-dataset consists of approximately N = 400 samples, with true underlying functions provided to enable evaluation within a well-defined dynamical framework. This dataset is thus critical for benchmarking performance under small- sample and nogginrtyin conditions.
two-dimensional dynamical
Black-box dataset. The Black-box dataset (Olson al. 2017) highlights the applicability of symbolic regression to complex real-world scenarios. This group is primarily sourced from the PMLB repository and has been widely adopted in SRBench. Unlike Feynman and Strogatz, the target functions in this dataset unknown, mimicking real-world regression tasks. These datasets often contain higher levels of noise and uncertainty, making them effective test of robustness. To ensure comparability. input dimensionality is limited to D < 10. The collection comprises 57 sub-datasets, covering both real-world and synthetic cases, with dataset sizes ranging from just a few dozen to tens of thousands of samples. This wide variation tests not only the accuracy of algorithms but also their stability and generalization under noisy conditions.
summary, the Feynman dataset emphasizes physical interpretability, the Strogatz dataset captures the complexity of nonlinear dynamical systems, and the Black-box dataset stresses robustmess under real-world uncertainty and noise. Together, these three datasets form complementary and multi-faceted evaluation environment for symbolic regression. Under this benchmark, we are able to comprehensively examine trade-offs in accuracy, efficiency, and complexity, while highlighting the robustness and generality of our proposed approach across diverse tasks.
G.2 EVALUATION METRICS
To evaluate the performance of symbolic regression methods, we consider three key metrics that capture accuracy, efficiency, and interpretability:
22

---

Published as conference paper ICLR 2026
Coefficient of determination (R²)): This metric quantifies how well the discovered symbolic ex- pression explains the variance of the target data. It is formally defined as:
---00= (M: 0=)²
R²=
(2)
(3: 1 g)²
where yi denotes the ground-truth values, y: the predicted values, and y the mean of the ground- truth values. A higher It² indicates better predictive accuracy, with R2 representing a perfect fit and R² 0 implying improvement
the predictor.
This metric measures the wall-clock time required for the algorithm to complete the symbolic regression task and output final expression, reported in seconds (s). reflects the computational efficiency of the method, where shorter search times indicate faster convergence and reduced resource consumption.
Time complexity:
Equation complexity: This metric evaluates the interpretability and compactness of the discovered expressions. We define complexity as the total number of tokens in the final equation:
m2
Complexity(f) Σ I(t)).
where f denotes the discovered equation, [ty] are its tokens (including coefficients, operators, and variables), and I(.) is an indicator function that counts the presence of each token. Lower complexity corresponds to simpler and more human-readable expressions, which are less prone overfitting and closer to the underlying physical or mathematical laws.
G.3 BASELINES
Evolutionary algorithm- symbolic regression. We include several genetic programming variants:
GPlearn (Stephens et al. (2015), a scikit-leann-styye GP framework; AFP (Schmidt & Lipson| (2010), which maintains diversity via age-fitness Pareto optimization; ITEA (de Franca Aldeia, 2021), which exploits tranactionntaansfoonatiion operators; EPLEX (La Cava et al. 2016). based on epsilon-lexicase selection; MRGP (Amaldo et al., 2014), combining multiple regression with GP: SBP-GP (Virgolin et al. 2019). using semantic backpropagation with linear scaling; GP-GOMEA (Virgolin et al. 2021), model-based GP with efficient linkage learning; Operon (Burlacu et al. 2020), a high-performance C++ SR framework; and FEAT La Cava al. 2019), which evolves networks of trees learn concise regression models.
Deep/Neural and reinforcement learning-bass symbolic regression. We further consider neural and RL-driven methods: SPL (Sun et a..). which employs Monte Carlo tree search for physics discovery: DSR (Petersen et al.). a reinforcement learning method with risk-seeking policy gradients; RSRM (Xu et al., 2024), reinforcement symbolic regression machine; TPSR (Shojaee al. (023). which formulates sym- bolic regression as a Transformer- based planning problem; E2ESR (Kamienny et al. 2022). an end-to-end Transformer framework; SNIP (Meidani et al.). which unifies symbolic and numeric pretraining; RAG-SR (Zhang et al., 2025a), leveraging retrieval-augmented generation; MDL.frmer (Yu et al. 2025), guided by the minimum description length principle; and NeurSR (Biggio al., 2021). scalable neural symbolic regression approach.
Physics- and probabilistic-inspiee methods. We also evaluate AI Feynman (Udrescu al. 2020 which exploits physics priors and graph modularity for symbolic discovery. and BSR (Jin et al. 2020). Bayesian symbolic regression approach that provides uncertainty estimation.
23

---

Published as conference paper ICLR 2026
Conventional machine learning baselines. Finally, we compare against strong non-symbolic basclines, including tree-based ensemble methods: XGBoost (Chen & Guestrin, 2016). LightGBM (Ke et al. 2017) AdaBoost (Freund & Schapire 1997), and Random Forest (Rigatti) 2017).
G.4 FUNCTION FAMILY CLASSIHICATION TASK
0.849
0.8
0.723
ACC
0.7
0.661
0.6
E2ESR
SNIP
Ours
Method
Figure 9: Accuracy comparison of different methods on the classification task. GenSR (Ours) achieves the highest accuracy (0.849), significantly outperforming E2ESR (0.661) and SNIP (0.723).
To further examine the transferability the learned generative latent space beyond symbolic regression, we design a function family classification task. The goal of this experiment is to assess whether the latent space representations learned by GenSR capture not only structural continuity and numerical smoothness, but also discriminative features that support downstream classification.
Dataset construction. We synthesize a dataset covering five representative function families: exponential, logarithmic,, sine, cosine, and tangent. To increase diversity and difficulty, each family is instantiated under four different input dimensionalitics (2D, 3D, 4D, and 5D), leading to a total of 5 x 4 = 20 classes. For each class, randomly generate 2000 samples, where each sample consists of input points and their corre- sponding function values. The dataset is split into training, validation, and test sets with a ratio of 3:1:1 to ensure fair and stable evaluation.
Model and training setup. We adopt the GenSR encoder pretrained on symbolic regression, which pro- duces 51 512-dimensional latent vectors. On top of the encoder, attach a lightweight classification head composed of two fully connected layers: the first layer maps the latent vector to hidden representation with ReLU activation, and the second layer outputs 20-dimensional logits vector followed by softmax layer to obtain class probabilities. We fine-tune the entire network using cross-entropy loss with the Adam optimizer, a learning rate of 10-3 and a batch size of 128. Training is run for at most 100 epochs, with early stopping if no improvement is observed on the validation set for five consecutive epochs. At test time, the class with the highest probability is selected as the prediction, and overall accuracy is reported.
Comparison and results. For fair comparison, apply the same classification head design, data split, and training procedure to encoders of E2ESR (Kamienny et al.) (022) and SNIP (Meidani et al.). As shown in Fig. [ GenSR achieves significantly higher classification accuracy than both baselines. This result demonstrates that the latent space learned by GenSR is not only generative and structurally continuous, but also discriminative and generalizable, providing strong foundation for downstream tasks beyond symbolic regression. In other words, GenSR yields representations that are simultaneously useful for generation, search, and classification, highlighting its potential for broader applications.
24

---

Published as conference paper ICLR 2026
In the future, the latent space obtained by the CVAE also can be applied to a variety of domains, such as biological contexts (Duan et al. 2025), psychological studies (Jia et al. 2025), energy systems and meteo- rological forecasting where spatiotemporal patterns play critical role (uu et al. 2024 2025), and broader
areas of artificial intelligence where interpretability and robustness are of particular importance (Li et al.
2024 2023b
Huang al. 2026 Feng et al. 2026 2025
Kang et al. 2025 Zhang et al. 2025b Yang
et al.
G.5 COMPLETE NUMERICAL RESULTS
6 and 7 report the performance of different methods on the Strogatz (Strogatz 2018) and Feynman datasets (Udrescu al., 2020) under varying noise levels, while Table summarizes the results on the Black-box datasets (Olson et al.) 2017). Across all these benchmarks, our method demonstrates a consistently strong balance among the three key metrics-accurace (R²). search efficiency (time), and equation complexity. Notably, even under noisy conditions, our approach maintains high R scores with smaller variance, achieves shorter runtime, and yields simpler equations, highlighting robustmess and overall superiority compared to existing baselines.
Tables
25

---

Published as conference paper ICLR 2026
Table Complete experimental results of black-box dataset.
Method
Test R² t
Complexity 4
Time (s)
Operon
0.7945
65.69
2974
GP-GOMEA
0.7381
30.27
9636
MDL
0.6258
29.88
541.7
DSR
0.5625
9,465
36852
SBP-GP
0.7869
634
149344
FEAT
0.7621
82.49
6432
EPLEX
0.7372
53.14
15796
AFP-FE
0.6400
36.04
6184
AFP
0.6333
34.89
6033
GPlearn
0.5390
19.06
24254
Linear
0.4437
17.4
0.2447
NeurSR
0.1228
13.33
11.7
XGB
0.7496
20186
236
AdaBoost
0.6939
9481
65.12
LGBM
0.6410
5734
29.9
ITEA
0.6295
116.7
12183
SNIP
0.3335
38.91
3.286
BSR
0.2725
22.52
59822
RandomForest
0.6615
1517178
120.4
KemnllRide
0.5952
1824
39.19
FFX
0.5575
1562
244.3
E2ESR
0.3612
61.09
7.101
MRGP
0.5300
10802
165007
MLP
0.5238
3882
30.49
AlFeynman2
0.2110
2240
86854
SPL
0.5472
12.96
2047
RSRM
0.3324
9.23
1752
Ours
0.8416
35.05
2725
26

---

Published as conference paper ICLR 2026
Table 4: Complete experimental results at noise level e = 0.0. The results include test set R², search time, and formula complexity of different methods in the Feynman and Strogatz datasets.
Strogatz Dataset (e = 0.0)
Feynman Dataset (E = 0.0)
Method
R²1
Time (s)
Complx
R2t
Time (s)
Complx
FEAT
0.9210 (02004))))
1135
119(210)
0.9190
2417
205.3
NeurSR
0,5206
14.82
11.3
0.3958 (0003))
24.18 14u2)
13.28
E2ESR
0,5341 (20043)
3.729
32.49
0.8570
4.062
36,021z0a8)
SNIP
(.9.52(s0.02))
509
29
0.985
1203.5 (2
31.71
GPlearn
0.7689 (60.667)
1195
28.96
0.8809(cm. (2007)
3900
72.43
AFP
0.9248 (00041)
192.2(+106)
38.17 (
0.9590
3655
36.87
AFP-FE
0.9442 (60.045)
1104 I (814277))
46.16
0.9806 (+0
17817 (e530)
39.97
EPLEX
1.81225m (e0006))
548 (
50.09
0.9869
(00006)
12771 (+4908)
52.95
SBP-GP
0.98 12 (50006))
1759 (1176))
712.2
0.9945
(200001)
2890 (237)
489.4:a10)
GP-GOMEA
0.9925 (20.009)
(760g)258)
36.43
0.9956::0000 003)
5030ka96)
34.57
Operon
0.9878 ( (200023)
66.58ta10)
59.23
0.9889
(00066)
2174
69. (9.88
SPL
0.7390
322. 122.1aata (1100)
14.55
0.7073
(2001))
209
12.88
DSR
0.7602 (60.086)
1858
15.6
0.844
(s0091)
1733
14.86en
RSRM
0.5501 (40)033)
121.9 121.9(43)
13.09
0.8003
(20013)
116.3
13.17
AlFeynman2
0.6459 (e0.009)
762. (1424)
22.26
0.93 4 (6006)
854.3
124.5
BSR
0.8455 (60.044)
31380
38.98
0.6609
(e0 018)
29065 (5765)
25.5(H023)
MDL
0.9900(za09))
186.6 (1355)
14.07
0.917
(200055)
467.3 467.3
23.4
TPSR
0.9655
291 (453)
56.2(+1.4)
0.992
(=0
236
57.27
RAG-SR
0.99 119114
1477
46.5
0.9926
3156
46.5 ((114)
Ours
0.9918 (0911800mm)
33
20.53
0.9872
192.3(+21)
22.94
27

---

Published as conference paper ICLR 2026
Table Complete experimental results at noise level 0.001. The results include test set R2, search time, and formula complexity of different methods in the Feynman and Strogatz datasets.
Strogatz Dataset (e 0.001)
Feynman Dataset (e 0.001)
Method
R
Time (s)
Complx
R²4
Time (s)
Complx
FEAT
0.9244
(4181)
106.7wis
0,9207
1726
196.5ta1
NeurSR
0,5219
(s001) (s0031)
15.07
11.41a030
0.3979Gu0013
24.51 (4t7
13.24G0.1)
E2ESR
0.5105 (0.06))))
3,436:20.75)
33.83
0.8585
3.894
35.85
SNIP
0.9667 (s0022)
271.941184)
27.57
0.9904
1365.07 (0000
31.43
GPlearn
0.7955
913.3
29.59
0.8902
3316
60.49
AFP
0.9172 e00052)
143.4 4
38.75
0.9606 =0.00))
371 (645))
39.33
AFP-FE
0.9447 4 (20022)
8 8108
48.74ka3
0.9805 (20.007)
26160ta157)
46.47
EPLEX
0.8488
(0005))
416. 1 (181)
49.26
0.9866 +0.007)
1234 (1430)
56.03
SBP-GP
0.9879
(00015)
19596
820.5
0,9953 (200001)
28940
574.4tam
GP-GOMEA
0.9914 4 (0C)
804
41.14
0.9962 (0.0001)
2904
45.23
Operon
0.9843
3 (0033))
75.6 (5.68(144)
67.03
0,991 6 H0U015)
2195 (3404)
69.67 (2100
SPL
0.7526 ( (0LD4)
(4.4
0.7073 3 (+00116)
275.5
13.15
DSR
0.8067
500.8
18.66122.11
0.8764
830. (
16.04g030
RSRM
0.5447 (05)
129.76 (+8.1)
(2.23
0.8 I (44
128.2 (
13.04
AlFeynman2
0.6855
(00091
84.19 (1744)
25.64ka3)
0.9177 (+0.018)
638
130.6
BSR
0.8224
(21)
24299 (+378)
37.68
0.6538 (+00223)
30255
25.85
MDL
0.9965 (s0004)
171.5
21.38
0.9079 (s0.008)
428.9(221)
32.35(a1.))
TPSR
0.921 ((002))
193.22
55.64
0.992 (+0.001)
190.2
59.75
RAG-SR
0.9908
591.44438.4
49.29
0.9917 7.99977.0.000) (10.003)
3590
49,4 49.42(a1.)
Ours
0.9951 (#0007)
41.5600.1
20.45
0.9883
279.31
23.5
28

---

Published as conference paper ICLR 2026
Table 6: Complete experimental results noise level 0.01. The results include set Ri, search time, and formula complexity of different methods in the Feynman and Strogatz datasets.
Strogatz Dataset (e 0.01)
Feynman Dataset (e 0.01)
Method
R²1
Time (s)
Complx
R²t
Time (s)
Complx
FEAT
0.9244c000m)
(72.9(20)
95.61a1
0.9212
1464
167.1esss
NeurSR
0.5179k2000
15.48(a1.))
11.63ra0.0
0.3942
24.62/a1.1)
13.26
E2ESR
0.503 (50.04))
3.392
35.94
0.8345
4.274
40.07
SNIP
0.9844x20025)
237.6 (432)
29.35
0.987
1210. (4150)
32.62(a1.3)
GPlearn
0.7956 (e0 (50.059)
907.4 10)
30.59(143)
0.8890 (0,009)
3351 (1437)
60.07
AFP
0.9153 (-0,55)
152.2(H15)
38.62
0.961 (0.000)
4090
40.86x1.D
AFP-FE
0.9582 (19581sa000
8898
48.8 48.8:a54)
0.981 9 (40.008)
27763
46.92(22.))
EPLEX
0.8562
(+0.040)
437.6 (337.6+644
53.07
0.991 0 (++00022
11043 (
54
SBP-GP
0.9813
(af.015)
20783
850.9
0.9950,a..0m (0.001))
28954ka15)
595 (
GP-GOMEA
0.9783 (s0,029)
765.6x41005)
42.64u4.5)
0.9967 (.0.001)
3020
44.67(aL.
Operon
0.9829
(30.33))
94.92(s15)
81.68
0.9878 8 (00.000)
3165
87.96x+1.)
SPL
0.7388 (60.000)
413.3
14.71
0.7133
295. (e175)
13.43
DSR
0.8199 (+0.05))
492.9
18.51
0.8782 (0.00))
929.8 (
16. 16.2
RSRM
0.5969 (el.077)
142.2(
14.22(+15)
0.8092 (+0.15) (000155)
131.6;a14)
12.97
AlFeynman2
0.7753
(20047))
85.17(w55)
32.41caAm
0.8732(20.00)
629.4.as
155.2(48)
BSR
0.8127 (20,070)
23622
38.74u2
0.6734
3041 LE47LD
28.03
MDL
0.97 18 (20055))
505.
20.3 1 (1881)
0.9140
844.3
31.15(41.5)
TPSR
0.9798 (s0.008)
175.09
56. 4
0.991 (2000000)
217.06u30
64.27uu2)
RAG-SR
0.9867 (+0.019)
509(+am
49.43(+17)
(.9905
3241 (1417)
72.4 12.41(22)) (42.2)
Ours
0.99361200m
44.58
20.42(u0.))
0.9872:s0.002
486.834448
23.36
29

---

Published as conference paper ICLR 2026
Table 7:
Complete experimental results noise level e 0.1. The results include test set Ri, search time, and formula complexity of different methods in the Feynman and Strogatz datasets.
Strogatz Dataset (e 0.1)
Feynman Dataset (e 0.1)
Method
R2t
Time (s)
Complx
R²t
Time (s)
Complx
FEAT
0.9228
(173))
84.2/s19)
0.9195
777.9tmm
99.48ua..)
NeurSR
0.5054200ss
(..47
0.3823 12.74
25.8 Laio)
13.54
E2ESR
0.5152
5.62 (4449)
38.49u27)
0.77 (.0004)
6. 6.183
44.09
SNIP
0.9187 (20.012)
5650.7us3
39.35u2.3
0.99 7 (20.020)
65 6514.8uss
37.7 ((2.0)))
GPlearn
0.8228 (s0.05))
894.6(s108)
25.84
0.891 (s0.007)
2938
48.83
AFP
0.91 0 (00065)
161 44
(4.44
0.9577 (+.0007)
3886
40.79(H15)
AFP-FE
0.9496
10082
50.96
0.9826
28812
(8.87(414)
EPLEX
0.8822 (+00078)
53.46
0.9901 1 (+0.004)
10283
45.62(+1.)
SBP-GP
0.9323
(a00080)
21886
901. (434)
0.9905 (+0:007)
28932
621.9(c6.))
GP-GOMEA
0.9668
(.966) (30033))
402.9(us0))
43.7lG2.
0.9957 (A0.003)
3186
46.43
Operon
0.9380
97.13ke15)
(3.44
0.9847 (10.008)
3090
89.23
SPL
0.7715 (+0044)
355.3(s59)
13.91
0.7109 (00.013)
270
13.54
DSR
0.8086 (+00047)
500
18.51 (+2.2)
0.8779 (0.0002)
814.9
16.03
RSRM
0.5553 (+00577)
138.3(+6.4)
13.54(a1.1)
0.81 (55104(H0..00 (+0:016)
133.8
12.8 1 (H0.41)
AlFeynman2
0.3170 0.3170(20.082) (20.882)
65.97(219)
23,53
0.2248
aman))
710.7 (anan)
176.6 6
BSR
0.7190
23292
49.54
0.6567
(A0:02)) (10:024)
32497 (L7914)
28.77
MDL
0.9686
522
20.5
0.9097
(20.009)
962.9 (918))
30.86
TPSR
0.9707 (20.022)
159. 1a.3)
56.28
0,9836 (0.005)
184.0 (A8.3)
66.96va1.9
RAG-SR
(-0043))
385.3(e37)
46.23(+2.)
0.9849((0.006)
2957 (+366)
72.12
Ours
0.9773ta0.000
257.64
20.2 (4.3)
0. 0.9886.10.00s)
618.9
23.62
30

---

Published as conference paper ICLR 2026
H COMPARISON WITH THE STANDARD BAYESIAN OPTIMIZATION (BO) ALGORITHM ON FEYNMAN DATASET
Table 8: Comparison of CMA and standard BO under different & with 512-dimensional latent space
Method
R²²
Time (s)
Complx
CMA
0.9872 0.9274
192.3
22.94
k=256
standard BO
295.8
21.37
CMA
0.9675
179.1
23.26
k= k=128
standard BO
0.9161
237.5
22.55
CMA
0.9388
173.2
23.93
k=64
standard BO
0.8974
165.6
23.07
Table 8 reports the comparison between CMA-ES and standard Bayesian Optimization (BO) GenSR with the 112-dimensional latent space. our main experiments, GenSR adopts k 256, corresponding to the top-k principal latent dimensions with the largest variance. We therefore first conduct BO experiments under this same setting. As shown ES, but its performance in terms of accuracy (IC²) and optimization time is substantially worse. This may be attributed to the challenges that BO potentially faces when dealing with nigh-dimensional optimization problems.
table, BO attains slightly lower expression complexity than CMA-
To further investigate the effect of reducing the number of optimized dimensions, we additionally evaluate k 128 and k 64 (a setting already examined for CMA-ES in the k ablation study reported in Appendix 7p. Although decreasing & improves the runtime BO, its accuracy does not increase accordingly. At the same time, the accuracy of CMA-ES drops significantly under smaller & values. These results suggest that overly small & leads to an excessive loss of important latent dimensions, thereby limiting the expressiveness of the search space and degrading the performance of both optimization methods.
Table 9: Comparison of CMA and standard BO under different i with 256- 56--dimensional latent space
Method
R
Time (s)
Complx
CMA
0.6984
181.4
18.23
k=256
standard BO
0.6241
206.6
17.12
CMA
0,4935
174.5
19.08
k=128
standard BO
0.4797
178.3
18.67
We further evaluate another pretrained model used in Appendix E.1 namely the GenSR variant with the
256-dimensional latent space, and conduct experiments with 256 and = 128. As shown in Table 9 both CMA-ES and standard BO exhibit relatively low accuracy under this lower-dimensional latent rep- resentation, which may suggest that the expressive capacity of the latent space is notably reduced such settings. Moreover, the performance gap between BO and CMA-ES becomes less pronounced.
Overall, these observations imply that BO may face certain limitations when dealing with figh-dimensional optimization, whereas CMA-ES appears relatively more robust for task. Attempts improve BO by ducing the latent dimensionality or decreasing the number of optimized dimensions may still leave accuracy constrained, potentially due to insufficient latent expressiveness the omission of important dimensions.
31

---

Published as conference paper ICLR 2026
I ABLATION ANALYSIS OF THE POSTERIOR SYMBOLIC ENCODER BRANCH
Table 10: Comparison between GenSR and its disrupted variant
Method
R²1
Complx
GenSR
0.921
19.04
GenSR w/o Posterior
0.752
26.41
To examine the role of the dual-branch framework shaping the latent space, we construct variant that explicitly disrupt the posterior branch, denoted as GenSR w/o Posterior. While keeping the overall two- branch architecture intact, we pair numerical samples X with randomly sampled equations F as inputs, and replace the reconstruction loss of the posterior branch with cross-entropy loss computed against the target expression corresponding to X. In addition, to directly evaluate the inherent performance of the model under this setting, we remove the CMA-ES optimization step and rely solely the initial localization produced by the pretrained model.
As shown in Table disabling the posterior branch leads to substantial drop in performance (R² decreases from 0.921 to 0.752, and the complexity notably increases), indicating that the underlying latent-space struc- ture is significantly disrupted. Consequently, the initial localization becomes unreliable, demonstrating the critical role of the dual-branch framework in maintaining well-structured latent representation and overall model performance.
J COMPARISON WITH THE STANDARD BAYESIAN OPTIMIZATION (BO) ALGORITHM
KL Term and Weight (1) During Training
0.8
Weight (2)
KL Term
0.6
20-
8.
10-
0.2
0
25
75
100
125
150
175
200
Epoch
Figure 10: The horizontal axis denotes the epoch, while the left and right vertical axes correspond to the KL term and the KL weight, respectively. The KL weight annealing strategy used here increases the weight from 0 to 1.0 during the first 50% of the training epoches and keeps it constant thereafter.
32

---

Published as conference paper ICLR 2026
Posterior Mean Distribution
Prior Mean Distribution
210
.23
300
210
255-
Fremoncy
200
150
100
203
0.29
0.62
0.30
0.53
0.0
5.50
0.50
5.2
0.0
3.23
6.00
5.00
0.50
060
0.30
0.5
Moa
Mean
Posterior Variance Distribution
Prior Variance Distribution
Frequnacy
Fremmenc
100
0.2
0.5
2.0
2.6
022
8.0
20
1.5
Parlonce
Varlance
Figure 11: Estimated distributions of prior and posterior means and variances obtained from 5,000 randomly sampled inputs in GenSR.
As described in Appendix F.2 10 the KL weight gradually increases from 0 to I in the first half of training, while the KL term eventually settles into relatively stable range in the later stages, indicating that no collapse occurs.
adopt a KL weight annealing strategy during training. As shown in Fig.
addition, GenSR, the prior distribution is learnable rather than fixed to a standard Gaussian. During training, the continual updates of the prior prevent the posterior distribution from trivially collapsing toward it, which may help reduce the likelihood of posterior collapse. Fig. 11 shows the distributions the prior and posterior means and variances computed from 5,000 randomly generated samples using the final pretrained model. Clear discrepancies are observed in both the mean and variance distributions, indicating that the prior and posterior branches have collapsed into single distribution. If posterior collapse had occurred, the posterior would converge toward the prior; however, the differences shown here demonstrate that the model maintains non-degenerate and informative posterior.
K COMPARING GENERATIVE AND DISCRIMINATIVE LATENT SPACES
Contrastive learning methods (e.g., SNIP) construct a discrininative latent space that is fragmented and contains numerous "holes," such that latent vectors sampled between two valid points often decode empty or otherwise invalid strings. In contrast, our GenSR constructs a generative space (a continuous manifold), meaning we sample continuously and decode valid equations.
It learns the distribution of equations This allows for "smooth interpolation." which is a prerequisite for the fine-grained search (CMA-ES) we perform. Table K provides clear summary of the differences between the generative and discriminative
33

---

Published as conference paper ICLR 2026
latent spaces. Besides, the interpolation experiments (Appendix C.2 confirm that SNIP often falls into semantically invalid areas during search.
Table 11: Comparison between Generative latent space and Discriminative latent space
Property
Generative latent space (GenSR)
Discriminative latent space (contrastive)
Learning Objective
Equation distribution Supported High
Pointwise alignment
Continuous latent sampling Suitability for search
Not supported
Low: limited with invalid regions
34


Reviews:
[{"review_id": "AqaMciFXkU", "reviewer": "Reviewer_mtFG", "summary": "This paper presents GenSR, a generative SR method that trains a conditional generative model that, at test time (inference/runtime), starts from a good latent-space initialization and refines it with CMA-ES to produce a set of good-fitting equations. The paper claims to introduce a new search paradigm, a Bayesian perspective for SR, and experiments that demonstrate state-of-the-art performance.", "questions": "* Can the authors clearly delineate how their work is novel, new, and different, compared to the existing highly similar papers such as [1,2,3] stated above, and the broader aspect of VAEs for SR, and similarly other generative models for SR [3]. As it stands with the current submission, this paper does not propose anything novel in context to the existing literature [1,2]; and has severe experimental issues.\n* Given that the authors refine constants using BFGS as an inner optimization, have the authors considered the functions getting stuck in local minima and strategies to prevent this, as routinely done in other SOTA SR methods [4]?\n* Can you compare experimentally against [4,6,5] by running these competing methods against your method?", "limitations": "", "rating": null, "confidence": 5, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "* The paper is straightforward to read and follow.\n* SR is an essential problem for the scientific ML community.\n* The interpretability and comparison plots were informative, and a nice visual comparison.", "weaknesses": "* Core contributions are not novel, and the authors of GenSR have significant overlap with existing related work that was not discussed at all [1]. [1] proposes a generative SR method using a conditional VAE (CVAE), which encodes expressions into a smooth low-dimensional space, and then applies an evolutionary search algorithm to this latent space, to efficiently sample equations, and at the time of its publication claimed state-of-the-art. Given that the contributions of GenSR are exactly these, I find this missing key related work in direct conflict with this submission. I encourage the authors to: 1) clarify how their method differs, and 2) clearly delineate what is novel and new in their paper in the context of all existing VAE SR methods [1,2] and existing generative methods [3].\n  * Using a VAE for SR is not novel and has been extensively done in prior work [1,2], which is conveniently not discussed in this paper. Furthermore, the other missing key related work [2] already shows that using a VAE can improve equation generation under noisy settings, a result also reported in this paper.\n* ``GenSR provides a new Bayesian perspective for SR''; this contribution is lacking without appropriate citation and delineation to existing generative SR methods that also provide a Bayesian perspective for SR, such as [3].\n* Claim of \"State-of-the-art performance\" is not validated, as SRBench and corresponding results are old and outdated, dating from 2021/2022. I would encourage the authors to implement the most recent SOTA SR methods and compare them against those to make such a bold claim. Furthermore, the results appear to be cherry-picked, as SR Bench has many problem sets, and only the subset of test equations (Feynman) results are presented graphically in the main figure; whereas other newer SOTA SR methods on SRBench show the same figure across all the equations in the benchmark [5]. I encourage the authors to directly compare to [5] and more recent SR methods to make any claim of state-of-the-art. Furthermore, I encourage the authors to compare against standard problem sets such as Nguyen, Livermore, and similar problem sets as done in [6]. Moreover, the authors omit from the SRBench results the results for \"Symbolic Solution (%)\" and \"Accuracy Solution (%)\", leading to suspicion over the results, as these are computed as standard when running the benchmark, and choose to only present R^2 instead. I encourage the authors to be transparent and report all standard results from the benchmark, especially if making any state-of-the-art claims.\n* Appendix B seems not necessary, or novel, and appears to be a standard re-statement of the classic textbook ELBO derivation. I would encourage the authors to simply cite this, without the proof in the appendix, as it is not new.\n\nReferences:\n* [1] Mežnar, Sebastian, Sašo Džeroski, and Ljupčo Todorovski. \"Efficient generator of mathematical expressions for symbolic regression.\" Machine Learning 112.11 (2023): 4563-4596.\n* [2] Popov, Sergei, et al. \"Symbolic expression generation via variational auto-encoder.\" PeerJ Computer Science 9 (2023): e1241.\n* [3] Holt, Samuel, Zhaozhi Qian, and Mihaela van der Schaar. \"Deep Generative Symbolic Regression.\" The Eleventh International Conference on Learning Representations.\n* [4] Kamienny, Pierre-Alexandre, et al. \"End-to-end symbolic regression with transformers.\" Advances in Neural Information Processing Systems 35 (2022): 10269-10281.\n* [5] Landajuela, Mikel, et al. \"A unified framework for deep symbolic regression.\" Advances in Neural Information Processing Systems 35 (2022): 33985-33998.\n* [6] Mundhenk, T. Nathan, et al. \"Symbolic regression via neural-guided genetic programming population seeding.\" arXiv preprint arXiv:2111.00053 (2021).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wN6BYoFBT8", "reviewer": "Reviewer_sUgp", "summary": "This paper points out that \"similar structure ≠ similar value\" in discrete symbol space leads to noisy fitting error signal and lack of search direction, which is the source of inefficiency of most SR methods. Therefore, GenSR is proposed. The pre-trained dual-branch CVAE reparameterizes the equation into a generative latent space with \"global symbolic continuity + local numerical smoothness\", according to the \"map construction → rough positioning → fine search\" paradigm and improves CMA-ES search. SR is reformulated as a Bayesian problem and implemented by ELBO training.", "questions": "1. What are the advantages of methods trained via CVAE over methods trained via contrastive learning? \n2. What's the reason? It seems that contrastive learning can achieve similar results to the map building mentioned in this article, and I think it might even be better.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. This method seems to achieve very good experimental results\n2. This paper experiment is relatively complete, considerable workload\n\n3. In addition, the authors give a Bayesian perspective and reformulate SR as a probabilistic inference problem of maximizing P(F|X) trained by ELBO, which provides a clear theoretical framework and interpretability support for the rationality of the method.", "weaknesses": "The article is well written, but it seems to be only transferring the method in the multi-modal of graphics and text to the field of symbolic regression, and the innovation of the algorithm is not very high. Ask the author to emphasize his technical innovation", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "4ot29aaEyi", "reviewer": "Reviewer_CAH6", "summary": "The paper proposes GenSR, a Conditional-VAE-based symbolic regression model that encodes equations into a continuous latent space and performs CMA-ES optimization there. It claims improved performance over existing symbolic methods such as E2ESR and SNIP.", "questions": "1. The paper should clarify why encoding the symbolic expression is required. In a CVAE, the encoder is necessary because many outputs can correspond to the same condition, such as images. However, in symbolic regression, for a given dataset, there are usually only a few valid equations. Please consider adding an ablation where the symbolic encoder/posterior branch is removed, to test whether the bottleneck itself is the key contributor to performance.\n\n2. Figure 5 illustrates numerical continuity, but the paper does not explain how this continuity emerges, since no explicit loss term enforces it. Please provide an explanation of this phenomenon. While it is expected that a VAE can cluster similar functions, the observed continuity is an interesting behavior that deserves discussion.\n\n3. On page 21, one citation appears incorrect and is displayed as a question mark. Please correct this.\n\n4. The model weights and source code should be released for reproducibility.\n\n5. Some methods listed in Table 4 are not reported in Table 3. Please include their performance on the black-box dataset for completeness.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "A key strength of the paper is the clear comparison with SNIP, which highlights the limitations of contrastive learning in SNIP and demonstrates the advantage of using a VAE framework for symbolic regression. Although using autoencoders for symbolic regression has been explored in the genetic programming literature [1], this paper presents a meaningful improvement by extending the autoencoder to a pre-training scenario and showing the benefits of a conditional VAE design.\n\n[1]. Wittenberg, David, Franz Rothlauf, and Christian Gagné. \"Denoising autoencoder genetic programming: strategies to control exploration and exploitation in search.\" Genetic Programming and Evolvable Machines 24.2 (2023): 17.", "weaknesses": "The main weakness is the lack of an ablation study to justify the necessity of the posterior symbolic encoder branch.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tsd1BbPW4G", "reviewer": "Reviewer_A4iV", "summary": "This paper proposes GenSR, a novel symbolic regression framework that transforms the discrete equation search space into a continuous generative latent space using a dual-branch conditional variational autoencoder. The method follows a \"map construction → coarse localization → fine search\" paradigm, where the CVAE learns a latent space with global symbolic continuity and local numerical smoothness. At inference, a modified CMA-ES algorithm refines candidate equations in this latent space. The authors provide a Bayesian interpretation and demonstrate strong empirical performance on SRBench datasets.", "questions": "1. The interpolation experiment in Appendix C.2 is vague. Linear interpolation in high-dimensional spaces can be problematic. For a VAE with a standard Gaussian prior, most of the probability mass is concentrated in a thin shell far from the origin. Linear interpolation between two points is highly likely to traverse low-probability regions of the latent space.  Could the authors clarify the specific setup for this experiment?\n2. The paper reformulates the task as a Bayesian optimization problem. Have the authors considered applying Bayesian optimization algorithms directly in the latent space? I would be interested to see a comparison between CMA-ES and standard BO methods in terms of sample efficiency and the quality of the final solutions.\n\nPlease note that if the authors are able to address my questions during the rebuttal and discussion phase, I am willing to increase my score and confidence.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. This paper presents a highly promising idea by formulating the symbolic regression problem as an optimization task in a continuous latent space. This approach transforms the search space from a structured space and allows continuous optimization algorithms such as CMA-ES to be applied directly. While prior work like SNIP has explored latent space optimization, constructing a smooth and meaningful latent space remains an open question because it is crucial for effective latent space optimization.\n2. The authors propose a dual-branch conditional variational autoencoder to decouple and align symbolic and numerical information. This design is clearly motivated. Using KL divergence loss to align the prior and posterior distributions is a reasonable approach to promote local numerical smoothness.\n3. The experimental evaluation is thorough. The authors test on the SRBench benchmark and compare GenSR with 18 baseline methods. In addition to main results, the paper reports noise robustness, latent space t-SNE visualizations, interpolation experiments, and a detailed ablation study. Figures 4 and 5 demonstrate that GenSR achieves better structural disentanglement and numerical continuity in the latent space than competing approaches.", "weaknesses": "1. The architectural choice of using a VAE appears somewhat outdated, as the overall model seems to be a slightly modified version of E2ESR, with an additional layer appended to the encoder’s final hidden state to map to ($\\mu$, $\\sigma$). A well-known issue with VAEs is posterior collapse, where the latent variables become uninformative and the decoder effectively ignores them (this problem can be even more pronounced with high-capacity models such as transformers). It would be beneficial to provide empirical evidence (KL divergence during training) to demonstrate that posterior collapse is indeed being avoided.\n2. Although Section 3.3 provides an explanation for the desired symbolic continuity and numerical smoothness, its reasoning remains informal. For example, the paper assumes that the overlap of high-probability regions for similar equations enables smooth interpolation (lines 203–218). However, this is merely asserted with no empirical study of how frequently this occurs in practice. Similarly, the connection between minimizing the KL divergence and inducing local numerical smoothness (lines 219–234) is also asserted rather than formally established.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)