📊 statistics

Probabilistic Symbolic Regression for Equation Discovery via Operator-induced and Regularized Symbolic Forests

This paper introduces a probabilistic symbolic regression framework that utilizes operator-induced and regularized symbolic forests to simultaneously achieve high predictive accuracy, optimal expression complexity, and robust uncertainty quantification, while providing theoretical guarantees for posterior concentration under both exact and misspecified symbolic conditions.

Original authors: Somjit Roy, Pritam Dey, Bani K. Mallick, Debdeep Pati

Published 2026-07-28

📖 2 min read☕ Coffee break read

CC BY 4.0

Original authors: Somjit Roy, Pritam Dey, Bani K. Mallick, Debdeep Pati

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a detective trying to solve a mystery, but instead of looking for fingerprints or footprints, you are looking for the hidden mathematical rules that govern the universe. This is the world of scientific machine learning, a field where computers try to learn from data not just to make predictions (like guessing tomorrow's weather), but to discover the actual equations that explain how things work. Usually, computers are great at finding patterns, but they often give us answers that look like a tangled ball of yarn—super accurate but impossible for a human to read or understand. Symbolic regression is the special tool scientists use to untangle that yarn. It searches for simple, clean formulas (like $F=ma$) that fit the data, rather than complex, black-box computer models. The big challenge, however, is that there are so many possible formulas that finding the right one is like finding a single needle in a haystack the size of a galaxy, and it's easy to get tricked by "noise" (random errors in the data) into thinking a fake formula is the real deal.

Enter BayeSymX, a new method introduced by Somjit Roy and his team that acts like a super-smart, cautious detective for finding these mathematical laws. Instead of just guessing one answer and hoping it's right, BayeSymX treats the search as a probabilistic game. Imagine a forest where every tree represents a different possible equation. BayeSymX doesn't just pick the tallest tree; it grows a whole forest of candidate equations, weighs how likely each one is to be true based on the data, and uses a special "Occam's Razor" rule to prune away the overly complicated ones. The researchers found that this approach is incredibly effective: it not only predicts data better than other top methods but also recovers the correct, simple structure of famous physics equations (like those from Feynman's lectures) even when the data is noisy. In a real-world test involving new materials for clean energy catalysts, BayeSymX successfully identified compact, scientifically meaningful "descriptors" that other methods missed or buried under layers of unnecessary complexity. The team proved mathematically that their method is reliable, showing that as more data comes in, the method's confidence in the correct equation grows rapidly, effectively zeroing in on the truth even when the perfect formula isn't immediately obvious.

Technical Summary: Probabilistic Symbolic Regression for Equation Discovery via Operator-induced and Regularized Symbolic Forests

1. Problem Statement

Symbolic Regression (SR) aims to discover interpretable analytical expressions that govern input-output relationships directly from data, a central task in scientific machine learning. While existing SR methods (e.g., genetic programming, deep symbolic regression, and compressed sensing approaches) have shown promise, they face significant statistical and computational challenges:

Heuristic Reliance: Many methods rely on stochastic search heuristics that struggle to balance predictive accuracy with expression complexity, particularly in noisy, small-sample scientific datasets.
Uncertainty Characterization: Current approaches offer limited characterization of symbolic uncertainty, often returning a single "best" expression without quantifying the plausibility of alternative structural explanations.
Theoretical Gaps: There is a scarcity of theoretical treatments regarding posterior concentration rates for symbolic regression, particularly under conditions of misspecification or non-identifiability (where algebraically distinct expressions yield identical predictions).

The paper addresses these gaps by proposing a unified probabilistic framework that treats symbolic expressions as ensembles of trees, enabling full uncertainty propagation and rigorous theoretical guarantees.

2. Methodology: The BayeSymX Framework

The authors introduce BayeSymX (Bayesian Symbolic regression forests for eXpression discovery), a probabilistic framework that models the unknown regression surface $f$ as an affine combination of symbolic trees (a "symbolic forest").

2.1 Model Structure

The model assumes observations $y_i = f(x_i) + \epsilon_i$ , where:
$y_i = \beta_0 + \sum_{j=1}^K g(x_i; T_j)\beta_j + \epsilon_i$
Here, $g(x; T_j)$ represents the evaluation of the $j$ -th symbolic tree $T_j$ , and $\beta$ are outer regression coefficients. The trees are constructed recursively from a library of primary features and mathematical operators (unary and binary).

2.2 Prior Specifications

The framework employs a hierarchical Bayesian specification designed to regularize complexity and learn data-adaptive preferences:

Tree Topology Prior: A depth-dependent splitting probability $p_m = \alpha_0(1+m)^{-\delta_0}$ penalizes deep trees, enforcing a form of Occam's razor that favors simpler representations.
Operator and Feature Priors: Unlike fixed-weight approaches, BayeSymX uses Dirichlet priors on tree-specific operator and feature assignment weights. This allows the model to learn which operators and features are relevant for specific trees in a data-adaptive manner.
Regression Coefficients: Conjugate Normal-Inverse-Gamma (NIG) priors are placed on the outer coefficients $\beta$ and noise variance $\sigma^2$ , ensuring full uncertainty propagation across all model parameters.

2.3 Posterior Inference

Marginalization: The continuous parameters $(\beta, \sigma^2)$ are marginalized out analytically using NIG conjugacy, yielding a joint marginal posterior (JMP) over the discrete space of symbolic forests.
Sampling: A Metropolis-within-partially-collapsed Gibbs sampler is used to explore the symbolic expression space. The sampler employs seven local tree moves (grow, prune, subtree replacement, delete, insert, change feature, change operator) to navigate the discrete space.
Model Selection (Occam's Window): Instead of selecting a single best tree, BayeSymX utilizes an Occam's window approach. It retains a set of high-posterior-probability forests ( $J_r$ ) to capture uncertainty across multiple plausible symbolic models.
Refinement: A post-MCMC refinement step uses the Bayesian Information Criterion (BIC) to prune redundant trees and algebraically simplify the final expressions.

3. Key Contributions

3.1 Theoretical Guarantees

The paper establishes novel posterior concentration results for symbolic regression, a field previously lacking rigorous theoretical treatment:

Approximate Realizability: Under mild regularity assumptions, the authors prove that the posterior concentrates around the true data-generating function $f_0$ at a rate governed by the trade-off between empirical approximation error and a newly derived symbolic complexity scale ( $\mathfrak{C}_{K,S,n}$ ).
Near-Parametric Rates: In the case of exact finite symbolic representability, the framework achieves a near-parametric concentration rate of $O(n^{-1/2}(\log n \log \log n)^{1/2})$ .
Misspecification and Oracle Inequalities: Under symbolic misspecification (where $f_0$ is not in the model class), the paper establishes a sharp oracle concentration result. The posterior concentrates around the optimal population approximation error without requiring the existence of a finite set of Kullback-Leibler minimizers or specialized testing conditions typically needed in classical misspecification theory.
Non-Identifiability Handling: The guarantees are formulated at the level of predictive functions, acknowledging that multiple distinct symbolic structures may represent the same function.

3.2 Methodological Innovations

Operator-Induced Forests: The use of forests (ensembles) rather than single trees allows for additive scientific structures while maintaining interpretability.
Data-Adaptive Learning: The Dirichlet priors on operator/feature weights enable the model to adaptively learn structural preferences, avoiding the rigid constraints of fixed-weight priors found in previous Bayesian SR methods (e.g., BSR).
Uncertainty-Aware Summarization: The Occam's window strategy provides a principled way to report multiple competing scientific hypotheses rather than a single point estimate.

4. Empirical Results

The authors evaluate BayeSymX against state-of-the-art competitors (including gplearn, operon, PySR, DSR, QLattice, SISSO++, BMS, and BSR) on two distinct benchmarks:

4.1 Feynman Equations (SRBench)

Setup: Recovery of 5 physical laws from the Feynman Lectures under varying noise levels and structural complexities.
Findings: BayeSymX consistently achieved a superior balance between predictive accuracy (lowest test RMSE), symbolic parsimony (compact expressions), and exact structural recovery. Competing methods often either failed to recover the correct structure or produced overly complex expressions to achieve similar accuracy. BayeSymX demonstrated robustness to increasing noise levels where other methods degraded significantly.

4.2 Oxide Perovskite Catalyst Discovery

Setup: Discovery of "materials genes" (descriptors) linking catalyst composition to oxygen evolution reaction (OER) activity.
Findings: BayeSymX identified compact, scientifically interpretable descriptor expressions (26–40 nodes) that recovered known structure-activity relationships (e.g., involving tolerance factor $\mu$ , electronegativities $\chi_A, \chi_B$ ). In contrast, high-accuracy competitors like operon produced unwieldy expressions (90–104 nodes), while compact methods like PySR showed lower predictive performance. BayeSymX occupied the Pareto frontier of the accuracy-complexity trade-off.

5. Significance and Claims

The paper claims that BayeSymX represents a significant advancement in probabilistic symbolic regression by:

Unifying Structure and Uncertainty: Providing a framework that jointly learns symbolic structure, controls complexity via regularization, and quantifies uncertainty across multiple plausible models.
Theoretical Rigor: Offering the first posterior concentration guarantees for symbolic regression that handle both exact representability and misspecification, establishing near-parametric rates and sharp oracle inequalities.
Scientific Utility: Demonstrating that probabilistic approaches can outperform heuristic and deep learning-based methods in recovering interpretable scientific laws, particularly in noisy, small-sample regimes typical of materials discovery and physics.

The authors conclude that the framework is particularly well-suited for scientific discovery settings where domain knowledge guides feature selection, but the underlying functional form remains unknown and requires robust handling of structural uncertainty.