MESSY Estimation: Maximum-Entropy based Stochastic and Symbolic densitY Estimation
Authors: Tony Tohme
Mohsen Sadr
KAMAL YOUCEF-TOUMI
Nicolas Hadjiconstantinou
Venue: Accepted by TMLR ()

Abstract: We introduce MESSY estimation, a Maximum-Entropy based Stochastic and Symbolic densitY estimation method. The proposed approach recovers probability density functions symbolically from samples using moments of a  Gradient flow in which the ansatz serves as the driving force.  In particular, we construct a gradient-based drift-diffusion process that connects samples of the unknown distribution function to a guess symbolic expression. We then show that when the guess distribution has the maximum entropy form, the parameters of this distribution can be found efficiently by solving a linear system of equations constructed using the moments of the provided samples. Furthermore, we use Symbolic regression to explore the space of smooth functions and find optimal basis functions for the exponent of the maximum entropy functional leading to good conditioning. The cost of the proposed method for each set of selected basis functions is linear with the number of samples and quadratic with the number of basis functions. However, the underlying acceptance/rejection procedure for finding optimal and well-conditioned bases adds to the computational cost. We validate the proposed MESSY estimation method against other benchmark methods for the case of a bi-modal and a discontinuous density, as well as a density at the limit of physical realizability. We find that the addition of a symbolic search for basis functions improves the accuracy of the estimation at a reasonable additional computational cost. Our results suggest that the proposed method outperforms existing density recovery methods in the limit of a small to moderate number of samples by providing a low-bias and tractable symbolic description of the unknown density at a reasonable computational cost.

Published in Transactions on Machine Learning Research (01/2024)
MESSY Estimation: Maximum-Entrope based Stochastic and Symbolic densitY Estimation
Tony Tohme
tohme@hmit.cm edu
Massachusetts Institute of Technology, USA.
Mohsen Sadr"
msadr Assid edu
Massachusetts Institute of Technology, USA.
Paul Scherrer Institute, Switzerland.
Kamal Youmef-Toumi
youcef@mitt edu
Massachusetls Institute of Technology, USA.
Nicolas G. Hadjiconstantinou Massuchuselts Institute of Technology, USA.
righ nigh@mit. edu
Reviewed on OpenReview: https: // openreview. net/forunt forum? id- Y2ru0Lu@eS
Abstract
We introduce MESSY estimation, a Maximum-Entropy based Stochastic and Symbolic densit Y estimation method. The proposed approach recovers probability density functions symbolically from samples using moments of a Gradient flow in which the ansatz serves as the driving force. In particular, nects samples of the unknown distribution function to guess symbolic expression. We then show that when the guess distribution has the maximum entropy form, the parameters of this distribution can be found efficiently by solving a linear system of equations constructed using the moments of the provided samples. Furthermore, we use Symbolic regression to explore the space of smooth functions and find optimal basis functions for the exponent of the maximum entropy functional leading good conditioning. The cost of the proposed method for each set of selected basis functions is linear with the number of samples and quadratic with the number of basis functions. However, the underlying procedure for finding optimal and well-conditioned bases adds to the computationa cost. We validate the proposed MESSY estimation method against other benchmark methods for the case of a bi-modal and discontimuous density, as well as a density at the limit of physical realizability. We find that the addition of a symbolic search for basis functions improves the accuracy of the estimation at reasonable additional computational cost. Our results suggest that the proposed method outperforms existing density recovery methods in the limit of small moderate number of samples by providing a low-bias and tractable symbolic description of the unknown density at reasonable computational cost.
construct gradient-based drift-diifusion process that
crreptance/rjjecttin
1 Introduction
Recovering probability density functions from samples is one of the fundamental problems in statistics with many applications. For example, the traditional task of discovering the underlying dynamics governing the corresponding distribution function is strongly dependent on the quality of the densily estima (Rudy el al., 2017). Applications include particle physics (Patrignani et al., 2016), boundary conditions for multi-scale kinetic problems (Frezzotti et al., 2005; Kon al., 2014), and machine learning (Song al., 2020).
'Equal contribution.

---

Published in Transactions on Machine Learning Research (01/2024)
Broadly speaking, two categories of methods have been developed for this task: parametric and non- parametric estimators. While parametric methods assume restrictive ansatz for the underlying distribution function, non paraameeteei methods provide a more flexible density estimate by performing a kernel integra- tion locally using nearby samples. Although non-parametric methods do not need any prior knowledge of the underlying distribution, they suffer from the unclear choice of kernel and its support leading to bias and lack of moment matching. Examples of non-parametric density estimators include histogram and Kernel Density Estimation (KDE) (Rosenblatt, 1956; Jones et al., 1996; Sheather, 2004).
On the other hand, parametric density estimators may allow matching of moments while introducing model- ing error, since a guess for the distribution is required. Parametric distributions include Gaussian, orthogonal expansion with respect to Gaussian using Hermite polynomials (also known as Grad's ansatz in kinetic the- ory) (Hermite, 1864; Grad, 1949; Cai et al., 2015), wavelet density estimation (Donoho et al., 1996), and Maximum Entropy Distribution (MED) (Kapur, 1989; Tagliani, 1999; Khinchin, 2013; Hauck et al., 2008) function among others. Given only the mean and variance, information theory provides us with the Gaussian distribution function as the least biased density, which has been used extensively in the literature. However, including higher order moments in similar way, i.e. moment problem, raises further complications. For example in the context of kinctic theory, Grad proposed a closure that incorporates higher-order moments by considering a deviation from Gaussian using Hermite polynomials, Even though the information from higher moments is incorporated as the parameters of the polynomial expansion in Grad's ansatz, such a formulation suffers from not guaranteeing positivity of the estimated density along with the introduction of bias.
Among parametric density estimators, the Maximum Entropy Distribution (MED) function has been pro- posed in information theory as the least biased density estimate given a number of moments of the unknown distribution (Kapur, 1989). While MED provides the least biased density estimate, it suffers from two lim- itations. First, the distribution parameters (Lagrange multipliers) can only be found by solving a convex optimization problem with ill-conditioned Hcssian (Dreyer, 1987; Levermore, 1996). The condition number increases either by increasing the order the matching moments or approaching the limit of physical re- alizability which motivated the use of adaptive basis functions (Abramov, 2007; 2009). Second, MED only exists and is unique in bounded domains. While existence/ uniqueness is guaranteed for recovering the dis- tribution in the subspace occupied by the samples, the computational complexity associated with the direct computation of Lagrange multipliers has prevented researchers from deploying MED in practice.
Related methods. The problem of recovering a distribution function from samples has been investigated and studied before. We briefly review some of the work most relevant to our paper:
Data-driven marimum entropy distribution function: Several attempts have been made in the literature to speed up the computation of Lagrange multipliers for MED using Neural Networks (Sadr al., 2021; Porteous et al., 2021; Schotthöfer et al., 2022) and Gaussian process regression (Sadr al., 2020). Unfortunately, these approaches are data-dependent with support only on the trained subspace of distributions. Similar to the standard MED and other related closures, the data-driven MED can only handle polynomial moments as input, even though the data may be better represented with moments of other basis functions.
Learning on invertible map: The idea is to train an invertible neural network that maps the samples to a known distribution function. Then the unknown distribution function is found by inverting the trained map with the known distribution as the input. This procedure is called the normalizing flow technique (Rezende & Mohamed, 2015; Dinh et al., 2016; Kingma & Dhariwal, 2018; Durkan et al., 2019; Tzen & Raginsky, 2019; Kobyzev el al., 2020; Wang & Marzouk, 2022). This method has been used for re-sampling unknown distributions, e.g. Boltzmann generators (Noé et al., 2019), as well as density recovery such as Al-Feynmann (Udrescu & Tegmark, 2020; Udrescuu et al., 2020). We note that Al-Feynman does not obtain the density from the samples directly; instead it first fits a density to the samples using the normalizing flow technique, constructs an input/outpur data set, then finds a simpler expression using symbolic regression. While invertible maps can be used accurately predict densities, they can become expensive since for each problem one has to learn the parameters of the considered map via optimization.
2

---

Published in Transactions on Machine Learning Research (01/2024)
Diffusion map: Instead of training for an invertible map, the diffusion map (Coifman et al., 2005; Coifman & Lafon, 2006) constructs coordinates using eigenfunctions of Markov matrices. Using pairwise distances between samples, in this method a kernel matrix is constructed as a generator of the underlying Langevin diffusion process. As shown by Li & Marzouk (2023), one can generate samples of the target distribution using Laplacian-adjusted Wasserstein gradient descent (Chewi et al., 2020). Unfortunately, this approach can become computationally expensive since requires singular value decomposition of matrices of size equal to the number of samples.
The gradient flow method has gained attention in recent years (Villani, 2009; Song et al., 2020; Song & Ermon, 2020). In particular, class of sampling methods has been devised for drawing samples from a given distribution function using Langevin dynamics with the gradient of log-density as the driving force (Liu, 2017; Garbuno-Inigg et al., 2020a;b). Yet, this approach does not provide the density of the samples by itself. In our paper, we benefit. from this formulation to recover the parameters of a density ansatz.
Grodient flow:
Stein Variational Gradient Descent method (SVGD): Given a target density function, the SVGD method as deterministic and aon-parameterized Gradient method generates samples of the target by carrying out a dynamic system on particles where the velocity field includes the grad-log of the target density (similar to the Gradient flow) and a kernel over particles (Liu & Wang, 2016). SVGD has been derived by approximating a kernelized Wasserstein gradient flow of KL divergence (Liu, 2017). While further steps have been taken to improve this method, e.g. SVGD with moment matching (Liu & Wang, 2018), similar to Gradient flow, this class of method cannot be used for estimating the density itself from samples.
Wavelet and Conditional Renormalization Group method: One of the powerful methods in signal processing is the wavelet method (Mallat, 1999), which may be considered as an extension of the Fourier method and domain decomposition. The basic idea is to consider the data on a multiple grids/scales, where the contribution from the smallest frequencies are found on the coarsest mesh and the highest frequencies on the finest mesh. While this method has been extended to finding high dimensional probability density functions (Marchand et al., 2022; Kadkhodaie et al., 2023), is not clear how much bias is introduced by the orthogonal wavelet bases.
KDE via diffusion: In this method, the bandwidth of the kernel density estimation is computed using the minimum of mean integrated squared error and the fact that the KDE is the fundamental solution to a heat (more precisely Fokker-Planck) equation (Botev, 2007; Botev et al., 2010). While improvement has been achicved in this direction, we note that the KDE-diffusion method suffers from smoothing effects which introduce bias. Moreover moments of the unknown distribution are not necessarily matched.
Symbolic regression (SR) is a challenging task in machine learning that aims to identify analytical expressions that best describe the relationship between inputs and outputs of a given dataset. SR does not require any prior knowledge about the model structure. Traditional regression methods such least squares (Wild & Seber, 1989), likelihood-based (Edwards, 1984; Pawitan, 2001), and Bayesian regression techniques (Lee, 1997; Leonard & Hsu, 2001; Tohun et al., 2020) use a fixed parametric model structure and only optimize for model parameters. SR optimizes for model structure and parameters simultaneously and hence is thought to be Na-hard, i.e. Non-deterministic Polynomial-time hard, (Udrescu & Tegmark, 2020; Petersen et al., 2021; Virgolin & Pissis, 2022). The SR problem has gained significant attention over recent years (Orzechowski el. al., 2018; La Cava al., 2021), and several approaches have been suggested in the literature. Most methods adopt genetic algorithms (Koza & Koza, 1992; Schmidt & Lipson, 2009; Tohme et al., 2023). Lately, rescarchers proposed using machine lcarning algorithms (c.g. Bayesian optimization, nonlinear least squares, neural networks, transformers, etc.) to solve the SR problem (Sahoo et al., 2018; Jin et al., 2019; Udrescu et al., 2020; Cranmer et al., 2020; Kommenda et al., 2020; Burlacu et al., 2020; Biggio et al., 2021; Mundhenk et al., 2021; Petersen et al., 2021; Valipour et al., 2021; Zhang et al., 2022; Kamienny et al., 2022). While most SR methods are concerned with finding map from the input to the output, very few have addressed the problem of discovering probability density functions from samples (Udrescu al., 2020).
Symbolic regression:
3

---

Published in Transactions on Machine Learning Research (01/2024)
Our Contributions. Our work improves the efficiency in determining the maximum entropy result for the unknown distribution. We specifically develop new method for determining the unknown parameters (Lagrange multipliers) of this distribution without solving the optimization problem associated with this approach. This is achieved by relating the samples to the MED using Gradient flow, with the grad-log of the MED guess distribution serving as the drift. This results a linear inverse problem for the Lagrange multipliers that is significantly casier to solve compared to the aforementioned optimization problem. We also propose Monte Carlo search in the space of smooth functions for finding an optimal basis function for describing (the exponent of) the maximum entropy ansatz. As selection criterion, rate randomly created basis functions according to the condition number associated with the coefficient. (Hessian) matrix of the inverse problem for the Lagrange multipliers. This helps to maintain good conditioning, which allows us to incorporale more degrees of freedom and recover the unknown density accurately. Discontinuous density functions are treated by considering only the domain supported by data and using a multi-level solution process.
The paper is organized follows. In Section 2 we review the concept of Gradient flow with grad-log of known density as the drift. In Section 3, we show how parameters of guess MED may be found by computing the relaxation rates of the corresponding Gradient flow. Using the maximum entropy ansatz, in Section 4 we derive a linear inverse problem for finding the Lagrange multipliers without the need for solving an optimization problem. In Section 5, we propose a symbolic regression method for finding basis functions that can be used to increase degrees of freedom while maintaining good conditioning of the problem by construction. In Section 6, we propose a generalization of the maximum entropy ansatz that allows including further degrees of freedom in a multi-level fashion. Section 7 presents the complete MESSY algorithm. Tn Scction 8, we validate MESSY by comparing its predictions to those of benchmark density estimators in recovering distributions featuring discontinuities and bi-modality, as well distributions close to the limit of realizability. Finally, in Section 9, we offer our conclusions and outlook.
2 Gradient flow and theoretical motivation
Consider a set of samples of a random variable X from an unknown density distribution function f(x). Let our guess for this distribution function, the "ansatz", be denoted by f(x)..
Instead of constructing non-parametric approximation of the target density numerically from sam- ples of X (like histogram or KDD) and then calculating its difference from the guess density f, in this work we suggest measuring the distance using transport. In particular, we use the fact that the steady-state distribution of X(t which follows the stochastic differential equation (SDE)
dx = V [log (f) + VZaW.
(1)
is the distribution f. Here, W, is the standard Wiener process of dimension dim(x). We note that Eq. (1) is known as the gradient flow (or Langevin dynamics) with grad-log of density as the force. We note that this drift differs from the score-based generative model in Song et al. (2020) where the drift is a function of time, i.e. Va [[)](f(t)]]
The distance of from f may be measured by the time required for the SDE with X(I 0) to reach steady state. Alternatively, may compare the moments computed from the solution to Xt) against the input samples to measure this distance. Both these approaches are subject lo numerical and statistical noise associated with the numerical scheme deployed in integrating Eq. (1). In the next section, we derive an efficient way of computing the parameters of our approximation f based on these ideas. We also show that the transition from f to ff is monotonic.
4

---

Published in Transactions on Machine Learning Research (01/2024)
3 Ansatz as the target density of Gradient flow
According to Ito's lemma (Platen & Bruti-Liberati, 2010) the transition of f to f is governed by the Fokker- Planck equation
0f V f V [s/[]]
(2)
01
~- V log log (f)] V [r[
Proposition 3.1. The distribution function f(t) governed by the Pokker-Plance Eq. (2) converyes to f
as f + co. Purthermore, the cross entropy distance between f and f monotonically decreases during this transition.
Proof. Let us multiply both sides of Eq. (2) by log(f/f) and take the integral with respect to æ in order to obtain the evolution of the cross-entropy S log(f// fod(ff))dz. It follows that
ds
log(s/f) v [ [f dr
dt
V log / An dx - vallog(f/f) - fl
Ve f log(s/f /f) V2 log(f (ff do - [ollog(f//) & V
=0
dim((r)
2 / [VValloo((/)) dx 00.
= Σ
(4)
Here, use the regularity condition that flog(S/f) log(S/f) - 0 as |æ]| + Therefore, given any initial condition for at t= O, the cross-entropy distance between and f following the Fokker-Planck in Eq. (2) monotonically decreases until it reaches the steady-state with the trivial fixed point f - f as t co. For details, see (Liu, 2017).
Instead of soeking solutions of Eq. (2), our approach focuses on working with appropriate empirical moments of this equation, which can be evaluated from the available samples. As will be demonstrated below, this approach lends itself to a very effective method for determining
Let us denote vector of basis functions in ddim(a) by Hxx). By multiplying both sides of Eq. (3) by (x)) and integrating with respect to æ, obtain the evolution equation for the moments, also known the relaxation rates,
d [[HTda) - HV= [Valog(())]] HHVgoVVallghhhhf dx + HV3 [
Assuming that the underlying density is integrable in Rdim(t)) and H as + dc, which implied by the existence of moments, use integration by parts to obtain
dt [/ Hfdo
d
/V vallo((fffs
(6)
Given samples of f, one can compute the relaxation rates of moments represented by Eq. (6) as a measure of the difference between of and f. These relaxation rates can be used as the gradient in the search for parameters of a given ansatz, i.e.
d 4 (H(X(1)))= = d1
(V (V((H(((()) -,,00)/()= + +(V2
g(l)
(7)
5

---

Published in Transactions on Machine Learning Research (01/2024)
In the above, (x(X)) denotes the unbiased empirical mcasure for the expectation of ((X) which is computed using samples of X, for i L..., .N via (o(X)) Φे Em= ((x)).
In what follows develop approach that uses g(t) the gradient of optimization problem to bring computational benefits to the solution of the maximum entropy problem.
4 Maximum Entropy Distribution an ansatz for the gradient flow
In this work, use the maximum entropy distribution function as our parameterized ansatz for f, i.e.
cxp f(x) = Z exp (A H())
(8)
where Z exp(A (x))d is the normalization constant. The motivation for choosing this family of distributions is the fact that this is the least-biased distribution for the moment problem, provided the given moments are matched.
Definition 4.1. Moment problem
The problem of finding a distribution function f(x) given its moments J (())f(z)d = μ for the vector of basis functions H (a:) will be referred to as the moment problem.
In particular, the density in Eq. (8) the extremum of the loss functional that minimizes the Shannon entropy with constraints on moments z using the method of Lagrange multipliers, i.c.
f(x) = arg min C(F(x)) FeK
(9)
No
where C[F(x)) = / F(xx) log log(F(x))
E Ai (ATla))oo)daa-iit -i(())
(10)
Here K denotes the space of probability density functions with measurable moments; see (Kapur, 1989) and Appendix A for more details. In this paper, we denote the number of considered basis functions by No, while Nm denotes the highest order of these basis functions. For instance, in the case of traditional one-dimensional random variable where polynomial basis functions are deployed, ie. H = [x,z?,.... x Nm) we have Nm No.
Here, we the following definition for the growth rate of basis function. Definition 4.2. Growth rate of n-th order
A function (x)) has the growth-rate of n-th order [f((x)) < Cw" for all x 2 where E RR and xo E R. This often denoted by (xx)) = O(=").
We note that the moment problem for the MED is reduced to the following optimization problem by substi- tuting the extremum 8 back in the objective functional equation 10, see e.g. Abramov (2006).
Definition 4.3. Standard dual optimization problem of Mazimum entropy distribution function
Given moments μt the Lagrange multipliers À of the maximum entropy distribution are the solution to the following unconstrained optimization problem
X arg min log /eeppl Har))da - x xx. AERN
(11)
Clearly, the standard optimization problem for finding Lagrange multipliers is nonlinear and requires deploying iterative methods, such as the Rephon-Raphson method detailed in A.
Instead, using the Gradient flow, we find an alternative optimization problem for finding Lagrange multipliers which is linear and simple to compute from given samples.

---

Published in Transactions on Machine Learning Research (01/2024)
Substituting Eq. (8) for f in Eq. (7) results in the relaxation rate
dim( (+)
(im
g(t) =
V. (((()))
(x)))
(12)
where @ indicates the outer product. Let us define the matrix L ME as
dim((z)
L ME (1)
Σ (H((H(((())
V. V (((())))
(13)
Definition 4.4. Optimization problem of Marimum entropy distribution function via Gradient flow
Given samples X of the target distribution f, Lagrange multipliers of maximum entropy distribution estimate f is the solution to the following unconstrained optimization problem
λ A
dim(w)
(((()))) 2).
À arg min AERN
MP MB(())
Σ
(14)
where (:) denotes the Frobenius inner (or double dot) product, i.e. A B = Eis Aij for given two- dimensional tensors A and B.
We note that the matrix ME is the Hessian of the optimization problem with gradient given by Eq. (12) which is positive definite, making the underlying optimization problem convex.
Proposition 4.5. The Hessian matris ME is symmmelrie posilive definile. As a result, the oplimization problem with gradient given by Eq. (12) and Hession matria given by Eq. (13) is strictly conver.
Proof. Clearly, the Hessian matrix defined by Eq. (13) is symmetric, ie. LNJ ME Vi,j = I,.., We further note that this matrix is positive definite, i.c. for any non-zero vector w E N, we can write
dim((w)
w L (I)
V (()())) V., [H((())))
(15)
dim(w)
V., ([ ?) >0.
(16)
Given the Hessian is symmetric positive definite, we conclude that the underlying optimization problem is convex (Chong & Zak, 2013).
When the matrix ME is well-conditioned, we can directly compute the Lagrange multipliers using samples, i.e. a linear solution to the optimization problem in Def. 4.4. This can be achieved by solving Eq. (12) for the Lagrange multipliers
dim(ar)
L ME (E g(t) - Σ V2 ((((())))
(17)
for given relaxation rate g.
We proceed by noting that a convenient way for determining the parameters of f is set f f(t = 0) in the above formulation, or other words, require that the given samples are also samples of as given. This corresponds to the steady solution of Eq. (12), namely g - 0, which implics the remarkably simple result
dimn(w))
V2 H H H[{x(=0)))) 0)1))
A = -CZMEE-- -1
Σ
(18)

---

Published in Transactions on Machine Learning Research (01/2024)
which implies closed-form solution for the Lagrange multipliers through the above lincar problem.
While Eq. (18) analytically recovers the Lagrange multipliers À directly from samples of X, it still requires inverting the matrix L ME which may be ill-conditioned (Abramov, 2010; Alldredge et al., 2014). This means that the resulting Lagrange multipliers may become sensitive to noise in the samples and the choice of the basis functions. In order to cope with this issue, we propose computing as outlined below.
Orthonormalizing the basis functions. We construct an orthonormal basis function with respect X ~ using the modified Gram-Schmidt algorithm as described in Algorithm 1. We deploy the orthonor- mal basis functions from the Gram-Schmidt procedure to construct V [=] i.e. V=[H]] is the input to Algorithm 1, and by integration we obtain H This leads to well-conditioned matrix L LM, since the resulting matrix should be close to identity ME ~ I with condition number cond ME) x 1 subject to round-off crror. We note that the cost of this algorithm is quadratic with the number of basis functions and linear with the number of samples.
(E)
Algorithm 1: Modified Trran-Schmidd: Given a vector of basis functions &, this algorithm constructs an orthonormal basis functions o- with respect f such that (x*(X) x (X)) I using the modified Gram-Schmidt procedure (Giraud et al., 2002; Abramov, 2010).
Input: Ø
Initialize &² + Φ; for i= 1..,.,.imm(
do of = /
(x))²);
for j = i 1,.... dim(mimmm) (6) do
0 + 0j - (0) (x) (x)) (X) //
end
end
Return o-
4.1 Comparing the proposed formulation to standard Maximum Entropy Distribution
Here we point out several advantages of using the proposed loss function compared to the standard maximum entr closure.
A closed-form solution. By setting the relaxation rate of the moments zero, the Lagrange multipliers can be computed directly from samples X ~ f. i.e. by solving the system 18, without the need for the lino-scarch associated with the Newton method of solving the optimization problem of standard MED, i.e. Def. 4.3. This is significant improvement compared to standard MED, where Lagrange multipliers are estimated iteratively see Eq. (A.5) and Algorithm 4 in Appendix A.
Avoiding the curse of dimensionality in integration. In the proposed method, the com- putational complexity associated with the integration only depends on the number of samples, and not on the dimension of the probability space. The proposed method takes full advantage of having access to the samples the unknown distribution function. This is in contrast to the standard Newton- Raptson method of solving optimization problem 4.3 where samples of the initial guessed MED corresponding the initial guessed À is not available and one runs the curse of dimensionality in compulation of gradient and Hessian.
In particular, we compute the orthonormal basis function, gradient, and Hcssian using the samples of X. This use of the Monte Carlo integration method avoids the curse of high dimension- ality associated with the conventional method for computing Lagrange multipliers. By deploying the Law Large Numbers (LLN) in computing integrals, the proposed method benefits from the well-known result that the error of Monte Carlo integration is independent of dimension. In particular, given N independent, identically distributed d-dimensional random variables XI...,,X, with probability density f(x), where X f, the variance of the empirical estimator of moment
8

---

Published in Transactions on Machine Learning Research (01/2024)
((x) of f, i.. [+((( E/= (X,)/ is
(ar(E-(o(X))) (x))) Var()()))
(19)
N
which is independent of the dimension d= dim(z). Therefore, for a required ratio of variance in prediction and variance of the underlying random variable, i.c. Var(EA( ())))///((dbbb the cost of integration is (ssN). The factor denotes the cost of computing o( (X) for one sample.
This is a considerable advantage compared to the standard approach of finding the Lagrange multipliers of MED where the cost associated with integration of order O(N") where N is the number of discretization points in each dimension. We remind the reader that in the standard Newton- Raphssn method of finding the Lagrange multipliers, one updates the guessed by
(20)
where the gradient g and Hessian need to be computed during each iteration, see section and reference therein for details. Since samples of the guessed distribution, i.e. guessed X, are not available, computation of the gradient and Hessian can become expensive. Specifically, one has to either generate samples of the guess distribution in each intermediate step of the Newton-Raphan Raphson method which adds complexity, deploy a deterministic integration method with cost O(N4) where N is the number of discretization points in each dimension and d = dim(x).
We further point out that since in both the proposed solution and the standard iterative method a matrix needs to be inverted, the cost in both algorithms scales ((N3) with number of basis functions (moments) regardless of the cost associated with integration. For example, in case of using monomials up to 2nd order in all dimensions as basis functions i.c. H [I1, ... Fdd 27 2222, .... x2). there are H = d d(d + 1)/2 unique bases, leading to complexity O ((d + d(d + 1)/2)9) for solving the lincar system 18.
Relaxed existence requirements. Since the proposed approach directly finds the Lagrange mul- tiplicrs for the realizable moment problem lincarly using Eq. (18), it avoids the problem of possible non-realizable distribution estimate with intermediate Lagrange multipliers that is present in the iterative methods. This is another advantage compared to the standard MED optimization problem, Eq. (4.3), where the line search Eq. (20) may fail as the distribution associated with the intermediate may not exist (not integrable). This is a common problem when finding Lagrange multipliers for the moment problem close to the limit of realizability when the condition number of the Hessian bocomes large and the itcrative Newton-Raphsee method fails, c.g. see Alldrodge al. (2014).
5 Symbolic-Basedd Maximum Entropy Distribution
In the standard moment problem it is common to consider polynomials for the moment functions in H, i.e. H = x2..?..... even though other basis functions may better represent the unknown distribution. Additionally, such polynomial basis functions are notorious for resulting in ill-conditioned solution processes, For these reasons, we introduce symbolic regression approach to introduce some diversity and ultimately optimize over our use of basis functions. As we will see in the next section, adding the symbolic search to our MED description improves the accuracy, convergence, and robustness of the density recovery problem.
Before diving into the proposed method, we first briefly review the general task of symbolic regression. Definition 5.1. Symbolic Regression (SR) problem
Given a metric C and a dataset = (z,, ys y consisting of N independent identically distributed (t.i.d.) paired samples, where x; ddi((=) and ys E R, the SR problem searches in the space of functions S defined by a set of given mathemalical fanctions (e.g., cos, sin, exp, In) and arilhmetic operations (e.g. +, -, x, + for a function 60 (x) which minimizes EN C(y,, (2:)) where so E S.
9

---

Published in Transactions on Machine Learning Research (01/2024)
In order to deploy the SR method for the density recovery, need to restrict the space of functions S to those which satisfy non-negativity, normalization and existence of moments with respect to the vector of linearly independent (polynomial) basis functions R. The space of such distributions can be defined as
S/RR = { ((x)) ((()) ES f(x) 2 0 Vo E (ddm(a)) / Rdim(a)
[ax) dæ < to0].
(xx) dæ - 1,
(21)
Rddiu(te)
In order ensure non-negativity, motivated by the MED formulation, we consider to be exponential, i.c.
((x) x exp (9(x))
log (f(z)) x G(x))
(22)
where ((x) is an analytical (or symbolic) function of æ = [x1, T2...., Tdim(#) While the non-negativity is guaranteed, existence of moments needs to be verified when a test function for G(x)) is considered. As our focus in this paper on the maximum entropy distribution function given by Eq. (8), we consider (xx)
have the form
No
((x) = X H()) = Σ (ar)
(23)
Now we proceed to provide a modified formulation for SR tailored to our MED problem.
Definition 5.2. Symbolic Regression for the Marimum Entropy Distribution (SR-MED) problem
Given measure of difference between distributions C (e-g- KL. Divergence) and dataset DD (X.) i=l consisting of N i.i.d. sumples, where X; E (ddi(a), the SR-MED problem searches in the space SNo for No. basis functions subject to f SHR which minimizes L.
Here, we deploy continuous functions consisting of binary operators (e.g. +, -, x, +) or unary functions (c.g. cos, sin, exp, log) to fill the space 5 As in most of the SR methods, we encode mathematical expressions using symbolic expression trees, a type of binary tree, where internal nodes contain operators or functions and terminal nodes (or leaves) contain input variables of constants. For instance, the expression tree in Figure 1 represents x² cos(x)) In this paper, we perform Monte Carlo symbolic search in the space of smooth functions (by generating random expression trees) to find vector of basis functions H that guarantees acceptable ( ME), by rejecting candidates that do not satisfy this condition. In our search, we do not consider test. basis functions with odd growth rates which lead to non-realizable distributions.
cos
(-)²
Figure 1: Expression tree for x² x cos(z).
6 Multi-level density recovery
We further improve our proposed method by introducing a multi-level process that improves our prediction as the distribution becomes more detailed. The goal is to obtain a more generalized MED estimate with the form
Nr.
f(xx = Σ m [] (()(())))
(24)
1=1
where
f()(x)) = Z exp A..H!((.))
(25)
())l denotes the level index, d is the normalization factor of density at level 1, N,, the number of levels
considered and m!!) indicates the portion of total mass that covered by fil. We note that this multi-level approach is recursive and can be described follows:
10

---

Published in Transactions on Machine Learning Research (01/2024)
Step 1: Find MED estimate f(i at level 1. At level 1, first pick basis function H by solving the SR-MED problem detailed in Def. 5.2. Then, orthonormalize the basis function with respect to the distribution of the samples using Gram Schmidt's procedure as outlined in Algorithm 1.
Step 2: Removing subset of samples covered by f Hcre, attempt to find and remove a subset samples D mask representing a fraction of the mass, m = |D muskkl l/ that can be estimated by our estimated f at this level. To this end, we deploy mask from the remaining samples D
cepplaanncerrjecttinn with
probability flll f to find and remove D
Step 3: Repeat steps 1-2 for the next level I++ until almost no samples are left. Repeat steps 1-2 with the remaining uncovered samples (which constitutes the next level) until there are (almost) no uncovered samples. The resulting total distribution weighted sum of the estimates from each level.
In Algorithm 2, we detail a pseudocode for our devised multi-level process. As we will see in the next section, our proposed multi-level recursive mechanism improves overall performance, and elegantly describe details of multi-mode distributions.
Multi-level, symbolic and recursive algorithm for density recovery. Here, D denotes the set of samples at level I and u is random variable that is uniformly distributed in (0, 1), i.c. ~ (([0 1)).
Algorithm 2:
Input: D
(X,) No NL
for / 1, ,..g Nr. do
Sample random basis functions HIl that satisfies Def. 5.2 starting from polynomials in level Compute (z) given DIIl using Algorithm 1; Dl1
mask + (D(t) (x) fhist (X) > u) where u u( 1)):
ml - -DD11 1Dmssk// I/|D|:
if E Dli Duakk (D] then
D mask DH!:
// Mask available samples
m!ll_ + |D mask I/D]:
No.t + 1;
break;
// Terninate the process // The uncovered samples left for the next level
else
Di++1 Dill maski
end
end
Return
N No A
(xx
! // D
No.t total number of recursive calls
7 Algorithm for MESSY estimation
The complete MESSY estimation algorithm summarized in Algorithm 3. Within the iteration loop, follow- ing each application of the multi-level, symbolic, and recursive density recovery summarized in Algorithm 2, we introduce maximum-cross entropy distribution (MxED) correction step (see Appendix B for details) reduce any bias in our prediction for from the former.
11

---

Published in Transactions on Machine Learning Research (01/2024)
Finally, after completing the desired number of iterations, the algorithm returns the candidate density with the smallest KT. Divergence given by
KL(SI1 - [(x)) log
f(x)
dæ
(26)
f(x))
/
f(x) log (f(z)) + / | f(x) log (f(z)) dx
(27)
/ [())log (f(z))da
(log (f(x))
(28)
constant with respect to
In other words, we use - (log (f(x))) as our selection criterion.
Algorithm 3: Pseudocode of the proposed MESSY cstimation method. Herc, R the vector of lincarly independent (polynomial) basis functions used in the moment matching procedure of MxED. Here, for MESSY-S the number of basis functions N, is sampled uniformly from the sample space f2,, e.g. here we use S2 = (,.... 8) unless mentioned otherwise.
Input: = [ ON,, Nm, Nuers
Initialize f(() = for = Niters
for i i= to Niters do
if i> then
Sample No ~((N,): end
Find ff using multi-level, symbolic and recursive Algorithm 2 for density recovery;
Generate samples of Y f;
Apply boundary condition (bounded/unbonn /unbouded) to f:
Correct f using MxED (Algorithm 5) given samples Y as prior and B[R(X)) target moments;
f(0 + f:
end
JMESSY (());
JMESSY argmin ff( i21 Return fMESSY-P and JMESSY-S.
[LT(flI))
The MESSY algorithm comes in two flavors: MESSY-P, which considers only polynomial basis functions for H, and MESSY-S which includes optimization over basis functions using the SR algorithms outlined above. In fact, by convention, the SR algorithm in MESSY-S starts its first iteration using polynomial basis functions up to order Nm, as the sample space of smooth functions. In other words, MESSY-P a special case of MESSY-S with Niter = 1. In the remaining iterations of MESSY-S, we perform the symbolic search in the space of smooth functions of order Nm find N₂ bases that provide managcable cond L MB), as discussed in Section 6.
In addition, provide the option to enforce boundedness of f on the support that is specified by the user, i.c. letting f(xx) = 0 for all æ outside the domain of interest. This allows us to recover distributions with insyyn at the boundary which may have application in image processing.
We also provide an option to further reduce the bias by minimizing the cross-entropy given samples of bounded/unboung /nbounddd multi-level estimate as prior and moments of input samples as the target moments (sce Appendix B for more details on the cross-entropy caleulation). For this optional step, we generate samples of and match the moments of polynomial basis functions up to order N,, Since the solution at cach level of f is close the exact MED solution, the optimization problem associated with the moment matching procedure of the cross-entropy algorithm converges very quickly, i.e. in few iterations, providing us with
12

---

Published in Transactions on Machine Learning Research (01/2024)
correction that minimizes bias along with the weighted samples of our estimate the by-product. We note that in general the order of the randomly created basis function during the MESSY-S procedure may be different from the one used in the cross-entropy moment matching procedure.
8 Results
In this section demonstrate the effectiveness of the proposed MESSY cstimation method in recovering distributions given samples, using number of numerical experiments, involving a range of distributions ranging from multi-mode to discontinuous. For validation, we provide comparisons with the standard KDE with an optimal bandwidth obtained using K-fold cross-validation with K = 5. We consider a uniform grid with 20 elements for the cross-validation in a range that spans from 0.01 to 10 on logarithmic scale, We also compare with cross-entropy closure with Gaussian as the prior (MxED) using Newton's method. We note that while the standard maximum entropy distribution function differs from MxED as the latter incorporates a prior, we intentionally use MxED benchmark instead because the standard approach can be extremely expensive.
Unless mentioned otherwise, we report error, time, and KL Divergence by ensemble averaging over 25 for different sets of samples. Furthermore, in the case of MESSY-S we perform Nitors = 10 iterations, and we consider (+,-,X) x) operators and (cos, sin) functions. Here we report the execution time using a single-core CPU for each method. Typical symbolic expressions of density functions recovered by MESSY for the test cases considered here can be found in Appendix E. Furthermore, in Appendix C, we perform an ablation study that shows the importance of each component of the MESSY algorithm.
8.1 Bi-modal distribution function
For our first test case, we consider a one-dimensional bi-modal distribution function constructed by mixing two Normal distribution functions NNx o), i.c.
f(x) = (x (1) (1-aa)
(29)
with a = 0.5, means μ1 = -0.6 and #2 = 0.7, and standard deviations 0₁ = 0.3 and 02 = 0.5.
Figure compares results from MESSY, KDE and MxED for three different sample sizes, namely 100, 1000, and 10, 000 samples of f. For MxED and MESSY-P, we N₂ Nm 4. In the case of MESSY-S, randomly create No basis functions which are O(24) (where No is sampled uniformly within 8)). Both MESSY results are subject to a cross-entropy correction step with N, = 4 polynomial moments. Clearly, the MxED and MESSY methods provide a reasonable estimate when a small number of samples is available; where KDE may suffer from bias introduced by the smoothing kernel. As a reference for the reader, we also compared the oulcome density with a histogram in Appendix D.
In order to analyze the error further, Fig presents the relative error in low and high order mo- ments, KL Divergence, and single-core CPU time the measure of computational cost for considered methods. Given optimal bandwidth is found, the KDE error can only be reduced by increasing the number of samples. However, maximum entropy based estimators such as MxED and MESSY provide more robust estimate when less samples are available. We point out that the convergence of the cases where only moments of polynomial basis functions are considered, i.e. MxED and MESSY-P, relies on the degree of the polynomials and not the number of samples. On the other hand, the additional scarch associated with MESSY-S returns more appropriate basis functions for given upper bound on the order of the basis functions.
Next, we perform convergence study on 10,000 samples and show that the parametric description converges to the solution when its degrees of freedom are increased. In Fig. 4, we show that both MESSY-P and MESSY-S converge to the true solution by increasing either the order of polynomial basis function, or the number of basis functions, respectively. In the case of MESSY-S, we generated symbolic expressions that are O22)). The improved agreement compared to the MESSY-P case highlights the benefit derived from non-traditional basis functions that may better represent. the data.
13

---

Published in Transactions on Machine Learning Research (01/2024)
The PD7
0.7
Thue PDF
Thue PDF
60E
KDE
KDE
0.8
MitE
0.6
Metil
HoED HESSY0
HESSYP
MESSVP
0.4
0.4
0.4
0.1
0.1
9.1
(a) 100 samples
(b) 1,000 samples
(c) 10,000 samples
Figure 2: Density estimation using KDE, MxED, MESSY-P, and MESSY-S given (a) 100, (b) 1,000, and (c)
10,000 samples.
Relative Error in Two Higher Order Moments
10
KDE
Relalive Error first Four Moments
MESSYP
10-4
10
10:4
10-13
10
1
I
100
1000
10000
1.00
1000
10000
(a)
(b)
102
10a
KL Divergence
Time (secands)
201
10a
100
1000
10000
100
1000
10000
(c)
(d)
Figure 3: Comparing the relative error in (a) the first four moments, (b) two higher order moments (i.e. fifth and sixth moments), (c) KL Divergence, and (d) the execution time for KDE, MxED, MESSY-P, and MESSY-S recovering distribution function for different sample sizes. Here, the error bar (in black) corresponds to the standard error of the empirical measurements.
As shown in Fig. 5, the MESSY-S procedure results in belter-conditioned ME matrices than the MESSY-P for the same degrees of freedom. However, the search for good basis function increases the computational cost. In each iteration of the search for basis functions, the MESSY-S algorithm may reject symbolic basis candidates based on the condition number of the matrix MF, other words, the improved performance associated with MESSY-S comes at some increased computational cost.
14

---

Published in Transactions on Machine Learning Research (01/2024)
0.7
M
Ng 2
0.6
M
Ng = 4
Mi
N,=
0.5
M
M; 10
N, =
0.4
0.4
[(x)
0.3
0.3
0.2
0.1
0.1
0.0
-2
3
-3
2
2
37
(a) MESSY-P
(b) MESSY-S
Figure 4: Convergence of MESSY estimation to larget distribution function by (a) increasing the order of polynomial basis functions for MESSY-P or (b) increasing the number of randomly selected symbolic basis functions with Nm 2 for MESSY-S.
MESSY-P
10a
MESSY-S
Condition Number
10a
KL Divergence
102
Time (seconds)
10a
10a
10a
10"
10a
2
Na
N₂
N₂.
Figure KL Divergence, execution time, and condition number against the degrees of freedom, i.e. the order of polynomial basis functions for MESSY-P or the number of symbolic basis functions with Nm 2 for the MESSY-S estimate.
8.2 Limit of realizability
One of the challenging moment problems for maximum entropy methods the one involving distributions near the border of physical realizability. In the one-dimensional case with moments of the first four monomials [x, 2 23 21 as the input, the moment problem is physically realizable when
f (aA(())des (f²((o)an) +
(30)
The moment problem with moments approaching the equality in Eq. (30) is called limit of realizablity (McDonald & Torrilhon, 2013; Akhiczer & Kemmer, 1965). We consider samples from distribution in this limit as our Lest case here, since the standard MED cannot be solved due to an ill-oonditioned Hessian (see Abramov (2007); Alldredge et al. (2014)).
In Fig. 6, we depict the estimated density of bi-modal distribution in this limit given its samples with moments (X) 0, (x2) 1, (X3) -2100 and (X1) 5.42. Here, compare the density obtained
15

---

Published in Transactions on Machine Learning Research (01/2024)
using KDE, MxED, MESSY-P, and MESSY-S to the histogram of samples. In this example, obtained the MESSY-S estimate by searching in the space of smooth functions with N₂ E (2,..., 8) basis functions and compare polynomial and symbolic basis functions of order 2 and 4.
In Fig. 7, we compare the KL Divergence, the execution time, and the condition number for each method. While KDE suffers from over-smoothing and MxED MxED/MESSY /M/SS--P require least N₂, (and consequently Nm 4, resulting in a stiff problem with large condition number), MESSY-S can obtain accurate den- sity estimates by using unconventional basis functions with Nm 2, thus maintaining manageable condition number.
Samples
Samples
5
KDE
KDE
MxED
MxED
MESSY-P
MESSY-P
MESSY-S
MESSY-S
2
-3
2
-1
0
-3
-2
-1
1
2
Figure Estimating density for a case of distribution near the limit of realizability using KDE, MxED, MESSY-P, and MESSY-S. The solutions of MxED, MESSY-P, and MESSY-S are obtained using basis functions of second (left) and fourth (right) order.
10a
KL Divergence
Time (seconds)
Condition Number
107
10a
103.
10a
10a
10a
10a
10-1
KDE
MXED (2)
MXED (4)
MESSY-P (2)
MESSY-P (4)
MESSY-S
KDE
MxED (2)
MxED (4)
MESSY-P (2)
MESSY-P (4)
MESSY.S (2)
MxED (2)
MxED (4)
MESSY-P (2)
MESSY-P (4)
MESSY-S
Figure 7: Comparing KL Divergence, execution time, and condition number of KDE, MxED, MESSY-P, and MESSY-S for an unknown distribution near the limit of the realizability. Here, we consider polynomial basis functions of second and fourth order for MxED and MESSY-P denoted by MxED (2), MxED (4), MESSY-P (2) and MESSY-P (4), respectively. In MESSY-S, we consider symbolic basis functions of sccond order only which we denote by MESSY-S (2).
16

---

Published in Transactions on Machine Learning Research (01/2024)
8.3 Discontinuous distributions
We now highlight the benefits of using MESSY estimation with piecewise conlinuous capability for recov- ering distributions with discontinuity at the boundary. As an example, let us consider the exponential distribution with probability density function given by
ae
i otherwise
f(x)
(31)
0
with a 1.
Given 10, 000 samples of this distribution, in Fig. we compare KDE, MxED, and the proposed MESSY-P and MESSY-S methodologies. In the case of MxED and MESSY-P we consider second-order polynomial basis functions, and for MESSY-S we search the space of smooth functions for N₂ E (2.... 8) symbolic basis functions of order O(x²). For MESSY-P and MESSY-S, apply the boundary condition
f(z)= =0 Vx < min
(32)
By providing information about the boundedness of the expected distribution, we enable MESSY to ac- curately predict densities with discontinuity near the boundary. As it can be seen clearly from Fig. 8, in contrast MxED, both MESSY-P and MESSY-S provide accurate predictions by taking advantage of the information about the boundedness of the target density.
1.0
Thue PDF
KDE
MxED
MESSY-P
0.6
0.4
0.0
2
Figure 8: Estimating density of exponential distribution function from its samples using KDE, MxED, MESSY-P, and MESSY-S. For MxED, MESSY-P and MESSY-S, with Nm 2.
The KL Divergence score and execution time for each method is shown in Fig. 9. These figures show that MESSY-P and MESSY-S provide a more accurate description compared to the KDE estimate, albeit at a higher computational cost.
17

---

Published in Transactions on Machine Learning Research (01/2024)
376,26
10a
95,17
36,57
KL Divergence
Time (seconds)
10a
10a
0.98
10a
46.11
38.62
27.93
10-1
0.10
MxED
MESSY-P
MESSY-S
MxED
MESSYPP
MESSY-S
Figure 9: KL Divergence and execution time for KDE, MxED, MESSY-P, and MESSY-S estimation of exponential distribution function given 10,000 samples.
8.4 Two-dimensional distributions
In this section the proposed MESSY methodology is applied to the estimation of 2-dimensional distributions. Let us consider a multivariate Gaussian distribution in two dimensions, i.e. X = (X1, ~ N(u, 2) with probability density
1 ( 11) 10)
N(x; E) (2m)d/2 /2e)d/2/det(1)) (E) exp
(33)
where d 2, x = (x1,22) IR? is the mean vector, and Σ is the covariance matrix. In this example,
1 0.5
E =
(34)
0.5
Let us also consider a more complex 2D distribution by multiplying two independent 1D distributions, namely the exponential and Gamma distributions, i.c. X1 Exp(A), X2 r(k, 0), and X (X1, Exp(A) 0). In particular, we consider samples of the joint pdf given by
de 1e-
Am
6-200 T(kjo*
3 22,0222
f(11, X, 0) f(1;; ()5 (, k,
10
(35)
otherwise
where À > the rate parameter, the gammna function. In our example, we use À 2, = 3, and 0.5.
is the shape parameter, 0>> 0 is the scale parameter, and r(.) is
Fig. 10 compares the true distribution against the KDE, MESSY-P and MESSY-S estimate given 10 samples of the distribution. In cases of MESSY-P and MESSY-S, consider basis functions up to 2nd order in each dimension. The figure shows that KDE, MESSY-P and MESSY-S provide reasonable estimate for the multivariate normal distribution. However, among all considered estimators for the case of the Gamma-cxponentia distribution, MESSY-S scoms to provide a slightly better estimate, i.c. it captures the peak in the most probable region the distribution.
8.5 Cost of MESSY-P with respect to dimension
As discussed in section 4.1, for a given moment problem (fxxd vector of basis functions), the cost of evaluating Lagrange multipliers in the proposed method is linear with respect to number of samples and independent of dimension benefit of Monte Carlo intogration. However, similar to the standard MED approach, the proposed MESSY estimate requires the inversion of Hessian matrix addition to integration, MESSY cost estimate also includes solution of linear system of equations with
N,x,, Therefore,
18

---

Published in Transactions on Machine Learning Research (01/2024)
True Distribution
KDE
MESSY-P
MESSY-S
Gaussian dist.
0.168
2-
2
0.144
0.120
0
0.096
-2
-2
-2
0.024
0.000
-2
-4
-2
-2
True Distribution
KDE
MESSY-P
MESSY-S
Gamma-ena exmma dist.
Samples
4-
4-
0.90
0.75
22
3-
0.60
2
2
0.45
1
1
1
0.00
301
21
Figure 10: Estimating the density of samples in two dimensions using KDE, MESSY-P, and MESSY-S. We use 10,000 samples from a Gaussian (top) and the Gammaaexponential density defined in Eq. (35) (bottom). For better visualization, we only show 100 random samples. We also report that the KL-divergence for MESSY-S, MESSY-P and KDE are on the same order.
cost of order O(Ng)- In this example, we consider the MESSY-P estimator with basis functions H [I1,,, Td, 122, x for estimating target d-dimensional multivariate normal distribution Eq. (33) with mean μ
1 ~~N(0 1 n
(36)
and covariance matrix
Σ I + U
(37)
0
<<j otherwise
and S~N(0.
where
Un
(38)
E
Using 2nd order polynomial basis functions, the number of unique basis functions N, d(d + 1)/2, leading lo cost. of order O(d + d(d++()))))))))) illustrate this by showing the normalized execution time (normalized by the time for d = and given sample size) against the dimension of the target density given N E {100, 200, 400, 800, 1600, 3200} samples. Furthermore, in order to highlight the linear dependence of cost with number of samples N for a fixed moment problem, in Fig. 11-(b) we also report the normalized execution time of computing Lagrange multipliers (normalized by the time for N 100 and dimension of target density) as a function of number of samples.
1(d+))/2){) 1) for d-dimensional probability space. In Fig. 11-(a), we
Since the most expensive component of MESSY, i.c. the computation of Lagrange multipliers, has been carried out single computer core with a direct solver, we confirm that MESSY algorithm can be used for density recovery in large number of dimensions at feasible computational cost, regardless of the scaling with dimension. The limiting factor appears to be the storage associated with solving the linear system. In case a larger system than the one studied here are of interest, one may use itcrative solvers, see Saad (2003).
19

---

Published in Transactions on Machine Learning Research (01/2024)
d=1
= 200
10a
N = 400
d
N = 800
d ==6
N 1600
d=8
10a
3200
104
10
r/r(d-1)
[d'+ d(d 1)/21
100)
x
10a
10
10c
3
5
910
100
200
800
1600
3200
(a)
(b)
Figure 11: Relative execution time of computing Lagrange multipliers given N {100, 200, 400, 800, 1600, 3200} samples of d-dimensional multivariate normal distribution func- tion where d = 1,... 10. Left (a) shows the normalized execution time versus dimension and right (b) the normalized execution time versus number of samples. The execution time is computed on a single core and single thread of 2.3GHz Quad-Core Intel Core i7 processor and averaged over 5 ensembles.
9 Conclusion and Outlook
We present new method for symbolically recovering the underlying probability density function of a given finite number of its samples. The proposed method the maximum entropy distribution an ansatz thus guaranteeing positivity and least bias. We devise a method for finding parameters of this ansatz by considering Gradient flow in which the ansatz serves as the driving force. One main takcaway from this work is that the parameters of the MED ansatz can be computed efficiently by solving a linear problem involving moments calculated from the given samples. We also show that the complexity associated with finding Lagrange multipliers, the most expensive part of the algorithm, scales linearly with the number of samples all dimensions. On the other hand, the cost of inverting the Hessian matrix with dimension NoNN is O(N)), where the number of basis functions N₂ grows with dimension. Overall, our numerical experiments show that densities of dimension 10 can be recovered using single core of typical workstation. It is worth noting that since the whole probability space is represented with a few parameters, the MESSY representation can reduce the training cost associated with the PDE-FIND Rudy et al. (2017) and SINDy Brunton et al. (2016) method, where in the latter the regression is done on a data set that discretizes the whole space with a large number of parameters that need to be fit. Further detailed analysis the tradc-off between cost and accuracy inferring densities high-dimensional probability spaces will be addressed future work.
The second main takeaway from this work is that accurate density recovery does not necessarily require the use of high-order moments. In fact, increasing the number of complex but low-order basis functions leads to superior expressiveness and better assimilation of the data. For this reason, the proposed method is equipped with a Monte Carlo search in the space of smooth functions for finding basis functions to describe the exponent of the MED ansatz, using KL Divergence, calculated from the unknown-disttrbbttion samples, as an optimality criterion. Discontinuous densities are treated by considering piece-wise continuous functions with support on the space covered by samples.
We validate and test the proposed MESSY estimation approach against benchmark non- non-parametric (KDE) and parametric (MxED) density cstimation methods. In our experiments, we consider three canoni- cal Lest. cases; a bi-modal distribution, a distribution close to the limit of realizability, and a discontinuous
20

---

Published in Transactions on Machine Learning Research (01/2024)
distribution function. Our results suggest that MESSY estimation exhibits several positive attributes compared to existing methods. Specifically, while retaining some of the most desirable features associated with MED, namely non-negativity, least bias, and matching the moments of the unknown distribution, it outperforms standard maximum-ontropp-baase
approaches for two reasons. First, uses samples of the target distribution the evaluation of the Hessian, which has a linear cost with respect to the dimension of the random variable. Second, the resulting linear problem for finding the Lagrange multipliers from moments is significantly more efficient than the Newton line search used by the classical MED approach. Moreover, our multi-level algorithm allows for recovery of more complex distributions compared to the standard MED approach. Combining the efficient approach of finding maximum entropy density via linear system with the symbolic exploration for the optimal basis functions paves the way for achieving low bias, consistent, and expressive density recovery from samples.
Although the proposed MESSY estimate alleviates a number of numerical challenges associated with finding the MED given samples, it cannot expand the space of the MED solution. In other words, if the moment problem does not permit any MED solutions, e.g. (Junk, 1998), the MESSY estimate of MED fails to infer any densities since there is no MED solution to be found. Secondly, there are guarantees that the cocfficient of the leading term in the exponent of the MED will be negative. Hence, in unbounded domains the MESSY estimate may diverge in the tails of the distribution. Both of these issues can be resolved by regularizing the MED ansatz, e.g. via incorporating a Wasserstein distance from a prior distribution as suggested by Sadr et al. (2023).
The most important, perhaps, distinguishing characteristic of the proposed methodology from non- parametric approaches, such as KDE, is the ability to recover a tractable symbolic estimator of the unknown density. This can be beneficial in number of applications, such as, for example, finding the underlying dynamics of a stochastic process. In such a case, where it is known that the transition probability takes an exponential form, this method can be used for recovering the drift and diffusion termss of an Ito process given a sequence of random variables, see c.g. Risken (1989); Oksondal (2013). Furthermore, the proposed method may be used to find a relation between moments of interest which could be helpful in determining the parameters of a statistical model, such as a closure for hierarchical moment description of fluid flow far from equilibrium. That said, exploring the efficiency of our algorithm in high-dimensional contexts is critical next step, considering the complexity of modern machine learning datasets. Other possible direc- tions for future work include: (i) data-driven discovery of governing dynamical laws from samples; and (ii) applications to variance reduction.
Acknowledgments
We would like to acknowledge Prof. Youssef Marzouk and Evan Massaro for their helpful comments and suggestions. M. Sadr thanks Hossein Gorji for his stimulating and constructive comments about the core idcas presented in this paper. T. Tohme was supported by Al Ahdab Fellowship. M. Sadr acknowledges funding provided by the German Research Foundation (DFG) under grant number SA 4199/1-1.
References
Rafail Abramov. A practical computational framework for the multidimensional moment-constrained maxi- mum entropy principle. Journal of Computational Physics, 211(1):198 200199999999999999 2006.
Rafail V Abramov. An improved algorithm for the multidimensional moment-cconstriined maximum entropy problem. Journal of Computational Physics, 226(1):621 644, 2007.
Rafail V Abramov. The multidimensional moment-constrained maximum entropy problem: A BFGS algo-
rithm with constraint scaling. Journal of Computational Physics, 228(1):06-108, 2009.
Rafail V Abramov. The multidimennsional maximum entropy moment problem: A review of numerical
methods. Communications in Mathematical Sciences, 8(2):377-392, 2010.
Naum Il'chh Akhiezer and N Kemmer. The classical moment problem: and some related questions in analysis, volume 5. Oliver & Boyd Edinburgh, 1965.
21

---

Published in Transactions on Machine Learning Research (01/2024)
Graham W Alldredge, Cory D Hauck, Dianne P O'Leary, and André L Tits. Adaptive change of basis in entropy-based moment closures for linear kinetic equations. Journal of Compulational Physics, 258: 489 508, 2014.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. In International Conference Machine Learning, 2021.
Zdravko Botev. Nouparametric density estimation via diffusion mixing. 2007.
Zdravko I Botev, Joseph F Grotowski, and Dirk P Kroese. Kernel density estimation via diffusion. The Annals of Statistics, 38(5), 2010.
Steven Brunton, Joshua Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113 (15):3932-385 39339337-377,,, 2016.
Bogdan Burlacu, Gabriel Kronberger, and Michael Kommenda. Operon C++: An efficient genelic program- ming framework for symbolic regression. In Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion, GECCO '20, pp. 1562-1570, New York, NY, USA, 2020. Association for Com- puling Machinery. ISBN 9781450371278 doi: 011145/3377929333333
Zhenning Cai, Yuwei Fan, and Ruo Li. A framework on moment model reduction for kinetic equation. SIAM Journal on Applied Mathematics, (5(55:2001-32
2015.
Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, and Philippe Rigollet. SVGD as a kernelized Wasserstein gradient flow of the chi-squared divergence. Advances in Newral Information Processing Sys- tems, 33:2098 2109, 2020.
Edwin KP Chong and Stanislaw H Zak. An introduction to optimization, volume 75. John Wiley & Sons, 2013.
Ronald R Coifman and Stéphanc Lafon. Diffusion maps. Applied and computational harmonic analysis, 21 (1):5-30, 2006.
Ronald Coifman, Stephane Lafon, Ann Lee, Mauro Maggioni, Boaz Nadler, Frederick Warner, and Steven W Zucker. Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps. Proceedings of the national academy of sciences, (02(21):4426-443) 2005.
Miles Cranmer, Alvaro Gonzalz,G Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. arXiv preprint arXiv:2006. 11287, 2020.
Kristian Debrabant, Giovanni Samacy, and Przemyslaw Zielinski. A micro-macro acceleration method for the Monte Carlo simulation of stochastic differential equations. SIAM Journal on Numerical Analysis, 55 (6):2745-2P))) 2786, 2017.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real NVP. arXie preprinl arXiv:1605.08e. 2016.
David L Donoho, lain M Johnstone, Gérard Kerkyacharian, and Dominique Picard. Density estimation by wavelet thresholding. The Annals of statistics, pp. 508-539, 1996.
Wolfgang Dreyer. Maximisation of the entropy in non-equilibrium. Journal of Physics A: Mathematical and General, 20(15);6505, 1987.
Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios. Neural spline flows. Advances in neural informalion processing systems, 32, 2019.
Anthony William Fairbank Edwards. Likelihood. CUP Archive, 1984.
22

---

Published in Transactions on Machine Learning Research (01/2024)
Aldo Frezzotti, Livio Gibelli, and Silvia Lorenzani. Mean field kinetic theory description of evaporation of fluid into vacuum. Physics of Fluids, 17(1):012102, 2005. doi: 10. 10.1063/1... 1063/1.1824111.
Alfredo Garbuno-nigo,, Franca Hoffmann, Wuchen Li, and Andrew M Stuart. Interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler. STAM Journal Applied Dynamical Systems, 19(1): 412 441, 2020a.
Alfredo Garbuno-Inigo Nikolas Nusken, and Sebastian Reich. Affine invariant interacting Langevin dynamics for Bayesian inference. SIAM Journal on Applied Dynamicol Systems, 19(3):1633 1658, 2020b.
Luc Giraud, Julien Langou, and Miroslav Rozlozmı. On the round-off error analysis of the gram-schmidt algo- rithm with reorthooonalizatin.. Technical report, Technical Report TRR/PA/02/34 CERFACS, Toulouse, France, 2002.
Harold Grad. Note on N-dimensional Hermite polynomials. Communications on Pure and Applied Molhe- matics, 3(4))325-33 330, 1949.
Cory D Hauck, C David Levermore, and André I, Tits. Convex duality and entropy-based moment closures: Characterizing degenerate densities. SIAM Journal on Control and Optimization, 47(4):1977-2000 2008,
M. Hermite. Sur nouveau déucloppement en série des fonctions. Imprimerie de Gauthier- 1864.
Ying Jin, Weilin Fu, Jian Kang, Jiadong Guo, and Jian Guo. Bayesian symbolic regression. arXiv preprint arXiv:1910. 08892, 2019.
M Chris Jones, James S Marron, and Simon J Sheather. A brief survey of bandwidth selection for density ostimation. Jowrnal of the American statistical association, 91 (433)::01-4077 1996.
Michael Junk. Domain of definition of levermore's five-moment system. Journal of Statistical Physics, 93: 1143-1167, 1998.
Zahra Kadkhodaie, Florentin Guth, Stéphane Mallat, and Eero P Simoncelli. Learning multi-scale local conditional probability models of images. arXiv preprint arXiv:2303.02m 02984, 2023.
Fierre-Alexanee Kamienny, Stéphane A'scoll,, Guillaume Lample, and François Charton. End-to-end sym- bolic regression with transformers. arXiv preprint arXiv:2204.1055 10532, 2022.
Jagat Narain Kapur. Mazimum--nntropy models in science and engineering. John Wiley & Sons, 1989.
A Ya Khinchin. Mathematical foundations of information theory. Courier Corporation, 2013.
Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1xl convolutions. Advances neural information processing systems, 31, 2018.
Ivan Kobyzev, Simon JD Prince, and Marcus A Brubaker. Normalizing flows: An introduction and review current mothods. IEEE transactions on pattern analysis and machine intelligence, 43911):394,,3979 2020.
Michael Kommenda, Bogdan Burlacu, Gabriel Kronberger, and Michael Affenzeller. Parameter identification for symbolic regression using nonlinear least squares. Genelic Programming and Evolvable Machines, 21 (3):4715501, 2020.
Misaki Kon, Kazumichi Kobayashi, and Masao Watanabe. Method of determining kinetic boundary condi- tions in net evaporation/ condensation. Physics of Fluids, 26(7):072003, 2014.
John R Koza and John R Koza. Genetic programméng: the programming of computers by means of natural selection, volume 1. MIT press, 1992.
William La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabrício Olivetti de França, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason H Moore. Contemporary symbolic regression methods and their relative performance. Xix preprint orXiv:2107 14351, 2021.
23

---

Published in Transactions on Machine Learning Research (01/2024)
Peter M Lee. Bayesian statistics. Arnold Publication, 1997.
Thomas Leonard and John SJ Hsu. Bayesian methods: an analysis for statisticians and interdisciplinary researchers, volume 5. Cambridge University Press, 2001.
C David Levermore. Moment closure hierarchies for kinetic theories. Journal of statistioul Physics, 83(5-6): 1021-1065, 1996.
Fengyi Li and Youssef Marzouk. Diffusion map particle systems for generative modeling. arXiv preprint arXiv:2344000000 2023.
Qiang Liu. Stein variational gradient descent as gradient flow. Advances neural information processing systems, 30, 2017.
Qiang Liu and Dilin Wang- Stein variational gradient descent: A general purpose bayesian inference algo- rithm. Advances in neural information processing systems, 29, 2016.
Qiang Liu and Dilin Wang. Stein variational gradient descent as moment matching. Advances in Newral Information Processing Systems, 31, 2018.
Stéphane Mallat. wavelet tour of signal processing. Elsevier, 1999.
Tanguy Marchand, Misaki Ozawa, Giulio Biroli, and Stéphane Mallat. Wavelet conditional renormalization group. arXiv preprint arXiv:2207.04941 2022.
James McDonald and Manuel Torrilhon. Affordable robust moment closures for CFD based on the maximum- entropy hierarchy. Journal of Computational Physics, 251:500-523, 2013.
T Nathan Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio P Santiago, Daniel M Faissol, and Brenden K Petersen. Symbolic regression via neural-guided genetic programming population seeding. arXiv preprinl. arXiv:2111.00053 2021.
Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. Science, 365(6447) 2019.
Bernt Oksendal. Stochastic differential equations: an introduction with applications. Springer Science & Business Modia, 2013.
Patryk Orzechowski, William La Cava, and Jason H Moore. Where are now? large benchmark study of recent symbolic regression methods. In Proceedings of the Genelic and Evolutionary Compulation Conference, pp. 1183- 1190, 2018.
Claudia Patrignani, K Agashe, G Aielli, C Amsler, M Antonelli, DM Asner, H Baer, Sw Banerjee, RM Bar- nett, T Basaglia, et al. Review of particle physics. CHINESE PHYSICS C, 2016.
Yudi Pawitan. In all likelihood: statistical modelling and inference using likelihood. Oxford University Press, 2001.
Brenden K Petersen, Mikel Landajuela Larma, Terrell N. Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, 2021.
Eckhard Platen and Nicola Bruti-Liberati. Numerical solution of stochastic differential equalions with jumps in finance, volume 64. Springer Science & Business Media, 2010.
William A Porteous, M Paul Laiu, and Cory Hauck. Data-driven, structure-presertte approximations entropy-based moment closures for kinetic equations. arXiv preprint arXiv:2106.08977 2021.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International con-
ference on machine learning, pp. 1530 1538. PMLR, 2015.
24

---

Published in Transactions on Machine Learning Research (01/2024)
H. Risken. The Fokker-Plande Equation: Methods of solution and applications. Springer, 1989.
Murray Rosenblatt. Remarks on some nonparametric estimates of density function. The onnals of mathe- matical statistics, pp. 832-68 832-837, 1956.
Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Data-driven discovery of partial differential equations. Science advances, 1(022446
602614, 2017.
Yousef Saad. Iterative methods for sparse linear systems, volume 82. SIAM, 2003.
Mohsen Sadr, Manuel Torrilhon, and M Hossein Gorji. Gaussian process regression for maximum entropy distribution. Journal of Computational Physies, 418: 109644, 2020.
Mohsen Sadr, Qian Wang, and M Hossein Gorji. Coupling kinetic and continuum using inta-driven maximum entropy distribution. Journal of Computational Physics, 444:1 10542, 2021.
Mohsen Sadr, Nicolas G Hadjiconstantinou, and M Hossein Gorji. Wasserstein-penat penalized entropy closure: A case for stochastic particle methods, arviv preprint ar Xrx::22008.00 02600, 2023.
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and control. In International Conference on Machine Learning, pp. 4442 4450. PMLR, 2018.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324 (5923):81-85, 2009.
Steffen Schotthöfer, Tianbai Xiao, Martin Frank, and Cory D Hauck. Neural network-based, structure- preserving entropy closures for the Boltzmann moment system. arXiv preprint arXiv:2201.1086, 2022.
Simon Sheather. Density estimation. Statistical science, pp. 588-591 2004.
Yang Song and Stefano Ermon. Improved techniques for training score- based generative models. Advances in neural information processing systems, 33:12438 -24488 2020.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXio preprint arXiv:2011.13450, 2020.
Aldo Tagliani. Hausdorff moment problem and maximum entropy: unified approach. Applied Mathematics and Computation, 105(2-3):291-30, 1999.
Tony Tohme, Kevin Vanslette, and Kamal Youcef-Toumi. Toumi. generalized Bayesian approach to model cali- bration. Reliability Engineering & System Safety, 204:107141, 2020.
Tony Tohme, Dehong Liu, and Kamal Youcef-Tounis GSR: A generalized symbolic regression approach.
Transactions on Machine Learning Research, 2023. ISSN 2835-8856. URL https://openrevee/ net/
forumPid=-heUUx-DNNP.
Belinda Tzen and Maxim Raginsky. Neural stochastic differential equations: Deep latent Gaussian models in the diffusion limit. arXiv preprint arXiv:190500955 2019.
Silviu-Marian Udrescu and Max Tegmark. Al Feynman: A physics-inspired method for symbolic regression. Science Advances, 6616), 2020.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Al Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Neural Information Processing Systems, 33:4860-488 4871, 2020.
Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. Symboliegpt: A generative transformer model for symbolic regression. orXiv preprint arXiv:2106. 14131, 2021.
Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2009.
25

---

Published in Transactions on Machine Learning Research (01/2024)
Marco Virgolin and Solon P Pissis. Symbolic regression is NP-hard. arXiv preprint arXtt:2207.01013 2022.
Sven Wang and Youssef Marzouk. On minimax density estimation via measure transport. ar preprint wrXiv:2277.1031 2022.
CJ Wild and GAF Seber. Nonlinear regression. New York: Wiley, 1989.
Hengzhe Zhang, Aimin Zhou, Hong Qian, and Hu Zhang. PS-Tree: A piecewise symbolic regression tree. Swarm and Evolutionary Computation, 71:101061, 2022.
26

---

Published in Transactions on Machine Learning Research (01/2024)
A Maximum entropy distribution function
The maximum entropy distribution (MED) function finds the least-biased closure for the moment problem. Given N, realizable moments E RR' associated with polynomial basis functions H of the unknown distribution function, MED obtained by minimizing the Shannon entropy with constraint moments using the method Lagrange multipliers as Hauck et al. (2008)
No Σ 1=1
((WL) mi(())
C[F()]] = (log (F(x)) 1)dx
A
(A.1)
By taking the variational derivative functional A.1, the extremum is found as
No F(Φ) 2 exp Σ (()
N2 Σ Addi(ae) dxc
where Z = exp
(A.2)
which is referred to as the maximum entropy distribution function. The Lagrange multipliers ap- pearing in Eq. (A.2) can be found using the Newton-Raphson approach. As formulated in (Debrabant et al., 2017), the raced dual formulation D(A) provides us with the gradient g (D()) and
Hessian [A) (() as
No
9i = ui 1u--//
7
1
11, exp
Σ
Aetli
de for i I,,,, No
(A.3)
No. dx for i,j = I,.., No Hg
1 H exp
and Laj 7
(A.4)
1=1
Once the gradient and Hessian are computed, the Lagrange multipliers X can be updated via
A A
(A.5)
detailed in Algorithm 4.
Algorithm 4: Newton's method for finding Lagrange multipliers of MED given moments μ for given
tolerance
Input: μ,, A
Initialize À + Ao:
Compute gradient 9 and Hessian L, i.c. Eqs. (A.3)-(A.4);
while ||g|| > do
Updale À + A L -
Update gradient g and Hessian L with the new À via numerical integration of Eqs. (A.3)-(A.4).
end
Return A;
27

---

Published in Transactions on Machine Learning Research (01/2024)
B Maximum cross-entropy distribution function
The maximum cross-cntropy distribution function (MxED) finds the least-biased closure for the moment problem given Ng realizable moments E IR associated with polynomial basis functions H of the unknown distribution function along with the prior J as the input. In other words, addition to the moments of the target distribution, in the MxED method we also have access to a prior distribution JPPrior as well as N samples of the former, ie. (xprion fj=1 - MxED is obtained by minimizing the Shannon cross-entropy from the prior with constraint moments using the method of Lagrange multipliers via the functional
F(x)) J
No
(HH()) - me(a))
c[F(x)] := F(x)log
dæ +
Σ Ae 1=1
(B.1)
By laking the variational derivative of functional B.1, the extremum found to be
Na
Na (xH dx.
F(z) = Z
Prior (x) exp
Σ
X A.H(())
where Z = Fprior (x) exp
Σ
(B.2)
Similar to the maximum entropy distribution function, the Lagrange multipliers appearing in Eq. (B.2) can found by following the Newton- Newton-Raphson approach. As formulated in (Debrabant et al., 2017), the uncon- strained dual formulation D(A) provides us with the gradient g VD(A) and Hessian L(A) = D2D as
No Σ
1 /
= #4 Z
P Hi exp
H
dx for i = 1,..,. No
(B.3)
and Lij
1 2
prior H,H, exp
X Hk dæ for i,j 1,... No
(B.4)
Once the gradient and Hessian computed, the Lagrange multipliers X can be updated via
+ A
(B.5)
Since we have access to the samples of Prior FPrior, we use the given samples to compute the gradient and Hessian, i.c.
9i Rd bi - (Hi (x >w for i I,..,, No Lijj R (H, (H; Prior) H, (X Prlor)) ) for i,j 1,.,,,
(B.6)
(B.7)
where ((XPrl = exp (ΣN, \ H2 (xrroo)) denotes weights for calculating moments using importance sampling, L.e. (o(X)) > 5).
o W
W(X)). More details can be found below (Algorithm
Algorithm 5: Newton's method for finding Lagrange multipliers of MxED given moments a and samples of prior Prior ~ PPrior for a given tolerance
Input: μ, X Prior ~ .Prior
Initialize X + 0;
Compute gradient 9 and Hessian L, i.c. Eqs. (B.6)-(B.7): while [|g| > do
Update À + A L L-1 g;
Update gradient g and Hessian L with the new À using samples, i.c. Eq. (B.6)-(B.7);
end
Return A;
28

---

Published in Transactions on Machine Learning Research (01/2024)
C Ablation studies
In this section, we report on the results of an ablation study that shows the effects of cach component in the proposed MESSY solution. The four components of MESSY Estimation are: i) the symbolic component, ii) the multi-level component, iii) the orthogonalization, and iv) the cross-entropy step. As we have already compared MESSY-S and MESSY-P (MESSY without the symbolic regression component) throughout the paper, we conduct an ablation study on the remaining components, i.c. the multi-level, the orthogonalization, and the cross-entropy step.
C.1 Multi-level component of MESSY
First, order show the effect of the multi-level part of the MESSY algorithm, let us consider the bi-modal density near the limit of realizability, see 8.2. In particular, let us consider the MESSY-P estimate with 4th order polynomials as basis functions. In Fig. 12, show that by disabling the multi-level step, the estimated density covers only one of the peaks and entirely misses the other one. We believe that this is due to the high conditionality of the matrix ME, However, MESSY-P with multi-level step accurately recovers both peaks of the underlying density.
Samples MESSY-P multi-level)
MESSY (w/ multi-level)
f(a)
Figure 12: Ablation study the multi-level component of MESSY-P estimation for the case of target distribution near the limit of realizability.
C.2 Orthogonalization
Another important element of MESSY process is the orthogonalization basis functions given samples, i.c. the modified Gram-Schmidt Algorithm As an example, here we consider the MESSY-P estimate of larget bimodal distribution function presented in 8.1 test case and report the condition number of the matrix ME that needs to be inverted to compute the Lagrange multipliers in Eq. (18). As shown in Fig. 13, the condition number of ME increases with the polynomial order, making the linear system (18) stiff and difficult to solve. However, by orthogonalizing the basis functions, we can maintain low condition number for the outcome ME which makes the resulting linear system tractable.
C.3 Cross-entropy step
The final step in the MESSY-P procedure is the cross-entropy. The goal of this step is to reduce the bias that may have been introduced by the multi-level component of the algorithm. In the MESSY Algorithm 3, we take the outcome of the multi-level step as the prior, and correct Lagrange multipliers to match the target moments. In Fig. 14, we compare the MESSY-P density estimate of the discontinuous density from the test case presented 8.3 with and without cross-cntropy step. As expected, we observe that the cross-cntropy
29

---

Published in Transactions on Machine Learning Research (01/2024)
MESSY-P (w/o orthogonalization)
1044
MESSY-P (w/ orthogonalization)
Condition Number
107
10a
101
12
14
Figure 13: Ablation study the orthogonalization component of MESSY-P algorithm for the case of bimodal distribution. Relation between condition number of matrix L and the order of considered polynomial with and without the orthogonalization step.
step is essential in recovering distributions with discontimuities. This due to the fact that the solution obtained via gradient flow assumes a smooth target density. Hence, it fails to find the Lagrange multipliers for the target discontimuous pdf. However, the cross entropy step only assumes that the target distribution is integrable to the extent that the given moments exist. This weaker condition is essential to recovering discontinuous distributions, egg. the one showed in Fig. 14.
MESSY-P (w/ cress-entrepy)
MESSY (w/o cress entropy)
f(2)
0.4
0.2
0.0
-2
Figure 14: Ablation study on the cross-entropy step of MESSY-P estimation of a discontinuous density.
D Histogram estimate to the bi-modal distribution function
Figure 15 provides comparison between histogram, KDE, MESSY-P and MESSY-S estimates of the underlying bi-modal distribution function considered in 8.1 given N 100, 1000, 10000 samples. Clearly, MED estimates provide a better solution than non-parametric ones, i.e. histogram and KDE, when not many samples are available. However, as more samples are considered, the non-parametric estimators become more accurate than MESSY (or any parametric estimator) with fixed number of basis functions.
30

---

Published in Transactions on Machine Learning Research (01/2024)
Thue RDF
0.7
Tue PDF
Tue PDP
Histogram
Histogram
Hizoogrm
0.8
0.6
HxED
MiED MESSY-P
0.4
0.4
0.4
0.1
0.1
9.1
(a) 100 samples
(b) 1,000 samples
(c) 10,000 samples
Figure 15: Density estimation using KDE, MxED, MESSY-P, MESSY-S and hislogram given (a) 100, (b)
1,000, and (c) 10,000 samples.
E
Solution found by MESSY for the considered test cases
Table 1: Density expressions recovered by our MESSY estimation method for several distributions.
Example
Expression
Binodal
MESSY-P
f(x) 0,288e 01:72 +0.1002" n0040 -0002" "o..00 20002/11.110000 -3..1722 -1722 sourts +2.0012 -.022
MESSY-8
[1x) 0.913e 18se
(5.00))
Limit Realizabiiity
MESSY-P
f(x)
1.91-
10-1110------- 12.8702 - -40 oa swwooms 022417
10 T.000e 1200022-2200000 22200222 -4.022" -20,5000
MESSY-8
[(0)
4.134
10%
21.855 sn cottee " (irr 2.0
Sts Del boy 3000 (231200021 (oe
1.000
0.000 2009
Discontiinous
MESSY-P
f(r)
otherwise
MESSY-S
f(x)
0.203e "1. po.m
022440 (Sk) 0fficos
If
otherwise
Gaussion 2D
MESSY-P
((n). = 0.18e -00000 40.5 t0000 0705
MESSY-8
f((() 0.182 -ww.:
[d.477e m05x700 Oiteee aodie -0..0-0000000mmmm
Coramn-oxp
MESSY-P
f(s()2g)
otherwise
[0.017e-A11 p.000em 1..1 71 (1 s0tt
(2.8kx)
MESSY-8
otherwise
31


Reviews:
[{"review_id": "xSiwXNMdWG", "reviewer": "Reviewer_tPUb", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The proposed approach appears to be significantly novel. \n- The paper provides justification for different components of the approach. For instance, to illustrate the benefit of symbolic regression, a comparision is made between MESSY-P and MESSY-S, algorithms utilizing a polynomial basis and a basis obtained using symbolic regression respectively.\n- The paper discusses different practical aspects of the approach such as the computational cost and conditioning of optimization problems.\n- The non-technical sections of the paper are well-written.\n\nWeaknesses:\n\n- Some terms in the paper are not clearly defined. For instance, in the line above equation 8 \"Here we note that the Hessian for this optimization\", the \"optimization\" problem isn't clearly defined or explained. Only after Proposition 4.3 and Eq. 16, it becomes apparent that \"g(t)\" is being interpreted as a gradient for an associated optimization problem on the Lagrange multipliers. This should be explained clearly with precise definitions before statements such as \"These relaxation rates can be used as the gradient in the search for parameters of a given ansatz\" above eq. 7.\n- Certain references to the literature are innacurate. The paper refers to Eq.1 as the \"score-based generative model\" in Song et al. (2020). However in the reverse-SDE utilized in Song et al. (2020), the \"score\" is the grad-log of the density at the given time t, not that of a fixed density $\\hat{f}$ corresponding to the stationary distribution. Instead, Eq.1 should be more appropriately referred to as the overdampled Langevin dynamics for $\\hat{f}$ .\n- The approach lacks bounds on sample complexity (with a fixed set of basis functions). In contrast, the sample rates for existing approaches such as Kernel density estimation (KDE) have been extensively analyzed.\n- The approach is evaluated only on a few simple one-dimensional distributions. This makes it difficult to assess the significance of the improvements.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a novel algorithm for density estimation. The approach considers the evolution of the moments of certain basis functions when a density follows a diffusion process (overdamped Langevin dynamics) having the estimated density function as its stationary distribution. This rate of evolution is used to construct an objective that allows estimation of the Lagrange multipliers corresponding to the maximum entropy distribution w.r.t a set of basis functions. The estimated multipliers depend on certain expectations w.r.t the empirical measure of samples from the true density. The objective relies on the observation that the rate of evolution of the expectations w.r.t the dynamics at the true density is $0$ when the estimated density matches the true density. The approach further utilizes symbolic regression to learn suitable basis functions. The paper discusses computational aspects of the proposed approach and verifies its effectiveness on certain one-dimensional densities.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Adjustments critical for securing recommendation for acceptance:\n\n- It should be clarified that Eq.1 is the overdampled Langevin dynamics corresponding to the density $\\hat{f}$ and differs from the score-based generative model in Song et al. (2020).\n- The optimization problem corresponding to $g(t)$ should be clearly defined and motivated.\n\nAdjustments that will strengthen the work:\n\n- Experiments on multi-dimensional distributions.\n- Discussion of sample and time complexity bounds for the proposed and existing approaches for density estimation.", "broader_impact_concerns": "The paper is primarily of a theoretical nature with simple toy experiments. Therefore, there are no concerns regarding its broader impact.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "e5xn7YNZBW", "reviewer": "Reviewer_h7mG", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\n- Despite the great amount of activity around symbolic regression in scientific machine learning over the last ten or so years, there are very few studies that address symbolic regression of probability density functions from finite samples, arguably a more difficult problem than \"simply\" fitting time series data. Given the central importance of density estimation in machine learning such efforts are clearly commendable.\n\n- Relying on the maximum entropy principle in density estimation, especially with explicit symbolic expressions, makes good sense. It is also a nice idea to look for moments which result in well-conditioned linear system for the Lagrange multipliers. This strategy results in an effective model in the reported computer experiments on simple bimodal univariate distributions. \n\n- In general, it is exciting to think about how symbolic regression could be used to facilitate the usually complicated problems of connecting score with density and in general work with the structured distributions we encounter in scientific ML.\n\n## Weaknesses\nI feel that the organization and motivation of the manuscript could be improved. More importantly, the computer experiments are lacking in many aspects.\n\n- Definition 5.1 is not directly used in the paper. Why not start directly from Definition 5.2? Perhaps there is actually no need to state these problems as formal mathematical definitions? Moreover, Definition 5.1 defines supervised learning, not symbolic regression, unless it says much more about how \\mathcal{S} is structured.\n\n- With the exception of Table I in the appendix there is no discussion about the used space of symbolic expressions. From Table I it seems that all expressions consist of compositions and products of sin, cos, and monomials. It is not clear what the advantage of generating symbolic expression at random is as opposed to, say, random Fourier features or random neural networks. The main reason to use things other than monomials appears to be the improved conditioning of the linear system for the multipliers and the expressions listed in Table I are hardly interpretable. Since symbolic regression is central to the proposed method this lack of information is strange.\n\n- **(Major) Mismatch Between Claims and Evidence**: The manuscript emphasizes its relevance in high-dimensional contexts. Quoting from page 5, \"With applications to high-dimensional problems in mind, instead of looking for solutions of Eq. 2, we choose to work with appropriate empirical moemnts of this equation, ...\" Yet the computer experiments focus solely on very simple one-dimensional scenarios. To make a strong case for high-dimensional applicability, it is crucial to provide experiments or discussions specific to these contexts. I imagine that this may work in dimensions 2 or 3 (although already with difficulties; see similar issues in symbolic regression methods like SINDy), but it would face serious issues going beyond that, let alone to what is usually considered high-dimensional in modern machine learning.\n\n- There are additional major concerns with computer experiments. Although the authors deal with various operators (e.g., +, −, ×, ÷) and unary functions (e.g., cos, sin, exp, log), the computer experiments only estimate a simple bimodal distribution, with a very simple (quadratic or polynomial) exponent. \n\n- Can the authors compare their results directly with histogram density? In Figure 2(c), I wonder how the histogram density looks with 10,000 samples. I suspect it may not be very far from the true pdf.\n\n- On a related note: I have some doubts about the comparison with KDE. It seems that Silverman's kernel bandwidth systematically overestimates the optimal bandwidth. How about using an oracle estimate for the bandwidth? Or a cross-validated / bootstrap estimate from the data? I believe that it would result in an excellent fit. Given how complex the proposed method is, a small search for the optimal bandwidth seems trivial.\n\n- The work should be properly situated in the literature and the setting should be made precise and explicit early on. Since you only work in 1D and only sample formulas randomly, the relevance for contemporary machine learning is not evident.\n\n- There is a huge amount of work on density estimation, especially in the high-dimensional setting, which seems highly relevant to this manuscript yet is not mentioned. An example that comes to mind is the recent work of Stèphane Mallat and collaborators on multiscale approaches inspired by the renormalization group. These in particular seem related to the multilevel scheme proposed in the manuscript, though rigorously motivated.\n\n- Many ideas in the first part of the paper come from Liu 2017. It is not clear for example why the proof of Proposition 1 is reproduced (but not in full). At any rate, one should acknowledge in the statement of the proposition or before that this is taken from earlier work.\n\n## Minor\n\n- I find the terminology \"moment conservation\" a bit strange in this context. I am used to seeing \"moment matching\" / \"method of moments\" (e.g., https://en.wikipedia.org/wiki/Method_of_moments_(statistics) and http://proceedings.mlr.press/v37/li15.pdf)\n\n\n- Please also specify what you mean by a \"least biased distribution\". I understand the relation to maximum entropy but if you are really referring to bias that should be properly defined.\n- After equation 4: \\bf{x}-> \\infty should be $\\| \\bf{x} \\|$\n- In the appendix, H should be \\bf{H}?\n- In the paragraph before equation (10), should the entropy be minimized or maximized?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a method for density estimation based on a combination of Langevin dynamics, maximum-entropy principle, moment matching and symbolic regression. They first define an exponential family ansatz motivated by entropy maximization. They then show how to estimate its parameters by computing relaxation rates. In the second part, they propose a Monte Carlo symbolic regression method to find closed-form moments beyond the usual monomials which result in a better conditioned linear systems for the related Lagrange multipliers. They apply the resulting method to simple bimodal univarate distributions and report good results.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "All of the following are critical in my opinion:\n\n- better organization of the manuscript, in particular with regard to positioning in the literature and a precise statement of the adopted setting, assumptions, notation, metrics, etc, early on in the manuscript\n- a better comparison with state of the art density estimation\n- comparisons with KDE with optimal data-driven kernel bandwidth choice\n- a comparison with a histogram estimate with an optimal bin size\n- (very critical) experiments on high-dimensional data\n- ablations to test the importance of the different design choices\n- more details about the used symbolic regression and the rationale for the related parameter choices\n- many more ground truth densities in various dimensions, including those coming from applications, and including real rather than synthetic datasets", "broader_impact_concerns": "--", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "02WTZixudn", "reviewer": "Reviewer_3RF5", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nThe manuscript elegantly introduces the idea of using gradient flow relaxation rates as an objective for optimizing a representation of the target density. In large part, the exposition is clear and not overly technical, while still mostly precise. Extending the density estimation problem to symbolic regression specialize the scope of the method, which in turn perhaps narrows the number of methods to which it should be compared. \n\nWeaknesses:\n \nConnections with the existing literature on density estimation, especially variational methods using the kernelized Stein discrepancy, are not really sufficient to determine the novelty of this approach beyond the symbolic regression problem.\n\nSome of the claims are, I think, overstated or made without a clear argument. Especially in Sec. 4.1.\n\nWhen are the relaxed existence requirements relevant?\n\nWhat does it mean to avoid the curse of dimensionality in this setting?\n\nFinally, I think the appropriate set of problems for this method could be more clearly highlighted. Much of the current ML literature on density estimation uses either nonparameteric density estimation to obtain convergence bounds or parametric estimators based on neural networks. If neither of these approaches are appropriate for the class of tasks that MESSY is intended to apply to, that could also be highlighted more cleary.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The central idea in this work is to draw iid samples from an unknown distribution $f$ and subsequently measure the distance between an ansatz $\\hat f$ and $f$ by evolving these samples using an SDE where the drift is the gradient of the log likelihood of $\\hat f$. Because an ensemble of particles evolving under this SDE converges weakly to the solution of an associated Fokker Planck equation, the relaxation time under this dynamics provides a measure of the quality of the estimator $\\hat f$. By introducing basis vectors and solving the linear problem when the relaxation rate is set to its limiting value (zero), the authors obtain a solution for $\\hat f$ without optimization. This procedure is subsequently used for symbolic regression. Numerical experiments on a few 1d examples are conducted and largely illustrate that MESSY obtains better solutions than convention KDE.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The related methods section should discuss Stein Variational gradient descent with a particular emphasis on moment matching variants of SVGD. \n\nSome discussion of sample problems for which SR would be natural would help motivate the objective of the paper.\n\nThe claims in Sec. 4.1 should be clearly justified. For example, I highly skeptical of the claim that the proposed method avoids the curse of dimensionality in any meaningful sense. This statement requires a clear argument, such as a convergence rate bound. \n\nFinally, I do not think the numerical experiments are sufficient to support the claims made in the paper. It is imperative that the authors consider a high-dimensional example. Furthermore, direct comparison to more modern methods (e.g., SVGD) should be made. If I understand correctly, the symbolic regression task could be implemented once any estimator is obtained.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper presents an algorithm for density estimation using overdamped Langevin dynamics and symbolic regression. This method, which focuses on the evolution of moments to determine Lagrange multipliers for maximum entropy distribution, assesses accuracy based on the alignment with true density dynamics. While the algorithm shows accuracy and introduces new ideas, there is uncertainty about its effectiveness in high-dimensional data, as the current research is limited to simple one-dimensional cases despite claims of potential high-dimensional applicability.
Decision: Accept (poster)