SymDE: Symbolic Density Estimation

• SymDE is a family of methods that extract explicit probability laws from samples by generating human-readable, closed-form expressions for densities or PMFs.
• It leverages techniques like orthogonal expansions, maximum-entropy formulations, and invertible symbolic flows to derive interpretable models with valid probabilistic structure.
• The approach bridges traditional density estimation and symbolic regression, employing decomposition and symmetry exploitation to tackle high-dimensional challenges.

Symbolic Density Estimation (SymDE) denotes a family of methods that recover explicit, human-readable probability laws from samples, typically as closed-form densities, probability mass functions, or analytically specified invertible maps rather than purely numerical kernel sums or large black-box networks. In recent formulations, the task is stated as recovering an interpretable, closed-form expression for an unknown probability density from samples, or automatically discovering a closed-form PMF from observed discrete data without prespecifying the family (Rajendram et al., 30 Mar 2026, Liu et al., 20 May 2026). Across the literature, the symbolic object varies: orthogonal expansions with coefficients given by sample averages, maximum-entropy densities with symbolically chosen basis functions, symbolic normalizing flows, and log-PMF expressions constrained by validity-aware search (Duda, 2017, Tohme et al., 2023, Tohme et al., 2024, Liu et al., 20 May 2026).

1. Scope, definition, and boundaries

SymDE is defined by two coupled requirements: the estimated law must be probabilistic, and the representation must remain explicit enough to be inspected as a mathematical expression. In continuous settings, this includes symbolic formulas for densities, symbolic surrogates of nonparametric estimates, and invertible analytic transports to a reference distribution (Duda, 2017, Tohme et al., 2024, Rajendram et al., 30 Mar 2026). In discrete settings, it includes direct discovery of closed-form log-PMFs using operator libraries tailored to combinatorial probability laws (Liu et al., 20 May 2026).

A recurring distinction in the literature is between symbolic structure and symmetry-aware structure. The Rao-Blackwell estimator for marginal probabilities in large symmetric probabilistic models exploits automorphism groups and orbit-conditioned probabilities to reduce MSE, but it is not primarily a full symbolic density learner (Niepert, 2013). Homothetic log-concave estimation exploits geometric symmetry by representing a multivariate density as

$f(x)=\exp\{\phi(\|x-\mu\|_K)\},$

thereby reducing a high-dimensional problem to estimation of a one-dimensional radial generator and a shape parameter $K$, but the output is described as a shape-constrained nonparametric MLE rather than an explicit symbolic formula per se (Xu et al., 2019). Lie PCA similarly learns a Lie algebra of manifold symmetries and uses it for symmetry-aware augmentation, not for direct algebraic expression discovery (Cahill et al., 2020).

A related but distinct line appears in density functional theory. SyFES searches over human-readable symbolic expressions for exchange-correlation functionals in Kohn-Sham DFT, representing a functional as a sequence of executable instructions and optimizing structure by regularized evolution (Ma et al., 2022). This is symbolic function learning over density-dependent scientific objects, but not ordinary probability density estimation.

2. Analytic expansions and maximum-entropy formulations

A direct antecedent of SymDE is the closed-form orthogonal-expansion method of rapid parametric density estimation. The core setup estimates a density $\rho_a:\mathbb{R}^D\to\mathbb{R}$ from samples $S=\{x^i\}_{i=1}^n$, optionally with weights $W(x)$ and fitting weight $w(x)$. After introducing a KDE-smoothed target and taking the zero-width kernel limit in an $L^2$ objective, the paper obtains

$$\min_a \ \langle \rho_a,\rho_a\rangle - \frac{2}{n}\sum_{x\in S} w(x)W(x)\,\rho_a(x),$$

with

$$\langle f,g\rangle := \int_{\mathbb{R}^D} f(x)\,g(x)\,w(x)\,dx.$$

For linear models

$\rho_a=\sum_{i=1}^m a_i f_i,$

the coefficients satisfy a linear system, and for an orthonormal basis they decouple completely:

$K$0

The resulting symbolic densities can be written as compact polynomials, Fourier series, or orthogonal-function expansions with no iterative optimization. The paper works out Legendre-polynomial models on finite domains, Fourier-series models on intervals or tori, and Hermite-polynomial models with a Gaussian envelope on $K$1. Normalization can be built in by fixing the zero-order coefficient, or corrected by a Lagrange-multiplier formula. Under $K$2, an orthonormal basis, and a realizable true density, the coefficient error is stated to drop like $K$3 (Duda, 2017).

The same paper also extends the framework beyond ordinary probability densities. Because the averaging formulas allow arbitrary weights $K$4, the estimator can represent mass distributions or other non-probabilistic densities, while negative or complex weights can be used for clustering and class separation. In the complex-weight construction, one fitted analytic function can encode multiple classes, with the zero set, sign, or phase acting as the classifier boundary (Duda, 2017).

MESSY estimation adopts a different symbolic route. It assumes a maximum-entropy ansatz

$K$5

where $K$6 is a vector of basis functions and $K$7 are Lagrange multipliers. The paper motivates the method through a drift-diffusion process

$K$8

whose associated Fokker–Planck equation implies monotone decrease of a cross-entropy-like divergence. Applying moment operators to this dynamics yields a linear system for $K$9, so the parameters of the maximum-entropy density are obtained from sample moments rather than from nonlinear dual optimization. MESSY-S then uses symbolic regression to search over basis functions $\rho_a:\mathbb{R}^D\to\mathbb{R}$0, accepting candidates that lead to well-conditioned linear systems. A multi-level recursive decomposition, bounded-support handling, and a maximum cross-entropy correction are used for multimodal and discontinuous densities (Tohme et al., 2023).

These formulations define several recurrent SymDE templates.

Framework	Symbolic object	Main estimation route
Rapid parametric density estimation	Orthogonal expansion	Coefficient averaging
MESSY	Maximum-entropy density	Moment linear system
ISR	Invertible symbolic map	Flow likelihood
AI-Kolmogorov	Symbolic density expression	Surrogate-density regression
SDE for discrete distributions	Log-PMF expression	Evolutionary WLS

3. Symbolic flows and decompositional continuous-density pipelines

ISR extends normalizing flows into a symbolic regime by replacing the subnetworks in affine coupling blocks with EQL subnetworks. For density estimation, the model learns a bijection

$\rho_a:\mathbb{R}^D\to\mathbb{R}$1

and uses the change-of-variables formula

$\rho_a:\mathbb{R}^D\to\mathbb{R}$2

Because each reversible block is an affine coupling construction, the Jacobian is triangular and the determinant is tractable. The scaling and translation functions are symbolic because they are implemented by EQL networks built from operators such as $\rho_a:\mathbb{R}^D\to\mathbb{R}$3, with the broader ISR framework also mentioning $\rho_a:\mathbb{R}^D\to\mathbb{R}$4. The density-estimation objective is the standard normalizing-flow NLL with a Gaussian latent. On four 2D target distributions—Gaussian, Banana, Mixture of Gaussians, and Ring—the paper states that ISR finds and generates samples with slightly better accuracy than INN, and for the Gaussian example it recovers a near-exact symbolic map (Tohme et al., 2024).

AI-Kolmogorov takes a decompositional approach. Its pipeline is explicitly multi-stage: problem decomposition through clustering and/or probabilistic graphical model structure learning, nonparametric density estimation, support estimation, and symbolic regression on the resulting density estimate. Additive decomposition is handled by clustering, with DBSCAN used as an example when the density is multimodal; multiplicative decomposition is handled by structure learning, with the PC algorithm used to identify conditional independences. The surrogate density target is supplied by KDE or, when KDE is inadequate, Neural Spline Flows. Support is estimated either through level sets

$\rho_a:\mathbb{R}^D\to\mathbb{R}$5

or by geometric heuristics such as convex hulls. PySR is then applied to regress a symbolic expression against the surrogate density on a grid over the support (Rajendram et al., 30 Mar 2026).

The reported examples span a 2D Gaussian mixture, a 4D Gaussian with independent blocks, a Rastrigin-based density, muon decay, and a dijet distribution. Decomposition materially changes the symbolic difficulty: for the Gaussian mixture, clustering improves MSE from $\rho_a:\mathbb{R}^D\to\mathbb{R}$6 to $\rho_a:\mathbb{R}^D\to\mathbb{R}$7, and for the 4D Gaussian, structure learning improves MSE from $\rho_a:\mathbb{R}^D\to\mathbb{R}$8 to $\rho_a:\mathbb{R}^D\to\mathbb{R}$9 (Rajendram et al., 30 Mar 2026). At the same time, the paper is explicit that exact recovery is not guaranteed, that symbolic expressions can exhibit small negative regions, and that sharp support boundaries remain hard.

4. Discrete distributions and validity-aware PMF discovery

The discrete-distribution variant of SymDE works directly in the log-PMF domain. Given samples $S=\{x^i\}_{i=1}^n$0, a grid $S=\{x^i\}_{i=1}^n$1, counts $S=\{x^i\}_{i=1}^n$2, and Laplace or symmetric Dirichlet smoothing, the method defines

$S=\{x^i\}_{i=1}^n$3

Because $S=\{x^i\}_{i=1}^n$4 is heteroscedastic, the search objective is weighted least squares,

$S=\{x^i\}_{i=1}^n$5

with weights derived from a Delta-method approximation to the variance of $S=\{x^i\}_{i=1}^n$6. The operator library is PMF-specific and includes arithmetic, elementary functions, and log-domain combinatorial primitives such as

$S=\{x^i\}_{i=1}^n$7

To control the search space, the method adds operator complexity profiles, grammar constraints restricting arguments of special functions to affine forms, and atomicity constraints on exponents. Search is population-based and evolutionary, followed by a validity-aware inference stage that screens candidates by loss, approximate normalization, bounded log mass, and complexity budget (Liu et al., 20 May 2026).

The same framework extends naturally to richer families. Finite mixtures are represented in log-space with

$S=\{x^i\}_{i=1}^n$8

while zero-inflated models use $S=\{x^i\}_{i=1}^n$9 as a primitive for the point mass at zero. The benchmark SDEBench contains 14 base discrete distributions across exponential-family, heavy-tailed, mixed, and additional structured families, together with binomial mixtures and zero-inflated models. The paper reports recovery of the correct symbolic form for all 14 base distributions, correct nested mixture structure for binomial mixtures, and correct zero-inflated forms for ZIP, ZINB, and ZIG, with parameter estimates typically within a few percent of the true values (Liu et al., 20 May 2026).

The real-data PBMC scRNA-seq experiment illustrates the interpretive aim of discrete SymDE. For Gene 4046 in PBMC3k, the discovered expression is

$W(x)$0

a nested logaddexp expression that the paper interprets as a finite-mixture-style model. On empirical log-PMF MSE, SDE reports $W(x)$1, compared with $W(x)$2 for PySR, $W(x)$3 for Pyro, $W(x)$4 for KDE, $W(x)$5 for MLE+AIC over 10 standard families, and $W(x)$6 for MoM (Liu et al., 20 May 2026).

5. Density-guided symbolization, symmetry exploitation, and adjacent extensions

A substantial adjacent literature uses density estimation to drive symbolic representations rather than to output a closed-form density directly. In time series, pSAX replaces SAX’s Gaussian/equiprobable discretization by a pipeline of PAA, KDE, and Lloyd-Max quantization, so that symbol intervals minimize expected distortion under an estimated density. cSAX uses PAA, KDE, and Mean-Shift clustering, allowing the alphabet size to be determined by density modes. The paper also proves that if PAA is applied to Z-normalized but correlated Gaussian samples, the resulting variance is

$W(x)$7

which is typically smaller than 1, and recommends normalizing after PAA rather than before it. Both pSAX and cSAX are reported to preserve the attractive lower-bounding properties of SAX while improving tightness, RMSE, or anomaly-detection performance in real-world datasets (Bountrogiannis et al., 2021).

edwSAX follows a closely related logic. It estimates a KDE for the transformed time series values, chooses breakpoints $W(x)$8 so that each interval has probability mass $W(x)$9, and computes density-aware centroids by splitting each interval into equal-probability halves. The discretization rule remains SAX-like,

$w(x)$0

but breakpoints and reconstruction values are learned from the estimated density rather than fixed by Gaussian quantiles. On 20 UEA/UCR datasets, the paper reports average reconstruction RMSE improvements from $w(x)$1 to $w(x)$2 at alphabet size 5 and from $w(x)$3 to $w(x)$4 at alphabet size 10, with practical lower-bound tightness tradeoffs around 20–30 symbols and diminishing gains after about 40 (Kloska et al., 2022).

Symmetry-aware estimation supplies a different structural route. In large probabilistic models, the Rao-Blackwell estimator

$w(x)$5

uses automorphism-group orbits to lower MSE relative to the standard sample-frequency estimator and is linear-time computable for single-variable marginals (Niepert, 2013). For high-dimensional continuous densities, homothetic log-concave estimation attains a worst-case $w(x)$6 risk bound independent of $w(x)$7 when $w(x)$8 is known, precisely because the symmetry reduction turns the statistical difficulty into a one-dimensional radial problem (Xu et al., 2019). On manifolds, Lie PCA estimates a symmetry algebra $w(x)$9 from local tangent information and generates synthetic points by

$L^2$0

improving density estimation relative to KDE and local PCA in several symmetric-manifold examples (Cahill et al., 2020). These methods are not symbolic in the narrow algebraic sense, but they embody the same principle of exploiting structure to compress density estimation.

Beyond probability, symbolic search over density-related objects appears in scientific computing. SyFES represents semilocal exchange-correlation functionals as short instruction sequences, uses regularized evolution with CMA-ES parameter optimization, reconstructs a B97-equivalent form from scratch on the TCE subset of MGCDB84, and evolves GAS22 from $L^2$1B97M-V with a reported 15% improvement in overall test WRMSD on the MGCDB84 test set (Ma et al., 2022). A plausible implication is that domain-specific operator libraries and structural priors, emphasized in discrete SymDE as well, are central whenever symbolic density-like objects are sought in scientific models.

6. Limitations, misconceptions, and open directions

A common misconception is that SymDE is simply generic symbolic regression applied to sample clouds. The literature rejects that equivalence. Standard symbolic regression assumes labeled targets, while density estimation starts from unlabeled samples; generic symbolic regression also does not know about PMF validity, non-negativity, normalization, or combinatorial structure (Rajendram et al., 30 Mar 2026, Liu et al., 20 May 2026). As a result, successful SymDE systems typically inject probabilistic structure explicitly: maximum-entropy ansätze in MESSY, tractable Jacobians and change-of-variables likelihoods in ISR, or validity-aware PMF screening in discrete SDE (Tohme et al., 2023, Tohme et al., 2024, Liu et al., 20 May 2026).

A second misconception is that symbolic output automatically guarantees a valid probability law. The early orthogonal-expansion approach notes that least-squares fits may produce negative values, which are undesirable for strict probability densities even if useful for classification (Duda, 2017). AI-Kolmogorov states directly that symbolic expressions are not guaranteed to be valid densities unless explicitly constrained, and its experiments on dijet distributions show sensitivity to negative predictions and sharp boundaries (Rajendram et al., 30 Mar 2026). MESSY mitigates this with a maximum-entropy form and cross-entropy correction, while discrete SDE adds explicit approximate-normalization and bounded-mass checks after search (Tohme et al., 2023, Liu et al., 20 May 2026).

High-dimensional scaling remains a central obstacle. Tensor-product polynomial or Fourier bases in orthogonal expansions grow combinatorially as $L^2$2 or $L^2$3 (Duda, 2017). AI-Kolmogorov addresses this through additive and multiplicative decomposition, but still treats high-dimensional SymDE as hard (Rajendram et al., 30 Mar 2026). Homothetic log-concave estimation avoids the usual curse only under strong structural assumptions (Xu et al., 2019). In discrete SymDE, grammar constraints and operator-cost profiles reduce a search space that otherwise explodes combinatorially (Liu et al., 20 May 2026). This suggests that future progress will continue to depend on decomposition, structural priors, symmetry, and validity-aware operator design rather than unconstrained formula search.

The field therefore spans several non-equivalent notions of “symbolic.” In some systems the symbol is an explicit density formula; in others it is a transport map, a maximum-entropy exponent, an orbit-conditioned estimator, or a density-guided symbolic discretization. What unifies them is the attempt to move from opaque estimation toward analyzable mathematical structure. The strongest recent formulations make that goal explicit: recover densities or PMFs in closed form, preserve probabilistic validity, and use symbolic structure not merely for post hoc explanation but as part of the estimation mechanism itself (Rajendram et al., 30 Mar 2026, Liu et al., 20 May 2026).