Symbolic Density Estimators for Unnormalized Distributions

# Symbolic Density Estimators for Unnormalized Distributions
Source: [https://arxiv.org/html/2610.10807](https://arxiv.org/html/2610.10807)
Vikas Kanaujiakvikas@iitk\.ac\.inRiyansha Singhriyansha@cse\.iitk\.ac\.inAffiliation:Department of Computer Science and EngineeringAffiliation:IIT KanpurShashank Sharmassharma23@iitk\.ac\.inAffiliation:Department of PhysicsAffiliation:IIT KanpurVipul Aroravipul\.arora@kuleuven\.beAffiliation:Department of Electrical EngineeringAffiliation:KU Leuven

###### Abstract

Estimating the symbolic or analytical form of probability density functions \(PDFs\) from observed samples is a fundamental challenge in statistical and computational modelling\. This process is critical for deriving interpretable and generalizable relationships characterizing the underlying phenomenon\. Traditionally, this estimation depends strongly on domain expertise and prior field\-specific knowledge, with experts selecting appropriate functional forms or parametric families based on empirical evidence and theoretical understanding\. The coefficients of these forms are then typically determined through parameter estimation\. In this paper, we develop a framework for estimating symbolic expressions of unnormalized distributions from observed samples using domain\-specific prior knowledge, such as the range of interactions and a predefined set of primitive functions\. We integrate deep generative models with symbolic regression \(SR\), incorporating inductive biases, such as factorizing large distributions, to keep the problem tractable\. The deep generative models we examine include likelihood\-based models, viz\., flow models, and score\-based models\. Experiments show the effectiveness of the proposed framework for estimating density functions for multivariate toy distributions as well as lattices from computational physics, namely, XY model andϕ4\\phi^\{4\}theory\. When applied to the renormalization problem inϕ4\\phi^\{4\}theory, the proposed framework estimates compact symbolic approximations of the hamiltonian function at different scales directly from samples, yielding expressions that may be challenging to derive using traditional perturbative or analytic approaches in nonperturbative settings\.

## 1Introduction

In the physical sciences, understanding complex phenomena fundamentally depends on modelling the underlying system dynamics\. The probability distributions that characterize these dynamics, as Hamiltonian or energy in a Boltzmann distribution, are traditionally formulated by domain experts, often in an analytical form with a set of parameters\. These parameters can be estimated from the data using regression or data\-fitting techniques\([Shanahan et al\., 2018](https://arxiv.org/html/2610.10807#bib.bib16)\)\. While effective, this process is labor\-intensive and is heavily reliant on expert intuition\. Consequently, the problem of automatically inferring an analytical or symbolic representation of a probability density function directly from observed data presents a significant challenge\([Tohme et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib15)\)\. Addressing this problem lies at the intersection of machine learning and physics, offering the potential to automate model discovery and to enhance our understanding of complex physical systems\. This problem is analogous to the system identification problem, widely studied in mechanical and control systems\([Sahoo et al\., 2018](https://arxiv.org/html/2610.10807#bib.bib27)\), where the system dynamics is modeled in differential form \(using differential equations\) instead of an integral form \(using Hamiltonian formulation\)\.

Boltzmann distributions are abundantly used in physical sciences\([Akhound\-Sadegh et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib3);[Midgley et al\., 2023](https://arxiv.org/html/2610.10807#bib.bib28);[Hasenfratz and Niedermayer, 1994](https://arxiv.org/html/2610.10807#bib.bib2)\)\. They have the formp⁡\(𝐱\)∝exp⁡\(−H⁡\(𝐱\)\)p\(\\mathbf\{x\}\)\\propto\\exp\(\-H\(\\mathbf\{x\}\)\), whereH⁡\(𝐱\)H\(\\mathbf\{x\}\)is a scalar denoting energy of the system in state𝐱\\mathbf\{x\}\. Popular deep generative approaches to model these distributions include likelihood\-based approaches, such as normalizing flows\([Dinh et al\., 2017](https://arxiv.org/html/2610.10807#bib.bib12)\)and auto\-regressive models\([Uria et al\., 2016](https://arxiv.org/html/2610.10807#bib.bib6);[Van Den Oord et al\., 2016](https://arxiv.org/html/2610.10807#bib.bib7)\), and energy\-based approaches, such as score matching\([Hyvärinen, 2005](https://arxiv.org/html/2610.10807#bib.bib29);[Song and Ermon, 2019](https://arxiv.org/html/2610.10807#bib.bib32)\)\. These methods are able to generate samples from the target distribution while also providing information about the model likelihood i\.e\., log likelihood in case of likelihood based approaches and the gradient thereof in case of score matching approaches\. However, these models are black box estimators\. On the other hand, Symbolic regression \(SR\) methods are used to estimate the equations of the input\-output relations in a supervised way\. Popular approaches for SR include genetic algorithms\([Schmidt and Lipson, 2009](https://arxiv.org/html/2610.10807#bib.bib14);[La Cava et al\., 2016](https://arxiv.org/html/2610.10807#bib.bib18);[Pal and Wang, 2017](https://arxiv.org/html/2610.10807#bib.bib1);[La Cava et al\., 2019](https://arxiv.org/html/2610.10807#bib.bib19);[Cranmer, 2023](https://arxiv.org/html/2610.10807#bib.bib21)\), symbolic neural networks\([Martius and Lampert, 2017](https://arxiv.org/html/2610.10807#bib.bib26);[Sahoo et al\., 2018](https://arxiv.org/html/2610.10807#bib.bib27);[Kim et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib23)\)and large pre\-trained networks\([Biggio et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib25);[Kamienny et al\., 2022](https://arxiv.org/html/2610.10807#bib.bib24);[Vastl et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib10)\), which work well for a small number of variables \(typically up to 10\)\. There are a few approaches\([Cranmer et al\., 2020](https://arxiv.org/html/2610.10807#bib.bib22)\)to combine SR with deep learning methods, which scale well with the number of variables, to simplify the SR task\. While a few approaches\([Tohme et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib15)\), attempt to recover expressions for low\-dimensional unnormalized distributions, they rely on assuming the distribution’s form\.

In this paper, we propose a framework to integrate deep generative modelling with SR to estimate the density functions as symbolic expressions\. We consider flow\-based and score\-based models for density estimation, and subsequently use SR methods, including genetic algorithms\([Cranmer, 2023](https://arxiv.org/html/2610.10807#bib.bib21)\)and EQL\-based approaches\([Kim et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib23)\), to estimate the symbolic density function while incorporating task\-specific inductive biases\. The effectiveness of SR methods deteriorates notably with an increase in the number of variables or the dimensionality of the data, often leading to inaccurately estimated expressions\([Cranmer et al\., 2020](https://arxiv.org/html/2610.10807#bib.bib22);[Biggio et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib25)\)\. For this, we employ a density factorization approach that leverages local dependencies within the distribution\. This strategy effectively reduces the number of variables, thereby enabling more efficient and accurate expression estimation\. Deep generative models sample from high density regions ofp⁡\(x\)p\(x\), but effective SR needs samples from both high and low density regions; therefore, we introduce noisy perturbation of samples\.

We apply this framework to general density functions, such as multi\-variate Gaussian and many\-well distributions, as well as to problems in computational Physics, such as XY model\([Kosterlitz and Thouless, 1973](https://arxiv.org/html/2610.10807#bib.bib11)\)and scalarϕ4\\phi^\{4\}theory\([Singha et al\., 2023](https://arxiv.org/html/2610.10807#bib.bib38);[Albergo et al\., 2019](https://arxiv.org/html/2610.10807#bib.bib37)\), and estimate the energy functions,H⁡\(𝐱\)H\(\\mathbf\{x\}\)\. We also estimate the effective Hamiltonian arising from the renormalization of the strongly coupledϕ4\\phi^\{4\}theory, for which an analytical derivation is challenging in the non\-perturbative setting\. These experiments highlight the broader applicability of the proposed framework\.

We emphasize that the primary contribution of this work is not the development of a new symbolic regression \(SR\) algorithm\. Instead, we propose a framework that integrates deep generative modeling with symbolic regression to estimate interpretable energy functions or Hamiltonians from samples of unnormalized distributions, incorporating physically motivated inductive biases\. The central idea is to convert an unsupervised density estimation problem into a supervised symbolic regression problem by leveraging quantities obtained from learned flow\-based and score\-based models, such as log\-likelihoods and score fields\. Standard SR methods can then be directly applied as modular components to obtain analytical approximations of the underlying energy landscape\.

Within this framework, we introduce three key components that improve symbolic estimation in practice\. First, we employ a local factorization of the Hamiltonian, exploiting spatial locality in many\-body systems to reduce the effective dimensionality of the symbolic search space and enable scalable estimation in high\-dimensional settings\. Second, we incorporate task\-specific inductive biases through the choice of SR primitives and structural constraints, which encode prior knowledge about the underlying physical system\. Third, we introduce a multi\-scale perturbation strategy that augments samples from high\-density regions produced by generative models with controlled noise, thereby improving coverage of low\-density regions that are critical for stable symbolic estimation\.

Overall, the proposed framework is complementary to existing advances in symbolic regression and symbolic distillation, and is designed to be compatible with a wide range of SR backends, including genetic\-programming\-based methods such as PySR and neural symbolic approaches such as EQL\.

![Refer to caption](https://arxiv.org/html/2610.10807v1/figures/new_model_diagram.drawio.png)

Figure 1:Block diagram of the proposed framework\. \(A\) Flow\-based estimation\. \(B\) Score\-based estimation\. \(C\) Integration\-free score\-based estimation\.
## 2The Framework

Our framework consists of the following parts\. \(1\) Training a deep generative model using the given samples\. \(2\) Fitting an SR model using the given samples and the corresponding likelihood information \(log values or gradients\) provided by the deep generative model, and with density factorization conforming to the known nature of the system\. The symbolic expression thus estimated gives the energyH⁡\(𝐱\)H\(\\mathbf\{x\}\)of the system\. A block schematic of the framework is depicted in Fig\.[1](https://arxiv.org/html/2610.10807#S1.F1)\.

Lattice Systems:Apart from general distributions, such as multivariate Gaussian and many\-well distributions, this paper focuses on lattice systems common in statistical and particle physics\. Lattices are obtained by discretizing the continuous space on a regular grid\. Recent works\([Albergo et al\., 2019](https://arxiv.org/html/2610.10807#bib.bib37);[Singha et al\., 2023](https://arxiv.org/html/2610.10807#bib.bib38);[Kanaujia and Arora, 2026](https://arxiv.org/html/2610.10807#bib.bib20)\)predominantly use normalizing flows to model the lattice distributions because of the availability of exact likelihood from the flow model\. In addition,[Wang et al\. \(2023\)](https://arxiv.org/html/2610.10807#bib.bib43)have investigated diffusion\-based approaches for the same purpose\. In this work, we start with flow based models, and subsequently, introduce score\-based methods to model the lattice systems\.

Problem statement:Given the samples𝒟=\{𝐱i∈𝒳\}i=1N\\mathcal\{D\}=\\\{\\mathbf\{x\}\_\{i\}\\in\\mathcal\{X\}\\\}\_\{i=1\}^\{N\}from an unknown distributionp⁡\(𝐱\);𝐱∈ℝdp\(\\mathbf\{x\}\);\\mathbf\{x\}\\in\\mathbb\{R\}^\{d\}, estimate a symbolic expressionH^​\(𝐱\)∈ℋ\\hat\{H\}\(\\mathbf\{x\}\)\\in\\mathcal\{H\}such thatqθ​\(𝐱\)=1Z​e−H^​\(𝐱\)q\_\{\\theta\}\(\\mathbf\{x\}\)=\\frac\{1\}\{Z\}e^\{\-\\hat\{H\}\(\\mathbf\{x\}\)\}approximatesp⁡\(𝐱\)p\(\\mathbf\{x\}\)\.

Symbolic Regression:SR is a supervised learning approach that aims to estimate an explicit mathematical expression describing the relationship between input and output variables\. Given a dataset𝒟=\{\(𝐱i,yi\)\}i=1N\\mathcal\{D\}=\\\{\(\\mathbf\{x\}\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{N\}, where𝐱i∈ℝn\\mathbf\{x\}\_\{i\}\\in\\mathbb\{R\}^\{n\}denotes the input variables,yi∈ℝy\_\{i\}\\in\\mathbb\{R\}denotes the corresponding output, andNNis the number of observations, SR seeks to identify a functionf:ℝn→ℝf:\\mathbb\{R\}^\{n\}\\rightarrow\\mathbb\{R\}such that

y^i=f⁡\(𝐱i\),\\hat\{y\}\_\{i\}=f\(\\mathbf\{x\}\_\{i\}\),\(1\)where the functional form offfis a concise closed\-form mathematical expression\. Unlike conventional regression methods, which estimate the parameters of a predefined model, symbolic regression simultaneously searches for both the structure of the equation and its parameters\([Schmidt and Lipson, 2009](https://arxiv.org/html/2610.10807#bib.bib14);[Schmidt and Lipson, 2010](https://arxiv.org/html/2610.10807#bib.bib17)\)\. The search is typically performed over a predefined library of primitive operators and basis functions \(e\.g\.,\+\+,−\-,×\\times,÷\\div,sin\\sin,cos\\cos, polynomials\), from which candidate expressions are constructed\.

The objective of SR is to find an expression that accurately models the observed data while remaining interpretable\. Letℱ\\mathcal\{F\}denote the space of candidate symbolic expressions generated from the chosen primitive library\. Given a loss functionℒ\\mathcal\{L\}that quantifies the discrepancy between the observed outputs and model predictions, SR seeks the optimal expression

f∗=arg⁡minf∈ℱ⁡ℒ⁡\(f\)\.f^\{\*\}=\\arg\\min\_\{f\\in\\mathcal\{F\}\}\\mathcal\{L\}\\\!\\left\(f\\right\)\.\(2\)
A common choice of loss function\(ℒ\)\(\\mathcal\{L\}\)is the mean squared error \(MSE\),

ℒMSE=1N​∑i=1N\(yi−f⁡\(𝐱i\)\)2\.\\mathcal\{L\}\_\{\\text\{MSE\}\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\left\(y\_\{i\}\-f\(\\mathbf\{x\}\_\{i\}\)\\right\)^\{2\}\.\(3\)
In practice, additional criteria such as expression complexity, sparsity, or computational cost may also be incorporated into the objective to favor simpler and more interpretable expressions\. The resulting symbolic model provides both predictive accuracy and an explicit analytical form, making SR particularly attractive for scientific discovery and equation identification tasks\([Dong and Zhong, 2025](https://arxiv.org/html/2610.10807#bib.bib44)\)\.

We use two approaches for SR, one is genetic programming based algorithm implemented using PySR library\([Cranmer, 2023](https://arxiv.org/html/2610.10807#bib.bib21)\)and the other is deep neural network based equation learner \(EQL\) approach\([Martius and Lampert, 2017](https://arxiv.org/html/2610.10807#bib.bib26);[Sahoo et al\., 2018](https://arxiv.org/html/2610.10807#bib.bib27);[Kim et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib23)\)\. PySR uses evolutionary algorithms to optimise a tree data structure that represents a mathematical expression in order to best fit the data\([Koza, 1994](https://arxiv.org/html/2610.10807#bib.bib13)\)in terms of mean squared error\. In the search for the underlying structure, model error is usually balanced with model complexity to ensure brevity of the learned equation and avoids overfitting\. PySR uses multi\-population evolutionary genetic algorithm with multiple evolutions performed asynchronously\. The other approach, EQL uses a symbolic neural network trained end\-to\-end through backpropagation\. The model minimizes the mean squared error betweenyyandy^\\hat\{y\}\. To enforce sparsity, a smoothed variant ofL0\.5L\_\{0\.5\}regularizer is used\. Additional details for both SR approaches are provided in Appendix[B](https://arxiv.org/html/2610.10807#A2)\.

The quality of convergence is indicated by the value of mean squared error\. Many times, SR fails to converge and this may be indicated by a high value of mean squared error\. One may use this to select the preferred solution from multiple SR models\.

Flow models:Flows model a target densityp⁡\(𝐱\)p\(\\mathbf\{x\}\)by applying a sequence of simple and learnable bijective transformations on samples from a known standard distributionq⁡\(𝐳\);𝐳∈ℝdq\(\\mathbf\{z\}\);\\mathbf\{z\}\\in\\mathbb\{R\}^\{d\}\([Papamakarios et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib33);[Kobyzev et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib35);[Rezende and Mohamed, 2015](https://arxiv.org/html/2610.10807#bib.bib34)\)\. In general,q⁡\(𝐳\)q\(\\mathbf\{z\}\)is chosen to be standard Gaussian\. Mathematically, a sample𝐳∼q⁡\(𝐳\)\\mathbf\{z\}\\sim q\(\\mathbf\{z\}\)is transformed to𝐱\\mathbf\{x\}with a bijective mapTθT\_\{\\theta\}, parameterized byθ\\theta\. The density of𝐱=Tθ​\(𝐳\)\\mathbf\{x\}=T\_\{\\theta\}\(\\mathbf\{z\}\)is given by

pθ​\(𝐱\)=q⁡\(Tθ−1​\(𝐱\)\)​\|det\(JTθ−1​\(𝐱\)\)\|p\_\{\\theta\}\(\\mathbf\{x\}\)=q\(T\_\{\\theta\}^\{\-1\}\(\\mathbf\{x\}\)\)\|\\det\(J\_\{T\_\{\\theta\}^\{\-1\}\}\(\\mathbf\{x\}\)\)\|\(4\)
Here,JTθ−1​\(𝐱\)=∂Tθ−1​\(𝐱\)/∂𝐱J\_\{T\_\{\\theta\}^\{\-1\}\}\(\\mathbf\{x\}\)=\\partial\{T\_\{\\theta\}^\{\-1\}\(\\mathbf\{x\}\)\}/\\partial\\mathbf\{x\}is the Jacobian ofTθ−1T\_\{\\theta\}^\{\-1\}\. In order to model the unknown target distributionp⁡\(𝐱\)p\(\\mathbf\{x\}\), using the samples fromp⁡\(𝐱\)p\(\\mathbf\{x\}\),θ\\thetais learnt by minimizing the KL divergence between the target densityp⁡\(𝐱\)p\(\\mathbf\{x\}\)and the model densitypθ​\(𝐱\)p\_\{\\theta\}\(\\mathbf\{x\}\)using gradient descent\.

Once trained, the flow model can be used to obtainy^=−log⁡pθ​\(𝐱\)\\hat\{y\}=\-\\log p\_\{\\theta\}\(\\mathbf\{x\}\)for any𝐱\\mathbf\{x\}, to form\{\(𝐱i,y^i\)\}\\\{\(\\mathbf\{x\}\_\{i\},\\hat\{y\}\_\{i\}\)\\\}pairs to train the SR models, which converts the unsupervised density estimation problem to a supervised SR problem\. An SR model estimates an expressiony^=f⁡\(𝐱\)\\hat\{y\}=f\(\\mathbf\{x\}\), which is related to the energy function as

f​\(𝐱\)=H^​\(𝐱\)\+log⁡Z^f\(\\mathbf\{x\}\)=\\hat\{H\}\(\\mathbf\{x\}\)\+\\log\\hat\{Z\}\(5\)The expression forH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)is obtained by dropping the terms inf⁡\(𝐱\)f\(\\mathbf\{x\}\)which are independent of𝐱\\mathbf\{x\}\.

Score Matching:Score function for any distributionp⁡\(𝐱\)p\(\\mathbf\{x\}\)is defined ass⁡\(𝐱\)=∇𝐱​log​p​\(𝐱\)s\(\\mathbf\{x\}\)=\\nabla\_\{\\mathbf\{x\}\}\\log p\(\\mathbf\{x\}\)\. An advantage of usings⁡\(𝐱\)s\(\\mathbf\{x\}\)\(overlog⁡p⁡\(𝐱\)\\log p\(\\mathbf\{x\}\)\) is that it is independent of the normalizing constantZZ\. In score matching\([Song and Ermon, 2019](https://arxiv.org/html/2610.10807#bib.bib32);[Song et al\., 2020](https://arxiv.org/html/2610.10807#bib.bib31)\), a deep modelsθ​\(𝐱\)s\_\{\\theta\}\(\\mathbf\{x\}\), parameterized byθ\\theta, estimates the target scores⁡\(𝐱\)s\(\\mathbf\{x\}\)\. Training is carried out by minimizing the expected square distance,

J⁡\(θ\)=12​𝔼p⁡\(𝐱\)​‖sθ​\(𝐱\)−s⁡\(𝐱\)‖2J\(\\theta\)=\\frac\{1\}\{2\}\\mathbb\{E\}\_\{p\(\\mathbf\{x\}\)\}\\\|s\_\{\\theta\}\(\\mathbf\{x\}\)\-s\(\\mathbf\{x\}\)\\\|^\{2\}\(6\)In the absence ofs⁡\(𝐱\)s\(\\mathbf\{x\}\),JJcan be simplified under certain regularity conditions as detailed in[Hyvärinen \(2005\)](https://arxiv.org/html/2610.10807#bib.bib29)by

J⁡\(θ\)=𝔼p⁡\(𝐱\)​\[Tr⁡\(∇xsθ​\(𝐱\)\)\+12​‖sθ​\(𝐱\)‖2\]J\(\\theta\)=\\mathbb\{E\}\_\{p\(\\mathbf\{x\}\)\}\[\\Tr\(\\nabla\_\{x\}s\_\{\\theta\}\(\\mathbf\{x\}\)\)\+\\frac\{1\}\{2\}\\\|s\_\{\\theta\}\(\\mathbf\{x\}\)\\\|^\{2\}\]\(7\)whereTr\\Tris the trace\. Once trained, the model can estimate the scoresθ​\(𝐱\)=∇𝐱​log​pθ​\(𝐱\)s\_\{\\theta\}\(\\mathbf\{x\}\)=\\nabla\_\{\\mathbf\{x\}\}\\log p\_\{\\theta\}\(\\mathbf\{x\}\)for any𝐱\\mathbf\{x\}to form pairs\{\(𝐱i,y^i=−sθ​\(𝐱i\)\)\}\\\{\(\\mathbf\{x\}\_\{i\},\\hat\{y\}\_\{i\}=\-s\_\{\\theta\}\(\\mathbf\{x\}\_\{i\}\)\)\\\}to train the SR models\.

The SR model estimates the expressiony^=f⁡\(𝐱\)\\hat\{y\}=f\(\\mathbf\{x\}\), which is related to the energy function as

f​\(𝐱\)=∇𝐱H^​\(𝐱\)f\(\\mathbf\{x\}\)=\\nabla\_\{\\mathbf\{x\}\}\\hat\{H\}\(\\mathbf\{x\}\)\(8\)which is free from the normalizing constant\. In principle, estimatingH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)fromf⁡\(𝐱\)f\(\\mathbf\{x\}\)requires the learned vector field to be conservative \(i\.e\., curl\-free\), such that it can be expressed as the gradient of a scalar potential\. In our setting, the score model is trained to approximate the gradient of the log\-density, which is theoretically conservative for the target distributions considered, although approximation errors may lead to small deviations in practice\. Finally,H^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)is estimated by integratingf⁡\(𝐱\)f\(\\mathbf\{x\}\)with respect to𝐱\\mathbf\{x\}\.

Denoising Score Matching:Score matching does not scale well to high dimensional𝐱\\mathbf\{x\}due to the high cost involved in computingTr⁡\(∇xsθ​\(𝐱\)\)\\Tr\(\\nabla\_\{x\}s\_\{\\theta\}\(\\mathbf\{x\}\)\)\. To address this,[Vincent \(2011\)](https://arxiv.org/html/2610.10807#bib.bib30)proposed denoising score matching approach that perturbs𝐱\\mathbf\{x\}with a pre\-defined noise distributionqσ​\(𝐱~\|𝐱\)q\_\{\\sigma\}\(\\tilde\{\\mathbf\{x\}\}\|\\mathbf\{x\}\), where,𝐱~=𝐱\+ϵ\\tilde\{\\mathbf\{x\}\}=\\mathbf\{x\}\+\\mathbf\{\\epsilon\}andϵ∼qσ​\(𝐱~\|𝐱\)=𝒩⁡\(𝐱~\|𝐱,σ2​I\)\\epsilon\\sim q\_\{\\sigma\}\(\\tilde\{\\mathbf\{x\}\}\|\\mathbf\{x\}\)=\\mathcal\{N\}\(\\tilde\{\\mathbf\{x\}\}\|\\mathbf\{x\},\\sigma^\{2\}I\)\. The score of the perturbed data distributionqσ​\(𝐱~\)q\_\{\\sigma\}\(\\tilde\{\\mathbf\{x\}\}\)is learned by minimizing

J⁡\(θ,σ\)=12​𝔼qσ​\(𝐱,𝐱~\)​\[‖sθ​\(𝐱~\)−∇𝐱~​log​qσ​\(𝐱~\|𝐱\)‖22\]J\(\\theta;\\sigma\)=\\frac\{1\}\{2\}\\mathbb\{E\}\_\{q\_\{\\sigma\}\(\\mathbf\{x\},\\tilde\{\\mathbf\{x\}\}\)\}\[\\\|s\_\{\\theta\}\(\\tilde\{\\mathbf\{x\}\}\)\-\\nabla\_\{\\tilde\{\\mathbf\{x\}\}\}\\log q\_\{\\sigma\}\(\\tilde\{\\mathbf\{x\}\}\|\\mathbf\{x\}\)\\\|^\{2\}\_\{2\}\]\(9\)For smallσ\\sigma, we getsθ​\(𝐱\)≈∇𝐱​log​p​\(𝐱\)s\_\{\\theta\}\(\\mathbf\{x\}\)\\approx\\nabla\_\{\\mathbf\{x\}\}\\log p\(\\mathbf\{x\}\)\. Training is carried out across multiple noise levels, which generates samples from low\-density regions of thep⁡\(𝐱\)p\(\\mathbf\{x\}\), thereby improving score estimation\.

Once trained, the score modelsθ​\(𝐱\)s\_\{\\theta\}\(\\mathbf\{x\}\)can be used to estimateH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)in the same way as described above in the score matching approach\.

Mitigating Integration Issues:In score\-based models, SR estimates the score function as a symbolic expression, which is subsequently integrated to estimateH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)\(see Eq\. \([8](https://arxiv.org/html/2610.10807#S2.E8)\)\)\. However, when SR yields a complex expression for score, its integration becomes analytically challenging\. To ameliorate this issue, we propose a reformulation of the modelling approach: instead of modelling the score function, we modelHθ​\(𝐱\)H\_\{\\theta\}\(\\mathbf\{x\}\)using a neural network\. The score is then obtained as the negative gradient obtained by autograd,sθ​\(𝐱\)=−∇𝐱Hθ​\(𝐱\)s\_\{\\theta\}\(\\mathbf\{x\}\)=\-\\nabla\_\{\\mathbf\{x\}\}H\_\{\\theta\}\(\\mathbf\{x\}\)\. This formulation enables training via standard objectives such as score matching or denoising score matching, as defined in Eqs\. \([7](https://arxiv.org/html/2610.10807#S2.E7)\) and \([9](https://arxiv.org/html/2610.10807#S2.E9)\)\.

Factorization:The effectiveness of SR methods reduces drastically with the increase in the size of search space \(the number of variables and operations\)\. The lattice models we apply our framework to are large\. However, most physical lattices have local interactions which can be leveraged to simplify the SR problem\. There are existing approaches leveraging this property for Monte Carlo sampling of lattices\([Kennedy and Pendleton, 1985](https://arxiv.org/html/2610.10807#bib.bib42);[Faraz et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib8)\)\. We use the same property for SR, however, allowing the local interactions to extend to any number of neighbours, makes our approach more generalizable\.

Given a lattice𝐱\\mathbf\{x\}\(a matrix, say\), the true energy function can be factorized as

H⁡\(𝐱\)=∑lh⁡\(xl\|xn⁡\(l\)\)\\displaystyle H\(\\mathbf\{x\}\)=\\sum\_\{l\}h\(x\_\{l\}\|x\_\{n\(l\)\}\)\(10\)where,hhis the local energy function,xlx\_\{l\}is the value of the field on lattice sitell, andn⁡\(l\)n\(l\)is the set of neighbors of sitell\. Fig\.[2](https://arxiv.org/html/2610.10807#S2.F2)\(a\) illustrates this\.

During training with SR,𝐱\\mathbf\{x\}is factorized into patches andH^\\hat\{H\}is estimated as

H^​\(𝐱\)=∑lh^​\(xp⁡\(l\)\)\\displaystyle\\hat\{H\}\(\\mathbf\{x\}\)=\\sum\_\{l\}\\hat\{h\}\(x\_\{p\(l\)\}\)\(11\)wherep⁡\(l\)p\(l\)is the set of lattice sites in the patch anchored at sitell\. Overall, SR estimatesh^\\hat\{h\}over variablesxp⁡\(l\)x\_\{p\(l\)\}which are small in number\. Fig\.[2](https://arxiv.org/html/2610.10807#S2.F2)\(b\) shows patches of different sizes\.

![Refer to caption](https://arxiv.org/html/2610.10807v1/figures/lattice.drawio.png)

Figure 2:Lattice Factorization\. \(a\) TrueH⁡\(𝐱\)H\(\\mathbf\{x\}\)can be factorized with a2×22\\times 2patch\. \(b\) While inference with SR, different patch sizes provide different inductive biases\. \(c\) Illustration for redundancy in representation\.It is to be noted that there is redundancy in this representation\. Multiple representations ofh^​\(xpl\)\\hat\{h\}\(x\_\{p\_\{l\}\}\)may lead to the sameH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)\. For instance, with a3×33\\times 3patch shown in Fig\.[2](https://arxiv.org/html/2610.10807#S2.F2)\(c\), anyh^​\(xp⁡\(l\)\)=∑i=08ai​g​\(xi\)\\hat\{h\}\(x\_\{p\(l\)\}\)=\\sum\_\{i=0\}^\{8\}a\_\{i\}g\(x\_\{i\}\), for any value of\{ai\}\\\{a\_\{i\}\\\}with∑i=08ai=a\\sum\_\{i=0\}^\{8\}a\_\{i\}=abeing constant, leads to the sameH^​\(𝐱\)=∑la​g​\(xl\)\\hat\{H\}\(\\mathbf\{x\}\)=\\sum\_\{l\}ag\(x\_\{l\}\)\. Here,ggcould be any univariate function\. In general,ggcould be multivariate; e\.g\.,h^′​\(xp⁡\(l\)\)=g⁡\(x0,x4,x6\)\\hat\{h\}^\{\\prime\}\(x\_\{p\(l\)\}\)=g\(x\_\{0\},x\_\{4\},x\_\{6\}\)andh^′′​\(xp⁡\(l\)\)=g⁡\(x1,x5,x7\)\\hat\{h\}^\{\\prime\\prime\}\(x\_\{p\(l\)\}\)=g\(x\_\{1\},x\_\{5\},x\_\{7\}\), with the same relative positions of lattice sites, lead to the sameH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)\.

In Fig\.[1](https://arxiv.org/html/2610.10807#S1.F1), methods \(A\) and \(C\), SR minimizes𝔼𝐱​\[‖y^​\(𝐱\)−∑lh^​\(xp⁡\(l\)\)‖2\]\\mathbb\{E\}\_\{\\mathbf\{x\}\}\[\\\|\\hat\{y\}\(\\mathbf\{x\}\)\-\\sum\_\{l\}\\hat\{h\}\(x\_\{p\(l\)\}\)\\\|^\{2\}\]to estimateh^\\hat\{h\}, wherey^​\(𝐱\)=−log⁡pθ​\(𝐱\)\\hat\{y\}\(\\mathbf\{x\}\)=\-\\log p\_\{\\theta\}\(\\mathbf\{x\}\)is the target estimated by the deep model\. It has no target for individualh^​\(xp⁡\(l\)\)\\hat\{h\}\(x\_\{p\(l\)\}\)but only matches the sum over the entire lattice\. On the other hand, in method \(B\), SR minimizes𝔼𝐱,l​\[‖sθ​\(𝐱\)​\[l\]−f⁡\(xp⁡\(l\)\)‖2\]\\mathbb\{E\}\_\{\\mathbf\{x\},l\}\[\\\|s\_\{\\theta\}\(\\mathbf\{x\}\)\[l\]\-f\(x\_\{p\(l\)\}\)\\\|^\{2\}\]to estimateff\(Eq\. equation[8](https://arxiv.org/html/2610.10807#S2.E8)\) with target as thelt​hl^\{th\}element of matrixsθ​\(𝐱\)s\_\{\\theta\}\(\\mathbf\{x\}\)estimated by the deep model\. Note that there is no redundancy of representation here as we are modelling site\-wise gradients and not the lattice likelihood\.

Learning from low density regions:Robust SR needs samples from both high and low values ofH⁡\(𝐱\)H\(\\mathbf\{x\}\)\. The original samples𝐱∼p⁡\(𝐱\)\\mathbf\{x\}\\sim p\(\\mathbf\{x\}\)represent regions with low energy, but high energy regions are conspicuous by their absence\. Hence, we append the original set𝒟\\mathcal\{D\}with noisy perturbations𝐱′=𝐱\+ϵ\\mathbf\{x\}^\{\\prime\}=\\mathbf\{x\}\+\\epsilonwhereϵ∼𝒩⁡\(𝟎,σ2​I\)\\mathbf\{\\epsilon\}\\sim\\mathcal\{N\}\(\\mathbf\{0\};\\sigma^\{2\}I\)\. Then, we estimatelog⁡pθ​\(𝐱′\)\\log p\_\{\\theta\}\(\\mathbf\{x\}^\{\\prime\}\)orsθ​\(𝐱′\)s\_\{\\theta\}\(\\mathbf\{x\}^\{\\prime\}\)using the deep models, and thus, form the input\-output pairs needed for SR\.

Algorithm[1](https://arxiv.org/html/2610.10807#alg1)outlines the training procedure for symbolic estimation of unnormalized distributions, combining deep generative modelling with SR\. The method first learns a density or score model from data, augments low\-density regions, and subsequently estimates an analytical expression for the Hamiltonian via SR\.

Algorithm 1Symbolic Function Estimation for Unnormalized Distributions1:Input:Data

𝒟=\{𝐱i\}i=1N\\mathcal\{D\}=\\\{\\mathbf\{x\}\_\{i\}\\\}\_\{i=1\}^\{N\}from distribution

2:Stage 1: Training the Generative Model

3:Initialise: Parameters

θ\\thetaof estimator

MM\(NF–based or score\-based\)

4:while

ℒM​\(θ\)\\mathcal\{L\}\_\{M\}\(\\theta\)not convergeddo

5:Sample a mini\-batch from

𝒟\\mathcal\{D\}
6:Compute loss

ℒM​\(θ\)\\mathcal\{L\}\_\{M\}\(\\theta\)
7:Update parameters:

θ←θ−η​∇θℒM​\(θ\)\\theta\\leftarrow\\theta\-\\eta\\nabla\_\{\\theta\}\\mathcal\{L\}\_\{M\}\(\\theta\)
8:endwhile

9:Stage 2: Low\-density Sample Augmentation

10:Generate perturbed samples

𝐱′←𝐱\+ϵ,ϵ∼𝒩⁡\(𝟎,σ2​I\)\\mathbf\{x\}^\{\\prime\}\\leftarrow\\mathbf\{x\}\+\\bm\{\\epsilon\},\\quad\\bm\{\\epsilon\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\sigma^\{2\}I\)
11:Construct augmented dataset

𝒟aug=𝒟∪\{𝐱i′\}i=1N\\mathcal\{D\}\_\{\\text\{aug\}\}=\\mathcal\{D\}\\cup\\\{\\mathbf\{x\}^\{\\prime\}\_\{i\}\\\}\_\{i=1\}^\{N\}
12:Stage 3: Label Construction

13:Evaluate the trained estimator

MMto obtain predicted quantities

y^i∈\{−NLL​\(𝐱i\),−s^​\(𝐱i\)\}\\hat\{y\}\_\{i\}\\in\\\{\-\\text\{NLL\}\(\\mathbf\{x\}\_\{i\}\),\\;\-\\hat\{s\}\(\\mathbf\{x\}\_\{i\}\)\\\}
14:Form regression dataset

𝒟′=\{\(𝐱i,y^i\)\}\\mathcal\{D\}^\{\\prime\}=\\\{\(\\mathbf\{x\}\_\{i\},\\hat\{y\}\_\{i\}\)\\\}
15:Stage 4: Symbolic Regression

16:Initialise: Parameters of the SR model

17:iflocal factorisable distributionthen

18:Get local patches from

𝒟′\\mathcal\{D\}^\{\\prime\}⊳\\trianglerightUse local Model Factorization

19:endif

20:Optimise SR parameters on

𝒟′\\mathcal\{D\}^\{\\prime\}using MSE loss

21:Recover analytical expression

y^=f⁡\(𝐱\)\\hat\{y\}=f\(\\mathbf\{x\}\)
22:Estimate Hamiltonian

H^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)from

f⁡\(𝐱\)f\(\\mathbf\{x\}\)
23:Return:

H^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)

## 3Case Studies

In this section we describe five specific examples where we apply our framework\.

Multivariate Gaussian density\.

p⁡\(𝐱,μ,𝚺\)\\displaystyle p\(\\mathbf\{x\};\\mathbf\{\\mu\},\\mathbf\{\\Sigma\}\)=1\(2​π\)d/2​\|𝚺\|1/2​e−12​\(𝐱−μ\)T​𝚺−1​\(𝐱−μ\)\\displaystyle=\\frac\{1\}\{\(2\\pi\)^\{d/2\}\|\\mathbf\{\\Sigma\}\|^\{1/2\}\}e^\{\-\\frac\{1\}\{2\}\(\\mathbf\{x\}\-\\mathbf\{\\mu\}\)^\{T\}\\mathbf\{\\Sigma\}^\{\-1\}\(\\mathbf\{x\}\-\\mathbf\{\\mu\}\)\}\(12\)Fixingμ=\[−1\.5,1\.5\]T\\mathbf\{\\mu\}=\[\-1\.5,1\.5\]^\{T\}, we use two variants of𝚺\\mathbf\{\\Sigma\}with𝐱∈ℝ2\\mathbf\{x\}\\in\\mathbb\{R\}^\{2\}\. In the first variant, we use a diagonal𝚺\\mathbf\{\\Sigma\}asII, resulting inH⁡\(𝐱\)=0\.5​\(x1\+1\.5\)2\+0\.5​\(x2−1\.5\)2H\(\\mathbf\{x\}\)=0\.5\(x\_\{1\}\+1\.5\)^\{2\}\+0\.5\(x\_\{2\}\-1\.5\)^\{2\}\. In the second one, we use a full𝚺\\mathbf\{\\Sigma\}as\[1\.00\.50\.51\.0\]\\begin\{bmatrix\}1\.0&0\.5\\\\ 0\.5&1\.0\\end\{bmatrix\}, leading toH⁡\(𝐱\)=0\.5​\(x1\+1\.5\)2\+0\.5​\(x2−1\.5\)2\+0\.5​\(x1\+1\.5\)​\(x2−1\.5\)H\(\\mathbf\{x\}\)=0\.5\(x\_\{1\}\+1\.5\)^\{2\}\+0\.5\(x\_\{2\}\-1\.5\)^\{2\}\+0\.5\(x\_\{1\}\+1\.5\)\(x\_\{2\}\-1\.5\)\. We apply our framework to the data sampled from these distributions, and expect it to recoverH⁡\(𝐱\)H\(\\mathbf\{x\}\)\.

Many Well Distribution \(dd\-dimensional\)\.It is a synthetic distribution given by the product ofd/2d/2copies of the 2\-D double well distribution\([Noé et al\., 2019](https://arxiv.org/html/2610.10807#bib.bib41);[Wu et al\., 2020](https://arxiv.org/html/2610.10807#bib.bib40)\)\. The Hamiltonian of 2\-D double well distribution is

H⁡\(𝐱\)=x14−6​x12−0\.5​x1\+0\.5​x22,𝐱∈ℝ2\\displaystyle H\(\\mathbf\{x\}\)=x\_\{1\}^\{4\}\-6x\_\{1\}^\{2\}\-0\.5x\_\{1\}\+0\.5x\_\{2\}^\{2\},\\qquad\\mathbf\{x\}\\in\\mathbb\{R\}^\{2\}\(13\)and that ofdd\-dimensional many well distribution is

H⁡\(𝐱\)\\displaystyle H\(\\mathbf\{x\}\)=∑i=1d/2\(x2​i−14−6​x2​i−12−0\.5​x2​i−1\+0\.5​x2​i2\),𝐱∈ℝd\\displaystyle=\\sum\_\{i=1\}^\{d/2\}\\left\(x\_\{2i\-1\}^\{4\}\-6x\_\{2i\-1\}^\{2\}\-0\.5x\_\{2i\-1\}\+0\.5x\_\{2i\}^\{2\}\\right\),\\qquad\\mathbf\{x\}\\in\\mathbb\{R\}^\{d\}\(14\)We conduct experiments ford=2,4,8,16,32,64d=2,4,8,16,32,64to evaluate the framework’s ability to estimateH⁡\(𝐱\)H\(\\mathbf\{x\}\)\.

XY Model\.It is a theoretical model used in statistical mechanics\([Kosterlitz and Thouless, 1973](https://arxiv.org/html/2610.10807#bib.bib11)\)to describe two\-dimensional systems of interacting spins and to study phase transitions\. It is widely employed to model Boltzmann distributions on lattices\([Kanaujia et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib39);[Singh et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib36);[Beach et al\., 2018](https://arxiv.org/html/2610.10807#bib.bib9)\)\. Here,𝐱∈\[0,2π\)d\\mathbf\{x\}\\in\[0,2\\pi\)^\{d\}represents the spin configuration on a square 2\-D grid withddsites, and the energy of the system is given by

H\(𝐱\)=−0\.5λ∑l=1d∑l′∈n⁡\(l\)cos\(xl−xl′\)H\(\\mathbf\{x\}\)=\-0\.5\\lambda\\sum\_\{l=1\}^\{d\}\\sum\_\{l^\{\\prime\}\\in n\(l\)\}\\cos\(x\_\{l\}\-x\_\{l^\{\\prime\}\}\)\(15\)Here,λ\\lambdais the coupling constant that quantifies interaction strength andn⁡\(l\)n\(l\)is a set of two neighbors \(right and below, cf\. Fig\.[2](https://arxiv.org/html/2610.10807#S2.F2)\(a\)\) of sitell\. We consider a grid of size8×88\\times 8, i\.e\.,d=64d=64, withλ=1/1\.4\\lambda=1/1\.4\. A patchppof size2×22\\times 2is used to factorizeH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)\.

Scalarϕ4\\phi^\{4\}theory\.It is a widely used in computational Physics\([Albergo et al\., 2019](https://arxiv.org/html/2610.10807#bib.bib37);[Singha et al\., 2023](https://arxiv.org/html/2610.10807#bib.bib38)\)to study scalar fields\. With𝐱∈ℝd\\mathbf\{x\}\\in\\mathbb\{R\}^\{d\}as a scalar field on a square 2\-D lattice withddsites, the energy function \(or action\) is given by

H⁡\(𝐱\)=∑l=1d\(λ1​xl4\+λ2​xl2\+2​∑l′∈n⁡\(l\)\(xl2−xl​xl′\)\)\\displaystyle H\(\\mathbf\{x\}\)=\\sum\_\{l=1\}^\{d\}\\left\(\\lambda\_\{1\}x\_\{l\}^\{4\}\+\\lambda\_\{2\}x\_\{l\}^\{2\}\+2\\sum\_\{l^\{\\prime\}\\in n\(l\)\}\(x\_\{l\}^\{2\}\-x\_\{l\}x\_\{l^\{\\prime\}\}\)\\right\)\(16\)wheren⁡\(l\)n\(l\)denotes the two neighboring sites \(right and below; see Fig\.[2](https://arxiv.org/html/2610.10807#S2.F2)\(a\)\)\. Here,λ1\\lambda\_\{1\}is the coupling constant andλ2\\lambda\_\{2\}represents the squared mass of the scalar particle associated with the field\. We conduct experiments on lattices of size8×88\\times 8, i\.e\.,d=64d=64, at two different values ofλ2\\lambda\_\{2\}, namely−4\-4and11, whileλ1\\lambda\_\{1\}is fixed at 4\. Here, a2×22\\times 2patchppis sufficient to recoverH⁡\(𝐱\)H\(\\mathbf\{x\}\)but we also experiment with patches of different sizes\. Larger patches allow recovering distant interactions but increase the number of variables, making it difficult for SR to estimate the true equation\.

Renormalization\.Renormalization in statistical Physics is akin to zooming out in image processing\. It enables one to study the phenomena on larger scales by reducing the resolution of lattices \(called coarse\-graining\) for computational efficiency\. Empirically, a finer lattice is coarsened by some kind of pooling; generally average pooling is used\. A fine lattice of sized×d\\sqrt\{d\}\\times\\sqrt\{d\}, on pooling with a2×22\\times 2filter, reduces to a coarser lattice of sized/2×d/2\\sqrt\{d\}/2\\times\\sqrt\{d\}/2\.

Traditionally, the effective actionHeff​\(𝐱\)H\_\{\\text\{eff\}\}\(\\mathbf\{x\}\)at the coarser scale is estimated analytically by integrating out terms in the original actionH⁡\(𝐱\)H\(\\mathbf\{x\}\)\. The effectiveness of this method is substantially limited by the extent of available analytical techniques\. Our framework, however, allows one to estimateHeff​\(𝐱\)H\_\{\\text\{eff\}\}\(\\mathbf\{x\}\)directly from coarse\-grained lattice samples\.

We test our framework to findHeff​\(𝐱\)H\_\{\\text\{eff\}\}\(\\mathbf\{x\}\)for highly correlatedϕ4\\phi^\{4\}theory withλ1=4\\lambda\_\{1\}=4andλ2=1\\lambda\_\{2\}=1, which lies in the non\-perturbative regime where perturbative analysis fails\. We start with64×6464\\times 64lattice samples from Eq\. \([16](https://arxiv.org/html/2610.10807#S3.E16)\), i\.e\.,d=4096d=4096, and coarse\-grain them to32×3232\\times 32lattice samples by average pooling\. These coarse\-grain samples are then used as𝒟\\mathcal\{D\}to estimateHeff​\(𝐱\)H\_\{\\text\{eff\}\}\(\\mathbf\{x\}\)\. These samples are progressively coarse\-grained to16×1616\\times 16, and then, to8×88\\times 8to see howHeff​\(𝐱\)H\_\{\\text\{eff\}\}\(\\mathbf\{x\}\)evolves to those lattice scales\.

## 4Experiments and Results

Most existing methods for symbolic density estimation assume a relatively simple parametric form for the underlying probability density function\. To the best of our knowledge, MeSSY\([Tohme et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib15)\)is one of the few existing methods specifically designed for symbolic density estimation that is applicable to our setting, and we therefore use it as the primary baseline\. MeSSY constructs an ansatz for the energy of the target distribution as a weighted linear combination of predefined basis functions, i\.e\., a set of linearly independent functions \(e\.g\.,xx,x2x^\{2\},sin⁡\(x\)\\sin\(x\),cos⁡\(x\)\\cos\(x\)\)\. The corresponding weights are estimated by maximizing the log\-likelihood over the available samples\. While we also investigated scaling MeSSY to higher\-dimensional problems, its reliance on numerical integration and rapidly growing polynomial basis functions resulted in prohibitive computational costs, limiting its applicability to very low\-dimensional settings in our experiments\. Consequently, MeSSY is reported only for cases where computation remained tractable\.

In contrast, the proposed approach integrates deep generative models with symbolic regression to estimate the analytic expression of density, using a rich set of operations \(e\.g\.,×,sin,cos,log,exp\\times,\\sin,\\cos,\\log,\\exp\) and relatively simpler basis functions\. This enables modelling complicated energy functions \(e\.g\., XY model\), and is not limited to only linear combinations of basis functions\.

In this section, we present experimental results across a range of datasets, spanning from low\-dimensional distributions \(d=2d=2\) to high\-dimensional ones \(d=1024d=1024\)\. These experiments validate the effectiveness of the proposed framework, highlighting its ability to handle increasing complexity and scale while maintaining performance\. We employ normalizing flows and score\-based models as the underlying generative components, which are coupled with two distinct SR frameworks, viz\., PySR \(genetic algorithm driven\)\([Cranmer, 2023](https://arxiv.org/html/2610.10807#bib.bib21)\)and EQL \(neural network based\)\([Kim et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib23)\), to form the integrated pipeline\.

The models are evaluated from two complementary perspectives: \(i\) predictive accuracy of the estimated Hamiltonian and \(ii\) fidelity of the estimated symbolic expression\. To assess predictive accuracy, we compare the estimated HamiltonianH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)with the ground\-truth HamiltonianH⁡\(𝐱\)H\(\\mathbf\{x\}\)using the mean squared error \(MSE\) and the coefficient of determination \(R2R^\{2\}\)\. TheR2R^\{2\}score is defined as

R2​\(H,H^,𝒟\)=1−∑𝐱∈𝒟\(H⁡\(𝐱\)−H^​\(𝐱\)\)2∑𝐱∈𝒟\(H⁡\(𝐱\)−H¯\)2,\\displaystyle R^\{2\}\(H,\\hat\{H\};\\mathcal\{D\}\)=1\-\\frac\{\\sum\_\{\\mathbf\{x\}\\in\\mathcal\{D\}\}\\left\(H\(\\mathbf\{x\}\)\-\\hat\{H\}\(\\mathbf\{x\}\)\\right\)^\{2\}\}\{\\sum\_\{\\mathbf\{x\}\\in\\mathcal\{D\}\}\\left\(H\(\\mathbf\{x\}\)\-\\bar\{H\}\\right\)^\{2\}\},\(17\)
whereH¯=𝔼𝐱∼𝒟​\[H​\(𝐱\)\]\\bar\{H\}=\\mathbb\{E\}\_\{\\mathbf\{x\}\\sim\\mathcal\{D\}\}\\left\[H\(\\mathbf\{x\}\)\\right\]denotes the empirical mean of the true Hamiltonian over the dataset𝒟\\mathcal\{D\}\. Lower MSE and higherR2R^\{2\}values indicate closer agreement betweenH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)andH⁡\(𝐱\)H\(\\mathbf\{x\}\)\. We additionally present the estimated symbolic expressions for qualitative comparison\.

While these metrics quantify predictive performance, the primary objective of the proposed framework is to estimate an interpretable symbolic representation of the underlying Hamiltonian\. As MSE andR2R^\{2\}alone do not fully characterize symbolic estimation quality, as distinct symbolic expressions may achieve similar predictive accuracy despite differing in structure, coefficient values, or the presence of spurious interactions\. To directly evaluate the quality of symbolic expressions, we report the following additional metrics\.

#### Relative Coefficient Error \(RCE\)\.

For symbolic terms that are correctly estimated, we measure coefficient accuracy using the average relative coefficient error\. Letℳ\\mathcal\{M\}denote the set of matched symbolic terms,cic\_\{i\}the ground\-truth coefficient, andc^i\\hat\{c\}\_\{i\}the estimated coefficient\. The RCE is defined as

RCE=1\|ℳ\|​∑i∈ℳ\|c^i−ci\|\|ci\|\.\\mathrm\{RCE\}=\\frac\{1\}\{\|\\mathcal\{M\}\|\}\\sum\_\{i\\in\\mathcal\{M\}\}\\frac\{\|\\hat\{c\}\_\{i\}\-c\_\{i\}\|\}\{\|c\_\{i\}\|\}\.\(18\)
An RCE of zero corresponds to exact coefficient estimation\.

#### Precision \(P\)\.

Precision measures the fraction of estimated symbolic terms that correspond to true terms in the target Hamiltonian,

Precision=TPTP\+FP,\\mathrm\{Precision\}=\\frac\{\\mathrm\{TP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FP\}\},\(19\)
where TP denotes the number of correctly estimated terms \(true positives\) and FP denotes the number of spurious terms \(false positives\)\.

#### Recall \(R\)\.

Recall measures the fraction of ground\-truth symbolic terms that are successfully estimated,

Recall=TPTP\+FN,\\mathrm\{Recall\}=\\frac\{\\mathrm\{TP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FN\}\},\(20\)
where FN denotes the number of ground\-truth terms that have not been estimated\.

#### F1 Score\.

The F1 score summarizes precision and recall through their harmonic mean,

F1=2​P​RP\+R\.\\mathrm\{F1\}=\\frac\{2PR\}\{P\+R\}\.\(21\)
An F1 score of one indicates perfect symbolic estimation\.

#### Number of Spurious Terms \(ST\)\.

We also report the number of symbolic terms present in the estimated expression that do not appear in the ground\-truth Hamiltonian\. This metric provides a direct measure of symbolic overfitting and highlights cases where good predictive performance is achieved through physically irrelevant interactions\.

Together, these metrics provide a comprehensive assessment of symbolic expression quality by quantifying predictive accuracy, coefficient estimation accuracy, structural correctness, and the presence of spurious interactions\.

We next describe the experimental setup and present the results separately for each dataset\.

Multivariate Gaussian density:We train both flow\-based and score\-matching models on2\.5×1042\.5\\times 10^\{4\}samples drawn from a two\-dimensional Gaussian distribution\. For symbolic regression \(SR\) with the EQL model, all2\.5×1042\.5\\times 10^\{4\}samples are used for training\. The proposed approach successfully estimates the true functional form ofH⁡\(𝐱\)H\(\\mathbf\{x\}\)\. Evaluation on an independent test set comprising2×1042\\times 10^\{4\}samples yields a mean squared error \(MSE\) close to zero and a coefficient of determination \(R2R^\{2\}\) close to one\. In comparison, PySR attains comparable performance using only10310^\{3\}samples\. These results hold consistently for both Gaussian distributions considered—those with diagonal and non\-diagonal covariance matrices𝚺\\mathbf\{\\Sigma\}\. Both methods estimate the trueH⁡\(𝐱\)H\(\\mathbf\{x\}\)without the need for noise perturbation to account for low\-density regions\. The quantitative results are presented in Table[1](https://arxiv.org/html/2610.10807#S4.T1)\.

Experiments with the baseline, MeSSY, also yield MSE values close to zero andR2R^\{2\}near one; however, the estimated expression deviates from the ground truth\. In particular, the coefficients of several polynomial terms differ from those in the trueH⁡\(𝐱\)H\(\\mathbf\{x\}\), and the estimated expression contains additional spurious multinomial terms\.

Table 1:Results and estimated expressionsh^​\(\[x1,x2\],λ\)\\hat\{h\}\(\[x\_\{1\},x\_\{2\}\];\\lambda\)for a 2\-D Gaussian distribution using various methods\. The mean vector is fixed atμ=\[−1\.5,1\.5\]T\\mu=\[\-1\.5,1\.5\]^\{T\}, with two different covariance matrices,Σ1\\Sigma\_\{1\}andΣ2\\Sigma\_\{2\}, as detailed in the table\. The target functions are: \(A\)0\.5​x12\+0\.5​x22\+1\.5​x1−1\.5​x20\.5x\_\{1\}^\{2\}\+0\.5x\_\{2\}^\{2\}\+1\.5x\_\{1\}\-1\.5x\_\{2\}, and \(B\)0\.67​x12−0\.67​x1​x2\+0\.67​x22\+3\.0​x1−3\.0​x20\.67x\_\{1\}^\{2\}\-0\.67x\_\{1\}x\_\{2\}\+0\.67x\_\{2\}^\{2\}\+3\.0x\_\{1\}\-3\.0x\_\{2\}\. Incorrect terms are shown in gray\. NF = Flow model, SM = Score matching\.Many Well Distribution \(dd\-dimensional\):We train flow, score\-matching, and denoising score\-matching models with10510^\{5\}samples from MW distribution\. For SR, we employ factorisation using patches of size2×12\\times 1to reduce the search space with this inductive bias\. Evaluation is performed on2×1042\\times 10^\{4\}test samples, and results are summarised in Table[2](https://arxiv.org/html/2610.10807#S4.T2)\.

For lower dimensions \(d=2,4,8d=2,4,8\), both flow and score\-matching models, as well as SR approaches \(PySR and EQL\), successfully estimate the trueH⁡\(𝐱\)H\(\\mathbf\{x\}\)\. EQL is trained using all10510^\{5\}samples, whereas PySR achieves comparable performance with only10310^\{3\}samples\. No noise perturbation is required\.

For higher dimensions \(d=16,32,64d=16,32,64\), the denoising score\-matching model outperforms the flow model, likely due to the highly multimodal nature ofp⁡\(𝐱\)p\(\\mathbf\{x\}\)with2d/22^\{d/2\}modes\. Flow models trained with forward KL divergence are susceptible to mode\-covering behavior\([Midgley et al\., 2023](https://arxiv.org/html/2610.10807#bib.bib28);[Kanaujia et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib39)\)\. In these cases, PySR also surpasses EQL, reflecting the advantage of the broader search capability of genetic algorithms over the gradient\-based optimization used in EQL\.

While MeSSY is able to estimate the general structure ofH⁡\(𝐱\)H\(\\mathbf\{x\}\)ford=2d=2, it fails to converge and estimate an expression ford\>2d\>2\. By “failure to converge,” we mean that the optimization procedure of MeSSY does not terminate successfully within a reasonable computational budget\. In our experiments, the algorithm repeatedly continued its search/optimization loop without producing a stable symbolic expression, even when trained using a very small number of samples for datasets with dimensions greater than two \(d\>2d\>2\)\. However, even ford=2d=2, the coefficients of several terms differ from those in the trueH⁡\(𝐱\)H\(\\mathbf\{x\}\), and the estimated expression includes additional multinomial terms \(cf\. Table[2](https://arxiv.org/html/2610.10807#S4.T2)\)\.

Table 2:Results for Many Well distribution with various dimensions,d=2,4,8,16,32,64d=2,4,8,16,32,64\. The true function ish⁡\(xpl\)=x14−6​x12−0\.5​x1\+0\.5​x22h\(x\_\{p\_\{l\}\}\)=x\_\{1\}^\{4\}\-6x\_\{1\}^\{2\}\-0\.5x\_\{1\}\+0\.5x\_\{2\}^\{2\}\. Estimated terms that do not match the true function are shown in gray\. NF = Flow model, SM = Score matching\.ddModelMSE\(↓\\downarrow\)R2\(↑\)R^\{2\}\(\\uparrow\)RCE\(↓\\downarrow\)ST\(↓\\downarrow\)P\(↑\\uparrow\)R\(↑\\uparrow\)F1\(↑\\uparrow\)Estimated Equationh^​\(xp⁡\(l\)\)\\hat\{h\}\(x\_\{p\(l\)\}\)2MeSSY0\.020\.980\.239620\.671\.000\.800\.96​x14\+0\.1​x13−5\.77​x12−0\.83​x1\+0\.61​x22−0\.23​x20\.96x\_\{1\}^\{4\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.1x\_\{1\}^\{3\}\}\-5\.77x\_\{1\}^\{2\}\-0\.83x\_\{1\}\+0\.61x\_\{2\}^\{2\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.23x\_\{2\}\}NF \+ PySR0\.000\.990\.000801\.001\.001\.00x14−6\.02​x12−0\.5​x1\+0\.5​x22x\_\{1\}^\{4\}\-6\.02x\_\{1\}^\{2\}\-0\.5x\_\{1\}\+0\.5x\_\{2\}^\{2\}NF \+ EQL0\.000\.990\.019601\.001\.001\.000\.99​x14−5\.95​x12−0\.49​x1\+0\.52​x220\.99x\_\{1\}^\{4\}\-5\.95x\_\{1\}^\{2\}\-0\.49x\_\{1\}\+0\.52x\_\{2\}^\{2\}SM \+ PySR0\.000\.990\.014201\.001\.001\.000\.98​x14−5\.9​x12−0\.51​x1\+0\.5​x220\.98x\_\{1\}^\{4\}\-5\.9x\_\{1\}^\{2\}\-0\.51x\_\{1\}\+0\.5x\_\{2\}^\{2\}SM \+ EQL0\.000\.990\.061701\.001\.001\.000\.99​x14−5\.9​x12−0\.59​x1\+0\.48​x220\.99x\_\{1\}^\{4\}\-5\.9x\_\{1\}^\{2\}\-0\.59x\_\{1\}\+0\.48x\_\{2\}^\{2\}4NF \+ PySR0\.000\.970\.029601\.001\.001\.001\.02​x14−6\.11​x12−0\.52​x1\+0\.52​x221\.02x\_\{1\}^\{4\}\-6\.11x\_\{1\}^\{2\}\-0\.52x\_\{1\}\+0\.52x\_\{2\}^\{2\}NF \+ EQL0\.000\.990\.046301\.001\.001\.001\.0​x14−6\.03​x12−0\.47​x1\+0\.46​x221\.0x\_\{1\}^\{4\}\-6\.03x\_\{1\}^\{2\}\-0\.47x\_\{1\}\+0\.46x\_\{2\}^\{2\}SM \+ PySR0\.010\.890\.018801\.001\.001\.001\.04​x14−6\.21​x12−0\.5​x1\+0\.5​x221\.04x\_\{1\}^\{4\}\-6\.21x\_\{1\}^\{2\}\-0\.5x\_\{1\}\+0\.5x\_\{2\}^\{2\}SM \+ EQL0\.480\.750\.096301\.001\.001\.001\.0​x14−6\.03​x12−0\.69​x1\+0\.50​x221\.0x\_\{1\}^\{4\}\-6\.03x\_\{1\}^\{2\}\-0\.69x\_\{1\}\+0\.50x\_\{2\}^\{2\}8NF \+ PySR0\.000\.900\.050801\.001\.001\.001\.0​x14−6\.02​x12−0\.55​x1\+0\.55​x221\.0x\_\{1\}^\{4\}\-6\.02x\_\{1\}^\{2\}\-0\.55x\_\{1\}\+0\.55x\_\{2\}^\{2\}NF \+ EQL0\.230\.960\.290040\.430\.750\.750\.92​x14−5\.58​x12−0\.86​x1\+0\.14​x12​x22\+0\.11​x12​x2\+0\.24​x1​x2−0\.2​x20\.92x\_\{1\}^\{4\}\-5\.58x\_\{1\}^\{2\}\-0\.86x\_\{1\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.14x\_\{1\}^\{2\}x\_\{2\}^\{2\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.11x\_\{1\}^\{2\}x\_\{2\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.24x\_\{1\}x\_\{2\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.2x\_\{2\}\}SM \+ PySR0\.010\.990\.021301\.001\.001\.00x14−6\.03​x12−0\.46​x1\+0\.5​x22x\_\{1\}^\{4\}\-6\.03x\_\{1\}^\{2\}\-0\.46x\_\{1\}\+0\.5x\_\{2\}^\{2\}SM \+ EQL0\.420\.930\.005010\.750\.750\.750\.99​x14−0\.11​x13−5\.97​x12\+0\.5​x220\.99x\_\{1\}^\{4\}\-0\.11x\_\{1\}^\{3\}\-5\.97x\_\{1\}^\{2\}\+0\.5x\_\{2\}^\{2\}16NF \+ PySR0\.510\.960\.130401\.001\.001\.001\.0​x14−6\.01​x12−0\.67​x1\+0\.59​x221\.0x\_\{1\}^\{4\}\-6\.01x\_\{1\}^\{2\}\-0\.67x\_\{1\}\+0\.59x\_\{2\}^\{2\}NF \+ EQL0\.670\.940\.190820\.671\.000\.800\.78​x14−4\.78​x12−0\.53​x1\+0\.64​x22\+0\.11​x12​x2−0\.26​x20\.78x\_\{1\}^\{4\}\-4\.78x\_\{1\}^\{2\}\-0\.53x\_\{1\}\+0\.64x\_\{2\}^\{2\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.11x\_\{1\}^\{2\}x\_\{2\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.26x\_\{2\}\}SM \+ PySR0\.650\.940\.143310\.801\.000\.891\.26​x14−5\.80​x12−0\.36​x1\+0\.50​x22−0\.06​x161\.26x\_\{1\}^\{4\}\-5\.80x\_\{1\}^\{2\}\-0\.36x\_\{1\}\+0\.50x\_\{2\}^\{2\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.06x\_\{1\}^\{6\}\}SM \+ EQL0\.060\.990\.047101\.001\.001\.00x14−5\.95​x12−0\.41​x1\+0\.5​x22x\_\{1\}^\{4\}\-5\.95x\_\{1\}^\{2\}\-0\.41x\_\{1\}\+0\.5x\_\{2\}^\{2\}32NF \+ PySR4\.630\.800\.366301\.001\.001\.000\.71​x14−4\.47​x12−0\.20​x1\+0\.34​x220\.71x\_\{1\}^\{4\}\-4\.47x\_\{1\}^\{2\}\-0\.20x\_\{1\}\+0\.34x\_\{2\}^\{2\}NF \+ EQL5\.290\.770\.753950\.380\.750\.500\.35​x14−2\.93​x12\+1\.05​x22−0\.33​x13​x2−0\.33​x12​x22−0\.11​x12​x2\+0\.48​x1​x2\+0\.36​x20\.35x\_\{1\}^\{4\}\-2\.93x\_\{1\}^\{2\}\+1\.05x\_\{2\}^\{2\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.33x\_\{1\}^\{3\}x\_\{2\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.33x\_\{1\}^\{2\}x\_\{2\}^\{2\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.11x\_\{1\}^\{2\}x\_\{2\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.48x\_\{1\}x\_\{2\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.36x\_\{2\}\}SM \+ PySR0\.480\.980\.130410\.801\.000\.890\.79​x14−4\.73​x12−0\.45​x1\+0\.50​x22−0\.02​x20\.79x\_\{1\}^\{4\}\-4\.73x\_\{1\}^\{2\}\-0\.45x\_\{1\}\+0\.50x\_\{2\}^\{2\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.02x\_\{2\}\}SM \+ EQL0\.690\.990\.363430\.571\.000\.730\.87​x14−5\.63​x12−1\.13​x1\+0\.49​x22\+0\.11​x16−0\.18​x15\+0\.54​x130\.87x\_\{1\}^\{4\}\-5\.63x\_\{1\}^\{2\}\-1\.13x\_\{1\}\+0\.49x\_\{2\}^\{2\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.11x\_\{1\}^\{6\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.18x\_\{1\}^\{5\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.54x\_\{1\}^\{3\}\}64NF \+ PySR1\.670\.930\.746730\.571\.000\.730\.13​x14−2\.06​x12\+1\.16​x1\+0\.43​x22\+0\.04​x16\+0\.19​x15−1\.14​x130\.13x\_\{1\}^\{4\}\-2\.06x\_\{1\}^\{2\}\+1\.16x\_\{1\}\+0\.43x\_\{2\}^\{2\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.04x\_\{1\}^\{6\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.19x\_\{1\}^\{5\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}1\.14x\_\{1\}^\{3\}\}NF \+ EQL3\.180\.931\.037920\.671\.000\.80−0\.93​x14\+0\.13​x12−0\.21​x1\+0\.19​x22\+0\.12​x13​x2−0\.48​x1​x2\-0\.93x\_\{1\}^\{4\}\+0\.13x\_\{1\}^\{2\}\-0\.21x\_\{1\}\+0\.19x\_\{2\}^\{2\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.12x\_\{1\}^\{3\}x\_\{2\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.48x\_\{1\}x\_\{2\}\}SM \+ PySR0\.530\.990\.135830\.571\.000\.731\.37​x14−6\.56​x12−0\.47​x1\+0\.51​x22−0\.06​x16−0\.02​x12​x2\+0\.05​x21\.37x\_\{1\}^\{4\}\-6\.56x\_\{1\}^\{2\}\-0\.47x\_\{1\}\+0\.51x\_\{2\}^\{2\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.06x\_\{1\}^\{6\}\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.02x\_\{1\}^\{2\}x\_\{2\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.05x\_\{2\}\}SM \+ EQL0\.570\.990\.582930\.571\.000\.731\.03​x14−6\.01​x12−1\.64​x1\+0\.49​x22−0\.22​x15\+0\.64​x13\+0\.23​x1​x21\.03x\_\{1\}^\{4\}\-6\.01x\_\{1\}^\{2\}\-1\.64x\_\{1\}\+0\.49x\_\{2\}^\{2\}\-\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.22x\_\{1\}^\{5\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.64x\_\{1\}^\{3\}\}\+\{\\color\[rgb\]\{0\.5,0\.5,0\.5\}0\.23x\_\{1\}x\_\{2\}\}

XY Model:We experiment on lattices of size8×88\\times 8\(d=64d=64\) withλ=1/1\.4\\lambda=1/1\.4\. The approximate energy functionH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)is factorized using patches of size2×22\\times 2\(c\.f\. Fig\.[2](https://arxiv.org/html/2610.10807#S2.F2)\(a\)\)\. We use10510^\{5\}samples for training flow and score\-based models\. For score\-based models, we find that parameterizingHθ​\(𝐱\)H\_\{\\theta\}\(\\mathbf\{x\}\)with a neural network and obtaining∇𝐱Hθ​\(𝐱\)\\nabla\_\{\\mathbf\{x\}\}H\_\{\\theta\}\(\\mathbf\{x\}\)via automatic differentiation yields better performance than directly modellingsθ​\(𝐱\)s\_\{\\theta\}\(\\mathbf\{x\}\), as the latter approach suffers from integration issues due to nested trigonometric terms\.

Denoising score matching outperforms standard score matching, while overall the flow model achieves superior results compared to score\-based methods\. Among symbolic regression approaches, PySR consistently performs better than EQL\. The quantitative results are summarized in Table[3](https://arxiv.org/html/2610.10807#S4.T3)\. In contrast, MeSSY fails to converge and is unable to estimate the target expression\.

Table 3:Results for XY Model\. The true function ish⁡\(xp⁡\(l\)\)=−0\.71​cos⁡\(x0−x1\)−0\.71​cos⁡\(x0−x2\)h\(x\_\{p\(l\)\}\)=\-0\.71\\cos\(x\_\{0\}\-x\_\{1\}\)\-0\.71\\cos\(x\_\{0\}\-x\_\{2\}\)\. Estimated terms that do not match the true function are shown in gray\. NF = Flow model, DSM = Denoising score model\.Scalarϕ4\\phi^\{4\}theory:We use10510^\{5\}samples for training flow and score\-based models\. Here, we train the score\-based model with denoising score matching\. For SR, the distribution is factorised over patches: a3×33\\times 3patch is used for estimating∇𝐱H^​\(𝐱\)\\nabla\_\{\\mathbf\{x\}\}\\hat\{H\}\(\\mathbf\{x\}\)via the score model, and a2×22\\times 2patch for estimatingH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)via the flow model\. Results are summarised in Table[4](https://arxiv.org/html/2610.10807#S4.T4)\.

Both flow and score\-based models accurately estimate the target energy function, while PySR consistently outperforms EQL\. For SR, a2×22\\times 2patch provides the most effective inductive bias for theϕ4\\phi^\{4\}theory\. However, we also assess larger patch sizes, namely3×33\\times 3and4×44\\times 4, to evaluate the robustness of the method when the interaction range is not known a priori\. The corresponding results forλ1=4\\lambda\_\{1\}=4andλ2=1\\lambda\_\{2\}=1using the flow model are reported in Table[5](https://arxiv.org/html/2610.10807#S4.T5)\. Although larger patches substantially increase the dimensionality of the symbolic search space, both PySR and EQL continue to estimate the dominant interaction structure\.

Notably, PySR consistently estimates the exact local Hamiltonian across all patch sizes without introducing any spurious interaction terms\. Even when provided with larger than necessary neighborhoods, the estimated expression retains only the physically relevant nearest\-neighbor interactions, indicating that the method is able to suppress unnecessary long\-range dependencies rather than merely fitting expressions within a prescribed patch\. This behavior provides evidence that the framework can identify the underlying interaction structure from an enlarged hypothesis space\.

EQL also estimates the dominant Hamiltonian terms for all patch sizes, but introduces an increasing number of spurious interactions as the patch size grows\. For the2×22\\times 2,3×33\\times 3, and4×44\\times 4patches, EQL produces 4, 4, and 7 additional terms, respectively\. Many of these involve variables located farther from the central site \(e\.g\.,x5x\_\{5\},x8x\_\{8\}, andx9x\_\{9\}for the4×44\\times 4patch\), corresponding to interactions that are absent in the true nearest\-neighbor Hamiltonian\. This trend suggests that the larger symbolic search space makes EQL more susceptible to fitting irrelevant long\-range dependencies, although the coefficients of the physically meaningful terms remain close to their ground\-truth values\.

Overall, these experiments demonstrate that the proposed framework can distinguish relevant from irrelevant interactions even when provided with a larger\-than\-necessary local neighborhood\. In practical applications where the true interaction range may be unknown, larger patches can therefore be employed without substantially affecting the estimation accuracy, particularly when using PySR\. In contrast, MeSSY fails to converge and is unable to estimate the target expression\.

Table 4:Results forϕ4\\phi^\{4\}theory with variousλ1,λ2\\lambda\_\{1\},\\lambda\_\{2\}settings\. The true expressions are: \(A\) forλ1=4,λ2=−4\\lambda\_\{1\}=4,\\lambda\_\{2\}=\-4,h⁡\(xp⁡\(l\)\)=4​x04−2​x0​x1−2​x0​x2h\(x\_\{p\(l\)\}\)=4x\_\{0\}^\{4\}\-2x\_\{0\}x\_\{1\}\-2x\_\{0\}x\_\{2\}, and \(B\) forλ1=4,λ2=1\\lambda\_\{1\}=4,\\lambda\_\{2\}=1,h⁡\(xp⁡\(l\)\)=4​x04\+5​x02−2​x0​x1−2​x0​x2h\(x\_\{p\(l\)\}\)=4x\_\{0\}^\{4\}\+5x\_\{0\}^\{2\}\-2x\_\{0\}x\_\{1\}\-2x\_\{0\}x\_\{2\}\. Incorrect terms are highlighted in gray\. NF = Flow model, DSM = Denoising score model\.Table 5:Results for estimatingh^​\(xp⁡\(l\)\)\\hat\{h\}\(x\_\{p\(l\)\}\)forϕ4\\phi^\{4\}theory with different patch sizes\|p\|\|p\|for flow model\. Trueh⁡\(xp⁡\(l\)\)=4​x04\+5​x02−2​x0​x1−2​x0​x2h\(x\_\{p\(l\)\}\)=4x\_\{0\}^\{4\}\+5x\_\{0\}^\{2\}\-2x\_\{0\}x\_\{1\}\-2x\_\{0\}x\_\{2\}\. Incorrect terms are shown in gray\.Renormalization:Given the comparable performance of flow\-based and score\-based models for theϕ4\\phi^\{4\}theory, we employ only the flow model for the renormalisation study\. For experiments, we use10510^\{5\}samples of64×6464\\times 64lattice size\. These lattices are progressively coarsened to multiple scales, namely,d=32×32d=32\\times 32,16×1616\\times 16and8×88\\times 8by applying2×22\\times 2average pooling\. Separate flow models are trained at each lattice scale\. Subsequently, symbolic regression is applied using PySR and EQL, though we report only PySR results due to better convergence and lower mean squared error betweenHθ​\(𝐱\)H\_\{\\theta\}\(\\mathbf\{x\}\)\(from the flow model\) andH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)\(from SR\)\. A2×22\\times 2patch is used for lattice factorization\.

Since the exact ground\-truth expression for the effective Hamiltonian is unknown for the coarsened lattices at multiple scales, we evaluate the results qualitatively through trends in the estimated coefficients and quantitatively through observable statistics\.

As the lattice size decreases, perturbative theory predicts a reduction in the coefficient ofxl4x\_\{l\}^\{4\}and an increase in the coefficient ofxl2x\_\{l\}^\{2\}\. Our results exhibit the similar trend even in the non\-perturbative regime, as shown in Fig\.[3](https://arxiv.org/html/2610.10807#S4.F3)\. At smaller lattice sizes, the coefficient ofxl4x\_\{l\}^\{4\}becomes negligible and is dropped by the method, validating the consistency of the estimatedH^​\(𝐱\)\\hat\{H\}\(\\mathbf\{x\}\)\.

To further validate the estimated Hamiltonians beyond coefficient trends, we also perform a quantitative observable\-based evaluation\. For each lattice scale \(32×3232\\times 32,16×1616\\times 16, and8×88\\times 8\), we generate samples from the estimated symbolic Hamiltonians using Hamiltonian Monte Carlo \(HMC\) and compare them with the corresponding coarse\-grained data\. We evaluate three physical observables: magnetization \(MM\), absolute magnetization \(\|M\|\|M\|\), and magnetic susceptibility \(χ\\chi\)\. The agreement between the coarse\-grained and generated distributions is quantified using the percentage overlap between the corresponding histograms and the Earth Mover Distance \(EMD\) between the observable distributions\([Kanaujia et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib39)\)\. As shown in Table[6](https://arxiv.org/html/2610.10807#S4.T6), the histogram overlap exceeds98%98\\%for all observables and lattice scales, while the EMD remains very small, ranging from approximately10−410^\{\-4\}to10−310^\{\-3\}\. These results indicate that the estimated Hamiltonians consistently reproduce the observable distributions of the coarse\-grained systems\.

In addition, we compare the first and second order moments of the observables\. Table[7](https://arxiv.org/html/2610.10807#S4.T7)shows that the means and standard deviations computed from samples generated using the estimated Hamiltonians closely match those obtained from the coarse\-grained data across all lattice sizes\. For example, at lattice size32×3232\\times 32, the mean absolute magnetization is0\.00530\.0053for the coarse\-grained data and0\.00550\.0055for the generated data, with corresponding standard deviations of0\.00390\.0039and0\.00410\.0041, respectively\. Similar agreement is observed for magnetization and magnetic susceptibility at all renormalization scales\. Overall, these results provide complementary evidence that the estimated symbolic Hamiltonians not only reproduce the qualitative renormalization trends shown in Fig\.[3](https://arxiv.org/html/2610.10807#S4.F3), but also preserve key statistical and physical properties of the coarse\-grained system\.

![Refer to caption](https://arxiv.org/html/2610.10807v1/figures/rg_flow.jpg)

Figure 3:Renormalization results for non\-perturbative region,λ1=4\\lambda\_\{1\}=4andλ2=1\\lambda\_\{2\}=1\.Table 6:Quantitative comparison of physical observables computed from coarse\-grained data and samples generated from the estimated symbolic Hamiltonians at different renormalization scales\. We use 100,000 samples to compute the observables\.Table 7:Observable statistics of the coarse\-grained data and generated samples\. Mean and standard deviation are reported for magnetization, absolute magnetization, and magnetic susceptibility across different lattice sizes\.

## 5Discussion

Error propagation in two\-stage learning:The proposed framework includes learning at two stages, first the deep generative model, and then, the SR model\. Previous works as in[Cranmer et al\. \(2020\)](https://arxiv.org/html/2610.10807#bib.bib22)report that using SR alone performs much worse than using it in conjunction with deep models for large scale learning\. We find the same with our attempts to apply SR directly to maximise the likelihood\. Here, we investigate how errors propagate through these two stages in our framework\. Table[8](https://arxiv.org/html/2610.10807#S5.T8)compares the MSE after flow model and that after the SR stage for various datasets\. We use mean subtraction to eliminate the constant term because of normalizing constantZZ\. We see that after SR, with either PySR or EQL, the MSE is less as compared to the MSE just after the flow model\. A plausible explanation could be that SR, by its low complexity estimates, reduces the noise in the likelihood information obtained from the flow model\. These results further support the use of SR methods instead of only black box deep learning based models for generating samples\.

The performance of the proposed framework is primarily governed by the accuracy of the underlying generative model, the extent to which the available data covers the relevant state space, and the scalability of the symbolic regression \(SR\) algorithm\. Since the SR stage operates on densities or score fields estimated by a generative model, inaccuracies in the learned density are directly propagated to the estimated symbolic expression\. This effect becomes particularly pronounced for highly multimodal distributions, where normalizing flows exhibit mode\-covering behavior\([Kanaujia et al\., 2024](https://arxiv.org/html/2610.10807#bib.bib39);[Kanaujia and Arora, 2026](https://arxiv.org/html/2610.10807#bib.bib20)\), leading to biased density estimates and, consequently, larger coefficient errors and spurious symbolic terms in the estimated Hamiltonian\.

Adequate sampling coverage is equally important for accurate symbolic expression\. The framework assumes that the available samples adequately represent all regions that contribute significantly to the target distribution\. In systems containing many metastable states separated by energy barriers, obtaining representative samples from all relevant regions can be challenging\. As a result, both density estimation and symbolic expression may deteriorate when important regions of the state space are under\-sampled\. This phenomenon is evident in the MW\-64 benchmark, which combines high dimensionality with strong multimodality and contains approximately2322^\{32\}metastable modes\. The resulting difficulty in accurately learning the underlying distribution leads to reduced symbolic expression accuracy\.

A second limitation arises from the scalability of the SR stage\. Although generative models can often accommodate high\-dimensional input, the search space explored by SR grows rapidly with the number of variables\. This challenge becomes particularly important in systems with long\-range interactions, where larger patches are required to capture all relevant interaction terms\. Increasing the patch size substantially enlarges the symbolic search space, making the estimation more challenging and increasing the likelihood of coefficient inaccuracies and spurious terms\. Consequently, the overall scalability of the framework is often constrained by the SR component rather than by the generative modeling stage\.

An additional limitation of the framework is its reliance on domain knowledge in specifying the symbolic library and patch structure\. The proposed approach is not intended as a fully automated scientific discovery system; rather, it employs physically motivated inductive biases to make symbolic estimation tractable in high\-dimensional settings\. Consequently, the quality of the estimated expressions depends on the expressiveness of the chosen symbolic library and the suitability of the patch decomposition for the underlying system\. If important functional forms or interaction patterns are excluded from the search space, the estimated symbolic expressions may only provide approximate descriptions of the target Hamiltonian\. Accordingly, the estimated expressions should be interpreted as compact and interpretable symbolic approximations of effective Hamiltonians, rather than exact representations of the underlying physical system\.

The score\-based experiments on XY model dataset illustrate an additional challenge associated with non\-Euclidean state spaces\. The XY model evolves on the periodic manifold\(S1\)N\(S^\{1\}\)^\{N\}, whereas standard denoising score matching employs Gaussian perturbations and Euclidean score fields\. This mismatch does not fully respect the circular topology of the underlying state space and can lead to larger score estimation errors near periodic boundaries and in regions containing highly oscillatory trigonometric interactions\. Since the Hamiltonian is computed by integrating the estimated score field, these errors may accumulate and reduce the accuracy of symbolic expression\. Overall, the framework performs best when the generative model accurately represents the target distribution, the available samples provide sufficient coverage of the relevant state space, and the local Hamiltonian can be represented using a moderate number of variables\. In more challenging settings, the estimated expressions should be viewed as compact and interpretable symbolic approximations of effective Hamiltonians rather than exact symbolic form of the true underlying dynamics\.

We also conduct a comprehensive analysis of the robustness, repeatability, and overall reproducibility of the proposed framework\. These aspects are thoroughly examined through a series of experiments presented in the Appendix[F](https://arxiv.org/html/2610.10807#A6)\. In addition, we also investigate the impact of various hyperparameters on the performance of the method, with a detailed discussion provided in the Appendix[D](https://arxiv.org/html/2610.10807#A4)\.

Table 8:Comparison of MSE\(H,Hθ\)\(H,H\_\{\\theta\}\), whereHθH\_\{\\theta\}is computed by flow model, and MSE\(H,H^\)\(H,\\hat\{H\}\), whereH^\\hat\{H\}is estimated by SR methods \(PySR and EQL\) using flow model outputs\.
## 6Conclusion

We present a framework that combines deep generative models and symbolic regression to estimate the energy function or Hamiltonian of a system from observed samples using domain\-specific prior knowledge\. We demonstrate the use of both likelihood\-based and score\-based generative models in conjunction with symbolic regression while utilizing density factorization and sample perturbation\. Experiments support the effectiveness of the proposed approach in reliably estimating the expressions for the Hamiltonian of the XY model and theϕ4\\phi^\{4\}theory, as well as density functions for several popular toy distributions\. As an important contribution towards the challenges theorists face in discretizing continuous theories, we demonstrate the ability of the proposed framework to estimate compact symbolic approximations of effective hamiltonians from discrete samples of locally interacting systems\. Furthermore, in renormalization experiments where the exact coarse\-grained Hamiltonian is unknown and challenging to derive analytically in the considered non\-perturbative setting, the proposed method estimates symbolic expressions that remain consistent with the coarse\-grained dynamics across multiple scales, even in non\-perturbative regimes\. An important limitation of the current work stems from the inability of existing SR methods to deal with very large numbers of variables\. This hinders the approach in discovering long\-range or non\-local interaction terms in lattice models, as they may require very large patches for factorization\. More generally, the framework relies on physically motivated inductive biases through the choice of symbolic libraries and patch structures, and the quality of expression depends on the expressiveness of these choices\. Nevertheless, our experiments yield good results with patches of size up to4×44\\times 4\(16 variables\)\. Future work includes studying how the mode\-covering behavior of deep generative models affects the estimated symbolic expressions, and extending the proposed framework to non\-equilibrium systems, which evolve dynamically\.

## 7Impact Statement

The potential of this work extends to many areas in computational sciences, including Physics, Chemistry and Climate science\. A common concern for traditional scientists is the inability of machine learning based approaches to give insights into the phenomenon under study\. The proposed framework enables scientists to harness the power of machine learning for analyzing large datasets while simultaneously estimating interpretable symbolic equations for locally interacting systems using domain\-specific prior knowledge\. The paper showcases challenging systems whose governing equations can be challenging to derive analytically in the considered non\-perturbative setting, demonstrating that the proposed framework can nevertheless estimate their approximate symbolic forms from observed data\. We believe this work can be taken up for follow up studies in the following areas:

- •Condensed matter Physics to study phase transitions and renormalization theories
- •Nuclear Physics to derive interaction potentials from scattering data
- •Climate sciences to model extreme climatic conditions from weather data

## References

- Akhound\-Sadeghet al\.\(2024\)T\. Akhound\-Sadegh, J\. Rector\-Brooks, J\. Bose, S\. Mittal, P\. Lemos, C\. Liu, M\. Sendera, S\. Ravanbakhsh, G\. Gidel, Y\. Bengio,et al\.Iterated denoising energy matching for sampling from boltzmann densities\.InInternational Conference on Machine Learning,pp\. 760–786\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Albergoet al\.\(2019\)M\. S\. Albergo, G\. Kanwar, and P\.E\. ShanahanFlow\-based generative models for markov chain monte carlo in lattice field theory\.Physical Review D100\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p4.1),[§2](https://arxiv.org/html/2610.10807#S2.p2.1),[§3](https://arxiv.org/html/2610.10807#S3.p5.1)\.
- Beachet al\.\(2018\)M\. J\. S\. Beach, A\. Golubeva, and R\. G\. MelkoMachine learning vortices at the kosterlitz\-thouless transition\.Phys\. Rev\. B97,pp\. 045207\.Cited by:[§3](https://arxiv.org/html/2610.10807#S3.p4.1)\.
- Biggioet al\.\(2021\)L\. Biggio, T\. Bendinelli, A\. Neitz, A\. Lucchi, and G\. ParascandoloNeural symbolic regression that scales\.InInternational Conference on Machine Learning,pp\. 936–945\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§1](https://arxiv.org/html/2610.10807#S1.p3.1)\.
- Chib and Greenberg \(1995\)S\. Chib and E\. GreenbergUnderstanding the metropolis\-hastings algorithm\.The american statistician49,pp\. 327–335\.Cited by:[§A\.3](https://arxiv.org/html/2610.10807#A1.SS3.p1.1)\.
- Cranmeret al\.\(2020\)M\. Cranmer, A\. Sanchez Gonzalez, P\. Battaglia, R\. Xu, K\. Cranmer, D\. Spergel, and S\. HoDiscovering symbolic models from deep learning with inductive biases\.Advances in neural information processing systems33,pp\. 17429–17442\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§1](https://arxiv.org/html/2610.10807#S1.p3.1),[§5](https://arxiv.org/html/2610.10807#S5.p1.1)\.
- Cranmer \(2023\)M\. CranmerInterpretable machine learning for science with pysr and symbolicregression\. jl\.arXiv preprint arXiv:2305\.01582\.Cited by:[§B\.4](https://arxiv.org/html/2610.10807#A2.SS4.p1.1),[Appendix B](https://arxiv.org/html/2610.10807#A2.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§1](https://arxiv.org/html/2610.10807#S1.p3.1),[§2](https://arxiv.org/html/2610.10807#S2.p10.1),[§4](https://arxiv.org/html/2610.10807#S4.p3.1)\.
- Dinhet al\.\(2017\)L\. Dinh, J\. Sohl\-Dickstein, and S\. BengioDensity estimation using real NVP\.InInternational Conference on Learning Representations,Cited by:[§B\.1](https://arxiv.org/html/2610.10807#A2.SS1.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Dong and Zhong \(2025\)J\. Dong and J\. ZhongRecent advances in symbolic regression\.ACM Computing Survey57\(11\)\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p9.1)\.
- Farazet al\.\(2024\)A\. Faraz, A\. Singha, D\. Chakrabarti, S\. Nakajima, and V\. AroraImprovement of heatbath algorithm in lft using generative models\.Proc\. of Science, 41st Intl\. Symp\. on Lattice Field Theory, LATTICE 2024,pp\. 1–10\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p21.1)\.
- Hasenfratz and Niedermayer \(1994\)P\. Hasenfratz and F\. NiedermayerPerfect lattice action for asymptotically free theories\.Nuclear Physics B414,pp\. 785–814\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Hyvärinen \(2005\)A\. HyvärinenEstimation of non\-normalized statistical models by score matching\.Journal of Machine Learning Research6,pp\. 695–709\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§2](https://arxiv.org/html/2610.10807#S2.p16.2)\.
- Kamiennyet al\.\(2022\)P\. Kamienny, S\. d’Ascoli, G\. Lample, and F\. ChartonEnd\-to\-end symbolic regression with transformers\.Advances in Neural Information Processing Systems35,pp\. 10269–10281\.Cited by:[Appendix B](https://arxiv.org/html/2610.10807#A2.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Kanaujia and Arora \(2026\)V\. Kanaujia and V\. AroraScoreNF: score\-based normalizing flows for sampling unnormalized distributions\.InICASSP 2026 \- 2026 IEEE International Conference on Acoustics, Speech and Signal Processing \(ICASSP\),Vol\.,pp\. 4336–4340\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p2.1),[§5](https://arxiv.org/html/2610.10807#S5.p2.1)\.
- Kanaujiaet al\.\(2024\)V\. Kanaujia, M\. S\. Scheurer, and V\. AroraAdvNF: reducing mode collapse in conditional normalising flows using adversarial learning\.SciPost Physics16,pp\. 132–160\.Cited by:[§B\.1](https://arxiv.org/html/2610.10807#A2.SS1.p4.1),[§3](https://arxiv.org/html/2610.10807#S3.p4.1),[§4](https://arxiv.org/html/2610.10807#S4.SS0.SSS0.Px5.p20.1),[§4](https://arxiv.org/html/2610.10807#S4.SS0.SSS0.Px5.p8.1),[§5](https://arxiv.org/html/2610.10807#S5.p2.1)\.
- Kennedy and Pendleton \(1985\)A\. Kennedy and B\. PendletonImproved heatbath method for monte carlo calculations in lattice gauge theories\.Physics Letters B156,pp\. 393–399\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p21.1)\.
- Kimet al\.\(2021\)S\. Kim, P\. Y\. Lu, S\. Mukherjee, M\. Gilbert, L\. Jing, V\. Čeperić, and M\. SoljačićIntegration of neural network\-based symbolic regression in deep learning for scientific discovery\.IEEE Transactions on Neural Networks and Learning Systems32,pp\. 4166–4177\.Cited by:[§B\.3](https://arxiv.org/html/2610.10807#A2.SS3.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§1](https://arxiv.org/html/2610.10807#S1.p3.1),[§2](https://arxiv.org/html/2610.10807#S2.p10.1),[§4](https://arxiv.org/html/2610.10807#S4.p3.1)\.
- Kobyzevet al\.\(2021\)I\. Kobyzev, S\. J\.D\. Prince, and M\. A\. BrubakerNormalizing flows: an introduction and review of current methods\.IEEE Transactions on Pattern Analysis and Machine Intelligence43,pp\. 3964–3979\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p12.1)\.
- Kosterlitz and Thouless \(1973\)J\. M\. Kosterlitz and D\. J\. ThoulessOrdering, metastability and phase transitions in two\-dimensional systems\.Journal of Physics C: Solid State Physics6,pp\. 1181–1203\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p4.1),[§3](https://arxiv.org/html/2610.10807#S3.p4.1)\.
- Koza \(1994\)J\. R\. KozaGenetic programming as a means for programming computers by natural selection\.Statistics and computing4,pp\. 87–112\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p10.1)\.
- La Cavaet al\.\(2019\)W\. La Cava, T\. Helmuth, L\. Spector, and J\. H\. MooreA probabilistic and multi\-objective analysis of lexicase selection andε\\varepsilon\-lexicase selection\.Evolutionary Computation27,pp\. 377–402\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- La Cavaet al\.\(2016\)W\. La Cava, L\. Spector, and K\. DanaiEpsilon\-lexicase selection for regression\.InProceedings of the Genetic and Evolutionary Computation Conference 2016,pp\. 741–748\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Martius and Lampert \(2017\)G\. S\. Martius and C\. LampertExtrapolation and learning equations\.In5th International Conference on Learning Representations, ICLR 2017\-Workshop Track Proceedings,Cited by:[§B\.3](https://arxiv.org/html/2610.10807#A2.SS3.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§2](https://arxiv.org/html/2610.10807#S2.p10.1)\.
- Midgleyet al\.\(2023\)L\. I\. Midgley, V\. Stimper, G\. N\. Simm, B\. Schölkopf, and J\. M\. Hernández\-LobatoFlow annealed importance sampling bootstrap\.InThe Eleventh International Conference on Learning Representations,Cited by:[§A\.2](https://arxiv.org/html/2610.10807#A1.SS2.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§4](https://arxiv.org/html/2610.10807#S4.SS0.SSS0.Px5.p8.1)\.
- Neal \(2011\)R\. M\. NealMCMC using hamiltonian dynamics\.Handbook of markov chain monte carlo2,pp\. 2\.Cited by:[§A\.4](https://arxiv.org/html/2610.10807#A1.SS4.p1.1)\.
- Noéet al\.\(2019\)F\. Noé, S\. Olsson, J\. Köhler, and H\. WuBoltzmann generators: sampling equilibrium states of many\-body systems with deep learning\.Science365,pp\. 1147\.Cited by:[§3](https://arxiv.org/html/2610.10807#S3.p3.1)\.
- Pal and Wang \(2017\)S\. K\. Pal and P\. P\. WangGenetic algorithms for pattern recognition\.CRC Press\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Papamakarioset al\.\(2021\)G\. Papamakarios, E\. Nalisnick, D\. J\. Rezende, S\. Mohamed, and B\. LakshminarayananNormalizing flows for probabilistic modeling and inference\.Journal of Machine Learning Research22,pp\. 1–64\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p12.1)\.
- Rezende and Mohamed \(2015\)D\. Rezende and S\. MohamedVariational inference with normalizing flows\.InInternational conference on machine learning,pp\. 1530–1538\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p12.1)\.
- Sahooet al\.\(2018\)S\. Sahoo, C\. Lampert, and G\. MartiusLearning equations for extrapolation and control\.InInternational Conference on Machine Learning,pp\. 4442–4450\.Cited by:[§B\.3](https://arxiv.org/html/2610.10807#A2.SS3.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§2](https://arxiv.org/html/2610.10807#S2.p10.1)\.
- Schmidt and Lipson \(2009\)M\. Schmidt and H\. LipsonDistilling free\-form natural laws from experimental data\.science324\(5923\),pp\. 81–85\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§2](https://arxiv.org/html/2610.10807#S2.p4.2)\.
- Schmidt and Lipson \(2010\)M\. Schmidt and H\. LipsonSymbolic regression of implicit equations\.InGenetic Programming Theory and Practice VII,pp\. 73–85\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p4.2)\.
- Shanahanet al\.\(2018\)P\. E\. Shanahan, A\. Trewartha, and W\. DetmoldMachine learning action parameters in lattice quantum chromodynamics\.Phys\. Rev\. D97,pp\. 094506\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p1.1)\.
- Singhet al\.\(2021\)J\. Singh, M\. Scheurer, and V\. AroraConditional generative models for sampling and phase transition indication in spin systems\.SciPost Physics11,pp\. 043–071\.Cited by:[§3](https://arxiv.org/html/2610.10807#S3.p4.1)\.
- Singhaet al\.\(2023\)A\. Singha, D\. Chakrabarti, and V\. AroraConditional normalizing flow for markov chain monte carlo sampling in the critical region of lattice field theory\.Phys\. Rev\. D107,pp\. 014512\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p4.1),[§2](https://arxiv.org/html/2610.10807#S2.p2.1),[§3](https://arxiv.org/html/2610.10807#S3.p5.1)\.
- Song and Ermon \(2019\)Y\. Song and S\. ErmonGenerative modeling by estimating gradients of the data distribution\.Advances in neural information processing systems32\.Cited by:[§B\.2](https://arxiv.org/html/2610.10807#A2.SS2.p1.1),[§B\.2](https://arxiv.org/html/2610.10807#A2.SS2.p3.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§2](https://arxiv.org/html/2610.10807#S2.p16.1)\.
- Songet al\.\(2020\)Y\. Song, S\. Garg, J\. Shi, and S\. ErmonSliced score matching: a scalable approach to density and score estimation\.InUncertainty in Artificial Intelligence,pp\. 574–584\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p16.1)\.
- Tohmeet al\.\(2024\)T\. Tohme, M\. Sadr, K\. Youcef\-Toumi, and N\. G\. HadjiconstantinouMESSY estimation: maximum\-entropy based stochastic and symbolic density estimation\.Trans\. Mach\. Learn\. Res\.\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p1.1),[§1](https://arxiv.org/html/2610.10807#S1.p2.1),[§4](https://arxiv.org/html/2610.10807#S4.p1.1)\.
- Uriaet al\.\(2016\)B\. Uria, M\. Côté, K\. Gregor, I\. Murray, and H\. LarochelleNeural autoregressive distribution estimation\.Journal of Machine Learning Research17,pp\. 1–37\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Van Den Oordet al\.\(2016\)A\. Van Den Oord, N\. Kalchbrenner, and K\. KavukcuogluPixel recurrent neural networks\.InInternational conference on machine learning,pp\. 1747–1756\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Vastlet al\.\(2024\)M\. Vastl, J\. Kulhánek, J\. Kubalík, E\. Derner, and R\. BabuškaSymformer: end\-to\-end symbolic regression using transformer\-based architecture\.IEEE Access\.Cited by:[§1](https://arxiv.org/html/2610.10807#S1.p2.1)\.
- Vincent \(2011\)P\. VincentA connection between score matching and denoising autoencoders\.Neural Computation23,pp\. 1661–1674\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p18.1)\.
- Wanget al\.\(2023\)L\. Wang, G\. Aarts, and K\. ZhouGenerative diffusion models for lattice field theory\.arXiv preprint arXiv:2311\.03578\.Cited by:[§2](https://arxiv.org/html/2610.10807#S2.p2.1)\.
- Wuet al\.\(2020\)H\. Wu, J\. Köhler, and F\. NoéStochastic normalizing flows\.Advances in Neural Information Processing Systems33,pp\. 5933–5944\.Cited by:[§3](https://arxiv.org/html/2610.10807#S3.p3.1)\.

## Appendix AData Generation

In this section, we briefly describe the methods used to generate samples from the distributions considered in this study\.

### A\.1Multivariate Gaussian Distribution

We use the ‘torch\.distributions‘ library to generate samples for training and testing the various models\. With a Gaussian distribution having a mean of\[−1\.5,1\.5\]\[\-1\.5,1\.5\]and unit covariance, we generate 25,000 samples for training, 5,000 for validation, and 5,000 for testing the NF, Score, and SR models\. Additionally, another set of 20,000 samples is used for inference with the SR model to compute metrics\. Similarly, for a Gaussian distribution with mean\[−1\.5,1\.5\]\[\-1\.5,1\.5\]and full covariance matrix\[1\.00\.50\.51\.0\]\\begin\{bmatrix\}1\.0&0\.5\\\\ 0\.5&1\.0\\end\{bmatrix\}, the same number of samples are generated as detailed above for the earlier setting\.

### A\.2Many\-Well Distribution

We generate samples from the Many\-Well \(MW\) distribution using rejection sampling as detailed in[Midgley et al\. \(2023\)](https://arxiv.org/html/2610.10807#bib.bib28)\. We explore six variants of the Many\-Well distribution with varying dimensions \(2, 4, 8, 16, 32, and 64\)\. For each variant, we generate 100,000 samples for training, 10,000 for validation, and 10,000 for testing the models\. Furthermore, an independent set of 20,000 samples is used for inference with the SR model to compute metrics\.

### A\.3XY Model distribution

We use the Metropolis\-Hastings \(MH\) algorithm\([Chib and Greenberg, 1995](https://arxiv.org/html/2610.10807#bib.bib5)\)to generate square lattices of size8×88\\times 8\(64 dim\) for the XY model\. We setλ=1\.4−1\\lambda=1\.4^\{\-1\}and the coupling constant\(JJ\) is set to unity\. For proposal distribution, we use a uniform distributionq⁡\(z\)q\(z\)in the MH algorithm\. On account of burn\-in, around 30000 initial samples are discarded from the MCMC chain\. After every 320 MCMC steps, a configuration is added to the set to have a minimum correlation across samples in the training data\. We generate 100,000 samples for training, 10,000 for validation, and 10,000 for testing various models\.

### A\.4Scalarϕ4\\phi^\{4\}theory distribution

We use the HMC algorithm\([Neal, 2011](https://arxiv.org/html/2610.10807#bib.bib4)\)to generate samples from the scalarϕ4\\phi^\{4\}theory distribution\. Lattices of size8×88\\times 8\(64 dimensions\) are generated\. The step size is set to 0\.05, with 20 steps in the leapfrog integrator of the HMC algorithm\. The first 1,000 samples are discarded on account of thermalization\. We generateϕ4\\phi^\{4\}data for two different mass settingsλ2=\{−4,1\}\\lambda\_\{2\}=\\\{\-4,1\\\}, while keeping the coupling constantλ1\\lambda\_\{1\}fixed at 4\. For each setting, we use 100,000 samples for training, 10,000 for validation, and 10,000 samples for testing the various models\.

For renormalization experiments in scalarϕ4\\phi^\{4\}field theory, we generate64×6464\\times 64\(4096\-dim\) lattices using the HMC algorithm\. The simulations are performed in the non\-perturbative regime with a coupling constantλ1=4\.0\\lambda\_\{1\}=4\.0and mass parameterλ2=1\\lambda\_\{2\}=1\. The HMC integrator uses a step size of 0\.05 and 25 leapfrog steps per trajectory\. To ensure thermalization, the first 5000 samples are discarded\. From the remaining data, we use 100,000 samples for training and 25,000 samples each for testing and validation\.

## Appendix BModel and Training details

In this section, we provide an overview of the model architecture and hyperparameters used in our experiments\. We primarily focus on three models\. First, we describe the architecture employed in the flow modelling, which is used to compute the negative log\-likelihood of the samples\. Next, we briefly outline the score estimation models used in the score\-based approach\. Finally, we discuss the two SR modelling approaches—EQL\([Kamienny et al\., 2022](https://arxiv.org/html/2610.10807#bib.bib24)\)and PySR\([Cranmer, 2023](https://arxiv.org/html/2610.10807#bib.bib21)\)—along with their respective training hyperparameters\. We run all experiments on four NVIDIA GeForce RTX 3080 GPUs, each with 10 GB of VRAM\. Training time for the PySR model ranged from 10 minutes to 90 minutes, whereas the EQL model required significantly more time, varying between 2 to 10 hours depending on the dataset\. While the training time for the flow\-based and score\-based models ranged from 0\.5 to 8 hours, depending on the dataset\.

### B\.1Normalizing Flows

We implement the Flow model using the RNVP architecture as detailed in[Dinh et al\. \(2017\)](https://arxiv.org/html/2610.10807#bib.bib12)\. For the multivariate Gaussian distribution, we utilize 2 Affine coupling layers for the first variant with unit covariance matrix and 4 coupling layers for the second variant with full covariance matrix, as described in Appendix[A\.1](https://arxiv.org/html/2610.10807#A1.SS1)\. Each coupling layer is composed of two dense layers, each containing 256 neurons with ReLU activation\. The output consists of two layers: one for scaling and one for translation\. The scaling layer outputs two neurons with Tanh activation, while the translation layer outputs two neurons with linear activation\. The model is trained with a batch size of 512, using the Adam optimizer with an initial learning rate of1×10−51\\times 10^\{\-5\}\.

For the many\-well distribution, we use 10 affine coupling layers to model the distribution across all variants \(n=\{2,4,8,16,32,64\}n=\\\{2,4,8,16,32,64\\\}\)\. The coupling architecture and training details remain the same as described above\.

For the XY model dataset, we use 10 affine coupling layers to model the distribution\. Each coupling layer consists of 2 convolutional layers with 256 filters of kernel size3×33\\times 3and ReLU activation, using circular padding\. The output layer is made up of 2 components: a scaling layer with Tanh activation and a translation layer with linear activation\. The model is trained with a batch size of 128 using the Adam optimizer with an initial learning rate of1×10−51\\times 10^\{\-5\}\.

Since we use a Gaussian distribution as the base distribution for the NF model, which has support over𝐑\\mathbf\{R\}, the model learns better when the dataset also lies in Euclidean space\. However, the XY dataset exists on a circular manifoldΘ∈\[0,2π\)N×N\\Theta\\in\[0,2\\pi\)^\{N\\times N\}\. To address this, we transform the dataset into Euclidean space using a sigmoid transformation, as detailed in[Kanaujia et al\. \(2024\)](https://arxiv.org/html/2610.10807#bib.bib39)\.

For the Scalarϕ4\\phi^\{4\}theory distribution, we use 12 affine coupling layers to model the distribution for both variants\. The coupling architecture and training details are the same as those used for the XY dataset, as discussed above\.

For the renormalization experiments in scalarϕ4\\phi^\{4\}field theory, we generate lattices of size64×6464\\times 64\. These are progressively coarsened via average pooling with a2×22\\times 2filter, reducing the spatial resolution at each step\. After the first pooling operation, the lattices are reduced to32×3232\\times 32in size\. At this resolution, we model the field distribution using RNVP flow architecture comprising 16 affine coupling layers, each employing 512 convolutional filters\. Upon further coarsening, the lattices are downsampled to16×1616\\times 16and subsequently to8×88\\times 8\. At the16×1616\\times 16level, we use 12 affine coupling layers, and at8×88\\times 8, we use 10 affine coupling layers\. For both of these resolutions, each coupling layer uses 256 convolutional filters\. The coupling layer architecture remains consistent across scales, with changes only in the number of filters to accommodate the reduced spatial complexity\.

### B\.2Score Model

We estimate the score from the data using a neural network that models the score functionsθ​\(x\)s\_\{\\theta\}\(x\)as followed in[Song and Ermon \(2019\)](https://arxiv.org/html/2610.10807#bib.bib32)\. The model is trained via score matching by minimizing the loss function as described in Eq\. \([7](https://arxiv.org/html/2610.10807#S2.E7)\)\.

For the multivariate Gaussian distribution, we use a multilayer perceptron \(MLP\) with the architecture \[2, 256, 256, 2\] to model the score function\. For the many\-well distribution, we model the score via an MLP with the architecture \[2, 256, 256, 256, 2\] forn=2n=2\. Forn=\{4,8,16,32,64\}n=\\\{4,8,16,32,64\\\}, we use an MLP architecture with 4 hidden layers, each consisting of 512 neurons\.

For the XY model and scalarϕ4\\phi^\{4\}distributions, we optimize the modeled score function using denoising score matching loss, as described in Eq\. \([9](https://arxiv.org/html/2610.10807#S2.E9)\)\. The score is modeled with the architecture provided by[Song and Ermon \(2019\)](https://arxiv.org/html/2610.10807#bib.bib32)\. The model, initially downloaded from their GitHub repository, has been adapted to fit the dataset configuration\. Training is carried out using the Adam optimizer, beginning with an initial learning rate of1×10−41\\times 10^\{\-4\}, which is subsequently reduced to1×10−51\\times 10^\{\-5\}\.

### B\.3EQL

This approach builds on the Equation\-based Learning \(EQL\) architecture\([Martius and Lampert, 2017](https://arxiv.org/html/2610.10807#bib.bib26);[Sahoo et al\., 2018](https://arxiv.org/html/2610.10807#bib.bib27);[Kim et al\., 2021](https://arxiv.org/html/2610.10807#bib.bib23)\), integrating it with other neural modules to allow end\-to\-end training via backpropagation\. To promote sparse and interpretable solutions, we employ a smoothed variant of theL0\.5L\_\{0\.5\}regularization, which prevents gradient singularities and improves convergence\. The regularization function is defined as:

L0\.5∗​\(w\)=\{\|w\|1/2,if​\|w\|≥a\(−w48​a3\+3​w24​a\+3​a8\)1/2,if​\|w\|<a\\displaystyle L^\{\*\}\_\{0\.5\}\(w\)=\\begin\{cases\}\|w\|^\{1/2\},&\\text\{if \}\|w\|\\geq a\\\\ \\left\(\-\\frac\{w^\{4\}\}\{8a^\{3\}\}\+\\frac\{3w^\{2\}\}\{4a\}\+\\frac\{3a\}\{8\}\\right\)^\{1/2\},&\\text\{if \}\|w\|<a\\end\{cases\}\(22\)
witha=0\.01a=0\.01in all experiments\. This penalty is selectively applied to the EQL\-specific weights when the model is embedded within a broader neural architecture\.

Training is conducted in two phases\. The first phase uses a higher learning rate of10−210^\{\-2\}, combined with a regularization weightλ=5×10−3\\lambda=5\\times 10^\{\-3\}for the multivariate gaussian and many\-well distributions and3×10−23\\times 10^\{\-2\}for theϕ4\\phi^\{4\}and XY model distributions\. This phase spans 20000 iterations, with small regularization weight promoting sparsity in the EQL layers\. At the end of this phase, EQL weights with magnitude less than a threshold \(α=0\.01\\alpha=0\.01\) are pruned and frozen\. In the second phase, training proceeds for 10000 iterations with a reduced learning rate of10−310^\{\-3\}and no regularization \(λ=0\\lambda=0\), enabling fine\-tuning of the pruned network\.

To balance model expressivity with interpretability, different sets of activation functions are used in the hidden and final layers\. For the many\-well distribution, we employ a functionally rich set of activations across all layers, including: Constant \(1 neuron\), Identity \(1 neuron\), Squareg2g^\{2\}\(1 neuron\), Power 4g4g^\{4\}\(1 neuron\), and Multiplicationg1⋅g2g\_\{1\}\\cdot g\_\{2\}\(1 neuron\)\. For the multivariate Gaussian distribution, a similar configuration to the many\-well distribution is used, with the exception of the Power 4 function, which is omitted\. In the case of theϕ4\\phi^\{4\}distribution, the Power 4 and Product functions are applied only in the hidden \(forward\) layers, while the final layer is limited to Constant, Identity, and Square functions to constrain the complexity of the output\. In the score\-based modelling ofϕ4\\phi^\{4\}, the Power 4 function is completely excluded from the primitive set to further reduce the search space and enhance tractability\.

For the XY model distribution, which exhibits a circular topology, we introduce an inductive bias by incorporating trigonometric primitives—sin\\sinandcos\\cos, each with 1 neuron—into both the hidden and final layers\. To prevent the nesting of trigonometric functions and reduce symbolic complexity, we adopt a feature transformation strategy in which the input kernel is explicitly expanded to include trigonometric variants\. For instance, for a kernel of size 2 with variables\[x1,x2,x3,x4\]\[x\_\{1\},x\_\{2\},x\_\{3\},x\_\{4\}\], the input is transformed to\[cos⁡\(x1\),…,cos⁡\(x4\),sin⁡\(x1\),…,sin⁡\(x4\)\]\[\\cos\(x\_\{1\}\),\\dots,\\cos\(x\_\{4\}\),\\sin\(x\_\{1\}\),\\dots,\\sin\(x\_\{4\}\)\]\. This approach is similarly extended to other kernel sizes\. This transformation not only helps suppress excessive equation complexity but also proves especially beneficial in score\-based modelling, as it simplifies integration tasks that would otherwise be intractable\.

### B\.4PySR

For the PySR experiments, we use the PySR\([Cranmer, 2023](https://arxiv.org/html/2610.10807#bib.bib21)\)library, which implements a multi\-population, tree\-based genetic algorithm where multiple evolutions are performed asynchronously\. We use PySR with a modified MSE loss, as described in Section[2](https://arxiv.org/html/2610.10807#S2), under factorization\. The experimental parameters are set as follows:

Sample Size:We train the model using varying sample sizes based on the convergence of the training loss\. For learning the equations with the flow\-based and score\-based approaches, either 512 or 1024 samples are used, depending on performance\. Experiments are conducted with sample sizes of 128, 256, 512, 1024, and 2048\. Beyond this range, we observe no further improvement in performance; hence, the aforementioned sample sizes are selected for the experiments\.

Maximum Complexity:This parameter controls the permissible size of symbolic equations, which corresponds to the number of nodes in the search tree\. It is set to 20 for the multivariate Gaussian and many\-well distributions, 30 for theϕ4\\phi^\{4\}distribution, and 40 for the XY model dataset\. For renormalisation experiments inϕ4\\phi^\{4\}theory, it is set to 30\. These values are chosen based on the convergence of the training loss\. Increasing the complexity further results in an exponential increase in search time\.

Population Size:The population size corresponds to the number of candidate solutions evolved in each generation\. It is set to 20 for the multivariate Gaussian and many\-well distributions and 30 for theϕ4\\phi^\{4\}and the XY model distributions, guided by training loss convergence\.

## Appendix CInductive Biases and Symbolic Search Space\.

The symbolic regression stage incorporates domain\-specific inductive biases through the choice of primitive operators, expression complexity, and patch size\. These design choices constrain the hypothesis space explored by the symbolic regressor and therefore directly influence the class of Hamiltonians that can be estimated\. To make these assumptions explicit, Tables[9](https://arxiv.org/html/2610.10807#A3.T9)and[10](https://arxiv.org/html/2610.10807#A3.T10)summarize the symbolic regression configurations used across all benchmarks, including the operator library, model complexity, patch size, true interaction range, and whether the primitive set contains the exact functional terms appearing in the target Hamiltonian\.

For EQL, symbolic expressions are represented using neural networks whose standard activation functions are replaced with symbolic operators such as identity, multiplication, polynomial, trigonometric, and exponential functions\. Expression complexity is controlled through the network depth and the number of neurons assigned to each operator type\. The corresponding configurations are reported in Table[9](https://arxiv.org/html/2610.10807#A3.T9)\. For PySR, a genetic\-programming\-based symbolic regression method, the inductive bias is determined by the primitive operator set, maximum expression complexity, and evolutionary population size\. The configurations used in our experiments are summarized in Table[10](https://arxiv.org/html/2610.10807#A3.T10)\.

It is important to emphasize that these choices incorporate domain knowledge rather than representing unconstrained symbolic estimation\. For example, quartic operators are included for the many\-well andϕ4\\phi^\{4\}benchmarks because quartic energy terms are known to be physically relevant, multiplicative operators are required to represent interaction terms, and trigonometric operators are included for the XY model because its Hamiltonian is naturally expressed through cosine interactions\. Similarly, the patch sizes were selected to match or exceed the known interaction range, ensuring that all relevant local interactions are accessible to the symbolic regressor\. Consequently, the estimated expressions should be interpreted as symbolic form within a physically motivated hypothesis space rather than completely assumptions\-free equation estimation\.

Table 9:EQL configurations and inductive biases used across benchmarks\. Numbers in parentheses denote the number of neurons associated with each symbolic operator\. The final column indicates whether the primitive library contains the exact functional terms present in the target Hamiltonian\.Table 10:PySR configurations and inductive biases used across benchmarks\.
## Appendix DImpact of Hyperparameter Choices on Method Performance

In this section, we discuss the sensitivity of the proposed framework to key hyperparameters governing both the generative models and the symbolic regression components\. The hyperparameters of the generative model are selected by monitoring performance on a validation set\. Similarly, the hyperparameters of the symbolic\-regression \(SR\) methods are adjusted based on the validation error and the model complexity, both of which should remain small\. In practice, PySR exhibits higher sensitivity to the number of training samples and to the complexity parameter, whereas EQL is more sensitive to the learning rate\.

i\. Number of samples for SR learning:EQL requires substantially more training data, typically on the order of 10k–100k samples, than PySR\. In PySR, we typically use 500–2000 samples\. This is in line with the general observation that genetic algorithms require less training data as compared to neural networks\. This number, however, increases with the data dimensionality and with the number of modes in the underlying distribution\. For example, achieving low complexity and low MSE for MW\-64 requires more samples than forΦ4\\Phi^\{4\}\(d = 64\), owing to the larger number of modes in MW\-64\.

ii\. Number of samples for NF or score\-model learning:The performance of both NF and score\-based models improves as the training sample size increases\. To analyse this effect on the overall framework, we conduct a study on the flow\-based pipeline \(NF with PySR\), trained on theϕ4\\phi^\{4\}distribution \(d=64d=64\) using varying numbers of training samples\. The corresponding results are summarised in Table[11](https://arxiv.org/html/2610.10807#A4.T11)below\.

Table 11:Effect of sample size on model performance\.iii\. Complexity control in PySR and EQL:For PySR, the complexity is typically set between 10 and 30\. In EQL, the model complexity is controlled by the architecture of the network, specifically by adjusting the number of layers and neurons in the primitive set of operators\. Deeper networks lead to increased model complexity\.

In both PySR and EQL, the method struggles to generate an interpretable equation when model complexity is too high or when the dataset fails to sufficiently capture all modes\. In multimodal settings, generative models such as NFs often exhibit mode\-covering behaviour\. This can lead to inaccuracies in likelihood estimation, which, when propagated to the SR stage, may result in the identification of spurious or extraneous terms\.

In PySR, both the dataset size and sample quality play a critical role\. By sample quality, we refer to the extent to which the data adequately represent all relevant modes of the distribution\. The complexity hyperparameter is equally important: as model complexity increases, the inferred symbolic expressions tend to include additional spurious terms\. While increasing the dataset size can initially suppress these extraneous terms, beyond a certain point the complexity constraint becomes the dominant factor determining the number of such terms\.

Similarly, in EQL, increasing the dataset size helps suppress additional terms in the inferred equation only up to a certain point\. Beyond this threshold, the intrinsic complexity of the model, as well as errors propagated from the generative model \(e\.g\., NF or score\-based model\), begin to significantly influence the accuracy of the SR\. These factors can ultimately lead to the appearance of extra terms in the final symbolic expression\.

## Appendix EMetrics

We outline the steps for calculating the overall MSE error andR2R^\{2\}score for the proposed approaches\. For each samplexix\_\{i\}, the corresponding target value isH⁡\(xi\)H\(x\_\{i\}\), while the symbolic regressor provides the predicted valueH^​\(xi\)\\hat\{H\}\(x\_\{i\}\)\. In the NF\-based and score\-based formulation, the estimated Hamiltonian is defined only up to an additive constant\. Besides integrating the predicted score field introduces an arbitrary integration constant\. Although this constant does not alter the underlying energy landscape or the associated probability distribution, it can artificially affect the MSE andR2R^\{2\}metrics during evaluation\. Therefore, before computing the metrics, we remove this offset by subtracting the empirical mean from both the target and predicted Hamiltonians:

- •Step1: Compute target meanH¯=1N​∑i=1NH⁡\(xi\)\\bar\{H\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}H\(x\_\{i\}\)
- •Step2: Mean normalize target values : H~​\(xi\)=H⁡\(xi\)−H¯\\tilde\{H\}\(x\_\{i\}\)=H\(x\_\{i\}\)\-\\bar\{H\}\(23\)
- •Step3: Compute predicted meanH^¯=1N​∑i=1NH^​\(xi\)\\bar\{\\hat\{H\}\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\hat\{H\}\(x\_\{i\}\)
- •Step4: Mean normalise predicted values: H^~​\(xi\)=H^​\(xi\)−H^¯\\tilde\{\\hat\{H\}\}\(x\_\{i\}\)=\\hat\{H\}\(x\_\{i\}\)\-\\bar\{\\hat\{H\}\}\(24\)
- •Step5: Compute MSE using mean normalized target values and mean normalized predicted values\. M​S​E=1N​∑i=1N\(H~​\(xi\)−H^~​\(xi\)\)2MSE=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\(\\tilde\{H\}\(x\_\{i\}\)\-\\tilde\{\\hat\{H\}\}\(x\_\{i\}\)\)^\{2\}\(25\)

Apply the same procedure to calculate theR2R^\{2\}score using the mean normalized target values and the mean normalized predicted values\. When the score is provided as input to the model, the symbolic regressor \(SR\) predicts the expression for−∇x\(H⁡\(x\)\)\-\\nabla\_\{x\}\(H\(x\)\)\. Upon integrating this expression, a constant term is introduced\. To evaluate the model, remove this constant and use the predictedH^​\(x\)\\hat\{H\}\(x\)without the constant\. Then, proceed with the same evaluation steps as before\.

## Appendix FRobustness Analysis

To ensure that the reported results are not artifacts of individual training runs, we investigate the repeatability of each method by performing multiple trials with distinct random seeds\. The corresponding performance metrics are aggregated and presented as seed\-wise distributions\. This analysis highlights the degree of variability across runs and provides a more comprehensive understanding of the reliability of the reported results\.

### F\.1Phi\-4 Benchmark

For the scalarϕ4\\phi^\{4\}theory benchmark withd=64d=64\(see Table[4](https://arxiv.org/html/2610.10807#S4.T4)in the main paper\), we report results aggregated over multiple runs with different random seeds\. To assess the robustness and stability of symbolic recovery, Tables[12](https://arxiv.org/html/2610.10807#A6.T12)and[13](https://arxiv.org/html/2610.10807#A6.T13)provide the symbolic recovery metrics and estimated expressions obtained for individual seed settings using the NF \+ PySR and NF \+ EQL approaches, respectively\. In addition, we also provide histogram\-based results aggregated over multiple runs w\.r\.t\. MSE andR2R^\{2\}metric to illustrate the overall variability and distribution\.

NF \+ PySR \(20 seeds\)

NF \+ EQL \(10 seeds\)

### F\.2Many\-Well Benchmark

For the Many\-Well distribution \(n=64n=64\) \(see Table[2](https://arxiv.org/html/2610.10807#S4.T2)in the main paper\), we report score\-based results aggregated over multiple runs with different random seeds\. Tables[14](https://arxiv.org/html/2610.10807#A6.T14)and[15](https://arxiv.org/html/2610.10807#A6.T15)present the symbolic recovery metrics and estimated expressions obtained for individual seed settings using the SM \+ PySR and SM \+ EQL approaches, respectively\. In addition, we provide histogram\-based results aggregated across all runs to illustrate the overall variability and distribution of performance w\.r\.t\. MSE andR2R^\{2\}metric\.

SM \+ PySR \(18 seeds\)

SM \+ EQL \(8 seeds\)\.

Table 12:Symbolic estimation results for theϕ4\\phi^\{4\}theory across 20 different random seed settings using NF \+ PySR\. The ground\-truth local Hamiltonian forλ1=4\\lambda\_\{1\}=4andλ2=1\\lambda\_\{2\}=1ish⁡\(xp⁡\(l\)\)=4​x04\+5​x02−2​x0​x1−2​x0​x2h\(x\_\{p\(l\)\}\)=4x\_\{0\}^\{4\}\+5x\_\{0\}^\{2\}\-2x\_\{0\}x\_\{1\}\-2x\_\{0\}x\_\{2\}\. Incorrect terms are highlighted in gray\. NF denotes the flow\-based model\.Table 13:Symbolic estimation results for theϕ4\\phi^\{4\}theory across 5 different random seed settings using NF \+ EQL\. The true expressions ish⁡\(xp⁡\(l\)\)=4​x04\+5​x02−2​x0​x1−2​x0​x2h\(x\_\{p\(l\)\}\)=4x\_\{0\}^\{4\}\+5x\_\{0\}^\{2\}\-2x\_\{0\}x\_\{1\}\-2x\_\{0\}x\_\{2\}\. Incorrect terms are highlighted in gray\.Table 14:Symbolic estimation results for the MW\-64 benchmark across 18 random seed settings using SM \+ PySR\. The true local Hamiltonian ish⁡\(xpl\)=x14−6​x12−0\.5​x1\+0\.5​x22h\(x\_\{p\_\{l\}\}\)=x\_\{1\}^\{4\}\-6x\_\{1\}^\{2\}\-0\.5x\_\{1\}\+0\.5x\_\{2\}^\{2\}\. Incorrect terms are highlighted in gray\.Overall, these seed\-dependent histograms demonstrate that the proposed methods achieve reliable performance across runs, with a consistently strong success rate\.

Table 15:Symbolic estimation results for the MW\-64 benchmark across 8 random seed settings using SM \+ EQL\. The true local Hamiltonian ish⁡\(xpl\)=x14−6​x12−0\.5​x1\+0\.5​x22h\(x\_\{p\_\{l\}\}\)=x\_\{1\}^\{4\}\-6x\_\{1\}^\{2\}\-0\.5x\_\{1\}\+0\.5x\_\{2\}^\{2\}\. Incorrect terms are highlighted in gray\.

## Appendix GAdditional Experiments on Long\-range Interaction Systems

To assess the applicability of the proposed framework beyond nearest\-neighbor systems, we consider a class of physical systems with long\-range interactions\. We consider a synthetic distribution inspired by the scalarϕ4\\phi^\{4\}model, where the interaction terms involve next nearest\-neighbor couplings rather than nearest\-neighbor interactions\. For an8×88\\times 8lattice, we factorize the Hamiltonian into overlapping5×55\\times 5patches\. The corresponding local Hamiltonianh⁡\(xp⁡\(l\)\)h\(x\_\{p\(l\)\}\), incorporating long\-range interactions, can be expressed in terms of patch elements as

x0x\_\{0\}x1x\_\{1\}x2x\_\{2\}x3x\_\{3\}x4x\_\{4\}x5x\_\{5\}x6x\_\{6\}x7x\_\{7\}x8x\_\{8\}x9x\_\{9\}x10x\_\{10\}x11x\_\{11\}x12x\_\{12\}x13x\_\{13\}x14x\_\{14\}x15x\_\{15\}x16x\_\{16\}x17x\_\{17\}x18x\_\{18\}x19x\_\{19\}x20x\_\{20\}x21x\_\{21\}x22x\_\{22\}x23x\_\{23\}x24x\_\{24\}
The Hamiltonian is factorized as

H⁡\(x\)=∑lh⁡\(xp⁡\(l\)\),H\(x\)=\\sum\_\{l\}h\(x\_\{p\(l\)\}\),\(26\)
where

h⁡\(xp⁡\(l\)\)=4​x124\+5​x122−x2​x12−x10​x12−x12​x14−x12​x22\.h\(x\_\{p\(l\)\}\)=4x\_\{12\}^\{4\}\+5x\_\{12\}^\{2\}\-x\_\{2\}x\_\{12\}\-x\_\{10\}x\_\{12\}\-x\_\{12\}x\_\{14\}\-x\_\{12\}x\_\{22\}\.\(27\)
Here,p⁡\(l\)p\(l\)denotes the set of lattice sites contained in the patch anchored at lattice sitell\.

Table 16:Results for a synthetic long\-range interaction system\. The true expression ish⁡\(xp⁡\(l\)\)=4​x124\+5​x122−x2​x12−x10​x12−x12​x14−x12​x22h\(x\_\{p\(l\)\}\)=4x\_\{12\}^\{4\}\+5x\_\{12\}^\{2\}\-x\_\{2\}x\_\{12\}\-x\_\{10\}x\_\{12\}\-x\_\{12\}x\_\{14\}\-x\_\{12\}x\_\{22\}Incorrect terms are highlighted in gray\. NF = Flow model, DSM = Denoising score model\.The results, summarized in Table[16](https://arxiv.org/html/2610.10807#A7.T16), indicate that the proposed framework remains effective when the interaction range extends beyond nearest neighbours, provided that the patch size is chosen to encompass the full interaction neighborhood\. In particular, NF \+ PySR and DSM \+ PySR estimate the correct symbolic structure with perfect precision and recall and no spurious terms, demonstrating that the framework can successfully identify interpretable long\-range interaction terms from samples alone\. However, when EQL is used as the symbolic regressor in both the NF and DSM pipelines, a significant number of spurious terms are estimated, especially in the NF \+ EQL pipeline\. Although most of the true Hamiltonian terms are estimated, as reflected by the relatively high recall values, the presence of numerous incorrect terms leads to poor precision\.

In physical systems with long\-range interactions, the primary challenge arises not from the generative modeling stage, but from the symbolic regression stage\. In the proposed framework, the first step uses either a normalizing flow \(NF\) or a score\-based model to convert the unsupervised density estimation problem into a supervised learning problem by generating suitable input\-output pairs\. The subsequent symbolic regression \(SR\) step is then used to estimate an interpretable expression for the underlying Hamiltonian\.

The scalability of the overall pipeline is therefore largely determined by the scalability of the SR method\. Most existing symbolic regression algorithms perform well only for a relatively small number of input variables \(typically fewer than 10\)\. In systems with long\-range interactions, larger patches are required to ensure that all relevant interaction terms are captured within the local Hamiltonian\. For example, interactions extending beyond nearest neighbors may require patch sizes of3×33\\times 3,4×44\\times 4,5×55\\times 5, or larger\. As the patch size increases, the number of variables presented to the symbolic regressor grows rapidly, significantly enlarging the search space of possible expressions\. This generally leads to reduced symbolic expression accuracy, larger coefficient errors, and an increased number of spurious terms\.

The results on synthetic long\-range interaction systems show that the proposed framework can still successfully estimate the correct interaction structure\. However, the difficulty of symbolic estimation increases as the number of variables increases\. In particular, we observed that PySR, which is based on genetic programming, remains comparatively robust in these higher\-dimensional settings, whereas EQL, a neural\-network\-based symbolic regressor, tends to exhibit a more pronounced degradation in performance, often producing additional spurious terms and less accurate coefficient estimates\.

Overall, the proposed framework is applicable to systems with long\-range interactions, provided that the patch size is large enough to capture the relevant dependencies\.

## Appendix HSymbolic Recovery Under Primitive Library Misspecification

To investigate the effect of library misspecification, we first consider the scalarϕ4\\phi^\{4\}benchmark and remove the explicit polynomial basis functions\{x,x2,x3,x4\}\\\{x,x^\{2\},x^\{3\},x^\{4\}\\\}from the primitive library, retaining only the basic arithmetic operators\{\+,−,×\}\\\{\+,\-,\\times\\\}together with the trigonometric operators\{sin,cos\}\\\{\\sin,\\cos\\\}\.

We note, however, that this setting should not be regarded as a strictly misspecified experiment\. Although explicit polynomial basis functions are removed, the target Hamiltonian remains representable because polynomial terms can still be estimated through repeated applications of the multiplication operator\. We therefore view this experiment as a test of robustness when important basis functions are removed, forcing the symbolic regressor to estimate them indirectly, rather than as a true misspecification experiment\. More generally, misspecification in symbolic regression is often not absolute, since many functions can be reconstructed or approximated through compositions of other operators\. For example, polynomial terms can be generated via repeated multiplication, while trigonometric functions admit polynomial approximations through Taylor expansions\. Consequently, the distinction between a specified and misspecified library often depends on the complexity of the resulting representation, rather than the strict representability\.

Table[17](https://arxiv.org/html/2610.10807#A8.T17)reports the results obtained using the flow model followed by symbolic regression\. Despite the removal of explicit polynomial basis functions from the primitive library, PySR estimates a highly accurate approximation of the target Hamiltonian, achieving anR2R^\{2\}value of0\.990\.99and estimating the dominant interaction structure\. In contrast, EQL exhibits a noticeable degradation in the performance and introduces additional spurious terms involving higher\-order polynomial combinations and nested trigonometric functions\. These results indicate that symbolic estimation becomes substantially more challenging when the symbolic hypothesis space \(primitive set\) is not well aligned with the true functional form\. Even when the target Hamiltonian remains representable, expression quality depends not only on the primitive set but also on the ability of the symbolic regressor to efficiently estimate the target expression from the available primitives\.

Table 17:Symbolic estimation on the scalarϕ4\\phi^\{4\}benchmark under primitive library misspecification\. Explicit polynomial primitives are removed from the symbolic library, leaving only arithmetic and trigonometric operators\.### H\.1Stronger Misspecification: XY Model Without Trigonometric Operators

To construct a more stringent misspecification setting, we consider the XY model, whose ground\-truth local Hamiltonian is

h⁡\(xp⁡\(l\)\)=−0\.71​cos⁡\(x0−x1\)−0\.71​cos⁡\(x0−x2\)\.h\(x\_\{p\(l\)\}\)=\-0\.71\\cos\(x\_\{0\}\-x\_\{1\}\)\-0\.71\\cos\(x\_\{0\}\-x\_\{2\}\)\.\(28\)
Unlike theϕ4\\phi^\{4\}experiment above, we remove all trigonometric operators \(sin\\sinandcos\\cos\) from the primitive library and retain only arithmetic operators\. In this setting, the target interaction is no longer directly representable within the primitive set, and the symbolic regressor must approximate the cosine interaction indirectly through polynomial compositions\.

Using the Taylor expansion

cos⁡\(z\)=1−z22\!\+z44\!−⋯,\\cos\(z\)=1\-\\frac\{z^\{2\}\}\{2\!\}\+\\frac\{z^\{4\}\}\{4\!\}\-\\cdots,\(29\)
the XY Hamiltonian can be approximated by a polynomial containing quadratic, quartic, and mixed interaction terms\. Consequently, successful estimation in this setting requires the symbolic regressor to construct an indirect polynomial approximation to the underlying trigonometric interaction\.

h⁡\(xp⁡\(l\)\)\\displaystyle h\(x\_\{p\(l\)\}\)≈−0\.71​\(1−\(x0−x1\)22\!\+\(x0−x1\)44\!\+⋯\)−0\.71​\(1−\(x0−x2\)22\!\+\(x0−x2\)44\!\+⋯\)\\displaystyle\\approx\-0\.71\\left\(1\-\\frac\{\(x\_\{0\}\-x\_\{1\}\)^\{2\}\}\{2\!\}\+\\frac\{\(x\_\{0\}\-x\_\{1\}\)^\{4\}\}\{4\!\}\+\\cdots\\right\)\-0\.71\\left\(1\-\\frac\{\(x\_\{0\}\-x\_\{2\}\)^\{2\}\}\{2\!\}\+\\frac\{\(x\_\{0\}\-x\_\{2\}\)^\{4\}\}\{4\!\}\+\\cdots\\right\)\(30\)≈−1\.42\+0\.71​x02\+0\.36​x12\+0\.36​x22−0\.71​x0​x1−0\.71​x0​x2\\displaystyle\\approx\-1\.42\+0\.71x\_\{0\}^\{2\}\+0\.36x\_\{1\}^\{2\}\+0\.36x\_\{2\}^\{2\}\-0\.71x\_\{0\}x\_\{1\}\-0\.71x\_\{0\}x\_\{2\}\+0\.06​x04\+0\.03​x14\+0\.03​x24−0\.12​x03​x1−0\.12​x03​x2\\displaystyle\+0\.06x\_\{0\}^\{4\}\+0\.03x\_\{1\}^\{4\}\+0\.03x\_\{2\}^\{4\}\-0\.12x\_\{0\}^\{3\}x\_\{1\}\-0\.12x\_\{0\}^\{3\}x\_\{2\}\+0\.12​x0​x13−0\.12​x0​x23\+0\.18​x02​x12\+0\.18​x02​x22\+⋯\\displaystyle\+0\.12x\_\{0\}x\_\{1\}^\{3\}\-0\.12x\_\{0\}x\_\{2\}^\{3\}\+0\.18x\_\{0\}^\{2\}x\_\{1\}^\{2\}\+0\.18x\_\{0\}^\{2\}x\_\{2\}^\{2\}\+\\cdots
The results are reported in Table[18](https://arxiv.org/html/2610.10807#A8.T18)\. Despite the absence of trigonometric operators, both PySR and EQL recover meaningful approximations to the underlying interaction, achievingR2R^\{2\}values of 0\.81 and 0\.85, respectively\. The estimated expressions primarily consist of low\-order polynomial interaction terms, consistent with the expected Taylor\-series approximation of the cosine function\.

Table 18:Symbolic estimation results for the XY model under primitive library misspecification\. The symbolic regressor approximates the cosine interaction indirectly through polynomial compositions\.Compared to theϕ4\\phi^\{4\}experiment, this setting constitutes a substantially stronger form of misspecification because the true functional form is absent from the primitive library\. Nevertheless, the quantitative results in Table[18](https://arxiv.org/html/2610.10807#A8.T18)show that both symbolic regression methods are able to estimate meaningful approximations of the underlying cosine interaction\. Although the estimated expressions do not match the true Hamiltonian exactly, they capture several dominant low\-order terms expected from the Taylor expansion of the cosine function, including quadratic self\-interaction terms and pairwise interactions such asx0​x1x\_\{0\}x\_\{1\}andx0​x2x\_\{0\}x\_\{2\}\. The resulting approximations achieve reasonably high coefficients of determination \(R2=0\.81R^\{2\}=0\.81for PySR andR2=0\.85R^\{2\}=0\.85for EQL\), indicating that a substantial fraction of the Hamiltonian variance is explained even under strong library misspecification\. Interestingly, EQL attains a slightly higherR2R^\{2\}and lower MSE than PySR in this setting, suggesting that lower\-order polynomial approximations may provide a more robust surrogate when the true trigonometric functional form is unavailable\. Overall, these results indicate that the proposed framework can still estimate useful symbolic approximations when the target Hamiltonian is not directly representable, although with a noticeable degradation in accuracy relative to the well\-specified setting\. More broadly, they highlight that symbolic expression quality depends not only on primitive library but also on the ability of the symbolic regressor to construct effective approximations within the restricted hypothesis space\.

Taken together, theϕ4\\phi^\{4\}and XY experiments highlight an important aspect of the framework: successful symbolic estimation depends not only on the quality of the learned generative model but also on the expressiveness of the primitive library\. When the target functional form is not directly available within the primitive set, the symbolic regressor may still estimate useful approximations through compositions of simpler operators, albeit with increased complexity and reduced accuracy\. Consequently, the estimated expressions should be interpreted as effective symbolic approximations rather than exact recoveries of the underlying Hamiltonian\.

### H\.2Error Propagation Under Primitive Set Misspecification

A central claim of the proposed framework is that symbolic regression can act as a denoising or compression step by reducing the approximation error present in the learned generative model\. However, when important basis functions are removed from the primitive library, it is not obvious whether this error\-reduction property continues to hold\.

To investigate this question, we compare the error of the learned flow model,MSE⁡\(H,Hθ\)\\mathrm\{MSE\}\(H,H\_\{\\theta\}\), against the error of the estimated symbolic Hamiltonian,MSE⁡\(H,H^\)\\mathrm\{MSE\}\(H,\\hat\{H\}\), for the restricted\-library experiments described above\. The results are reported in Table[19](https://arxiv.org/html/2610.10807#A8.T19)\.

Table 19:Error propagation analysis under primitive library misspecification\.MSE⁡\(H,Hθ\)\\mathrm\{MSE\}\(H,H\_\{\\theta\}\)denotes the error of the learned flow model relative to the true Hamiltonian, whileMSE⁡\(H,H^\)\\mathrm\{MSE\}\(H,\\hat\{H\}\)denotes the error of the estimated symbolic approximation\.The results demonstrate that error reduction is not automatic and depends on the ability of the symbolic regressor to construct an accurate approximation within the available hypothesis space\. In theϕ4\\phi^\{4\}experiment, PySR substantially reduces the error relative to the underlying flow model despite the removal of explicit polynomial basis functions, reducing the MSE from0\.220\.22to0\.010\.01\. This indicates that PySR is able to efficiently estimate the missing polynomial structure through compositions of the remaining primitives\. In contrast, EQL introduces several additional spurious terms, resulting in an error substantially larger than that of the underlying flow model\.

A different behavior is observed in the XY experiment, where all trigonometric operators are removed from the primitive library\. In this setting, the target cosine interaction is no longer directly representable and must be approximated indirectly through polynomial compositions\. Both PySR and EQL reduce the error relative to the learned flow model, with EQL achieving a lower symbolic estimation error\. Nevertheless, the reduction is considerably smaller compared to theϕ4\\phi^\{4\}experiment, reflecting the increased difficulty in approximating the functional form\.

In general, these results indicate that the error\-reduction property of symbolic regression is strongly dependent on the degree of misspecification in the primitive set and the expressive power of the underlying SR algorithm\. When the target Hamiltonian remains representable, symbolic regression can substantially reduce modeling errors by estimating compact symbolic representations\. Under stronger misspecification, where the true functional form is absent from the primitive set, symbolic regression can still yield meaningful approximations, although the achievable error reduction is more limited\.