Bayesian Symbolic Regression

• Bayesian symbolic regression is a probabilistic approach to equation discovery that uses structural priors, likelihoods, and posterior distributions to infer mathematical expressions from data.
• It employs advanced inference algorithms such as MCMC, RJMCMC, and variational methods to navigate the NP-hard, multimodal, and transdimensional model spaces.
• The method leverages structured priors and domain-specific constraints to encapsulate scientific knowledge, enabling robust model averaging and uncertainty quantification.

Bayesian symbolic regression is the formulation of symbolic regression as probabilistic inference over mathematical expressions. Instead of ranking candidate equations solely by heuristic fit-complexity criteria, it places priors on expression structure and parameters, combines them with a likelihood for noisy observations, and targets a posterior distribution over equations, coefficients, and predictions. In this formulation, the central objects are posterior model probabilities, marginal likelihoods, description lengths, and posterior predictive distributions rather than a single heuristic optimum (Guimera et al., 22 Jul 2025, Jin et al., 2019).

1. Probabilistic formulation

At its core, Bayesian symbolic regression treats a candidate expression $m_i$ or $f_i$ as a statistical model. A standard formulation assumes noisy observations

$y^k = m^*(x^k, \theta^*) + \epsilon^k,$

with Gaussian observational noise of zero mean and unknown variance $\sigma^2$. For a candidate model $m_i$ with parameters $\theta_i$,

$$p(m_i,\theta_i\mid D)=\frac{p(D\mid m_i,\theta_i)\,p(\theta_i\mid m_i)\,p(m_i)}{p(D)},$$

and the posterior over expressions is obtained by marginalizing parameters,

$$p(m_i\mid D)=\frac{p(m_i)}{p(D)}\int d\theta_i\, p(D\mid m_i,\theta_i)\,p(\theta_i\mid m_i).$$

The maximum-a-posteriori equation is

$\hat m = \arg\max_{m_i} p(m_i\mid D),$

which the physicist-oriented treatment identifies with the model minimizing description length (Guimera et al., 22 Jul 2025).

An equivalent description-length formulation appears in Bayesian Machine Scientist-style work: $$p(f_i\mid D) = \frac{1}{Z}\int_{\Theta_i} d\theta_i\; p(D\mid f_i,\theta_i)\,p(\theta_i\mid f_i)\,p(f_i) = \frac{\exp[-\mathcal{L}(f_i)]}{Z}.$$ Here the posterior over symbolic models is written as an exponential family over description lengths, making explicit that Bayesian symbolic regression is simultaneously a model-selection problem and a compression problem (Fox et al., 2023).

The formulation need not be restricted to a single tree. One influential Bayesian symbolic regression model assumes an additive structure,

$f_i$0

so that complex signals are represented as linear combinations of several concise symbolic trees rather than one large expression (Jin et al., 2019). This suggests a broad conception of Bayesian symbolic regression: the latent symbolic object may be a single formula, an additive ensemble of formulas, or a support set in a screened symbolic library, provided that inference is posterior-based rather than purely heuristic.

A defining implication of the probabilistic view is that the full predictive target is not only a MAP equation but a posterior predictive distribution. One explicit expression is

$f_i$1

which makes model averaging intrinsic rather than optional (Guimera et al., 22 Jul 2025).

2. Prior structure, representation, and scientific knowledge

Bayesian symbolic regression depends critically on priors over symbolic form. In tree-based formulations, a symbolic expression is represented as a tree whose internal nodes are operators and leaves are variables or constants. Complexity is then controlled through a depth-dependent branching prior. Early Bayesian symbolic regression used

$f_i$2

for the probability that a node at depth $f_i$3 is non-terminal, together with priors over operators and terminal features (Jin et al., 2019). Hierarchical Bayesian symbolic regression later generalized this idea through

$f_i$4

so that deeper nodes are less likely to split, thereby regularizing tree depth and placing exponentially small mass on dense trees (Roy et al., 24 Sep 2025).

Beyond generic parsimony, priors can encode structural regularities of scientific expressions. One approach uses an $f_i$5-gram LLM over expression trees. Instead of counting operators independently, it models the probability of local operator configurations along tree paths, so that structurally different expressions with identical operator counts need not receive the same prior mass. The same work combines this with a Fractional Bayes Factor treatment of numerical parameter priors to make evidence-based model comparison well-defined even under weak or improper parameter priors (Bartlett et al., 2023).

A more formal route represents prior knowledge with probabilistic Regular Tree Expressions and evaluates the resulting tree prior via probabilistic tree automata. In that framework, symbolic expressions are treated as trees over a ranked alphabet, syntactic correctness is guaranteed by construction, and independent scientific constraints can be combined through Boolean operations because regular tree languages are closed under union, intersection, difference, and complement (Schneider et al., 2023). This is especially relevant when prior knowledge is modular—for example, translational invariance, dimensional restrictions, or known motifs.

Domain priors can also be explicitly scientific rather than purely structural. In adsorption thermodynamics, a corpus-based prior in Bayesian Machine Scientist was replaced by a thermodynamic prior that penalized violations of constraints such as

$f_i$6

These constraints were incorporated as soft penalties in the Bayesian prior rather than as hard feasibility rules, because hard constraints “severely hindered the search process,” whereas soft constraints preserved connectivity of the search space (Fox et al., 2023). This use of priors illustrates a recurrent theme in Bayesian symbolic regression: background knowledge changes the posterior distribution itself, rather than merely modifying a downstream fitness score.

3. Inference algorithms and computational hardness

The posterior over symbolic expressions is transdimensional, combinatorial, and typically multimodal, so Bayesian symbolic regression has been developed through a range of approximate inference algorithms. Classical Bayesian symbolic regression uses Metropolis-Hastings and reversible-jump MCMC over symbolic trees, with proposal moves such as Grow, Prune, Delete, Insert, ReassignOperator, and ReassignFeature. Reversible-jump updates are required because the number of parameters changes when tree structure changes, for example when linear-transform nodes are inserted or removed (Jin et al., 2019).

Bayesian Machine Scientist-style systems also perform MCMC over expression trees, using detailed-balance-preserving moves such as node replacement, root addition/removal, and elementary tree replacement (Fox et al., 2023). In missing-physics workflows, RJMCMC has been used as the second-stage interpretability mechanism after training a neural network inside a universal differential equation, yielding posterior samples over symbolic trees instead of a single post-hoc equation (Strouwen, 16 Mar 2026). Sequential Monte Carlo has also been proposed for Bayesian symbolic regression, using adaptive tempering, probabilistic selection, stratified resampling, rejuvenation, and a normalized marginal likelihood score to approximate the posterior over expressions under noise (Bomarito et al., 11 Dec 2025).

Other Bayesian formulations exploit partial conjugacy or continuous relaxations. In polynomial neural networks and polynomial neural ODEs, the inference target is a posterior over weights and biases, which is then propagated to a posterior over recovered polynomial coefficients; the compared methods were Laplace approximation, Hamiltonian Monte Carlo and NUTS, and variational inference, with Laplace found to be the best method for that class of problems (Fronk et al., 2023). VaSST replaces discrete operator and feature assignments by Binary Concrete and Gumbel-Softmax relaxations in “soft symbolic trees,” then optimizes a variational ELBO with black-box variational inference (Roy et al., 27 Feb 2026). Bayesian SISSO keeps the expensive sure independence screening stage deterministic and reformulates only the sparsifying operator stage as Bayesian inference over the screened support space, enabling exact posterior enumeration after screening (Bhattacharjee, 2 Jun 2026).

The computational difficulty of all such methods is illuminated by the formal hardness result for symbolic regression. Symbolic regression is posed as

$f_i$7

with decision version

$f_i$8

The optimization problem is NP-hard, and the decision problem remains NP-complete even under the restricted primitive set $f_i$9, $y^k = m^*(x^k, \theta^*) + \epsilon^k,$0, and a single observation (Virgolin et al., 2022). This suggests that exact posterior normalization or exact MAP search over realistic expression spaces is generally intractable, so Bayesian symbolic regression is usually implemented through approximation, restriction of the search space, or both.

4. Uncertainty quantification, ensembles, and theoretical guarantees

A major distinction between Bayesian symbolic regression and classical symbolic regression is that uncertainty in model structure is treated as a primary output rather than a nuisance. The physicist-oriented formulation emphasizes that the posterior over equations is the “true object of interest,” and that the full predictive distribution requires averaging over sampled models rather than relying on a single winner. This is especially important when $y^k = m^*(x^k, \theta^*) + \epsilon^k,$1 is far below 1 and several distinct equations have similar posterior weight, forming a Rashomon set (Guimera et al., 22 Jul 2025).

This posterior view changes both prediction and interpretation. Bayesian-model-averaged prediction appears explicitly in Bayesian SISSO: $y^k = m^*(x^k, \theta^*) + \epsilon^k,$2 where $y^k = m^*(x^k, \theta^*) + \epsilon^k,$3 indexes screened symbolic supports. The same framework computes posterior inclusion probabilities

$y^k = m^*(x^k, \theta^*) + \epsilon^k,$4

which quantify descriptor stability and reveal when a deterministic symbolic descriptor is merely one representative of a broader family of nearly degenerate alternatives (Bhattacharjee, 2 Jun 2026).

Structural uncertainty is also explicit in RJMCMC-based missing-physics recovery. There, posterior samples showed robust support for predator-prey interaction terms in Lotka–Volterra dynamics, but also many syntactically different trees encoding the same mathematical idea, including redundant variants and equivalent forms. The posterior therefore distinguishes robust interaction structure from syntactic non-uniqueness (Strouwen, 16 Mar 2026).

The Bayesian literature on symbolic regression also includes formal performance guarantees. One account stresses consistency: if the model class contains the true law and the prior does not exclude it, Bayesian symbolic regression selects the true model with probability tending to one as $y^k = m^*(x^k, \theta^*) + \epsilon^k,$5 (Guimera et al., 22 Jul 2025). A stronger structural-learning result is the near-minimax posterior concentration theorem for HierBOSSS, where under a well-conditioned expression design and a sub-critical split rule $y^k = m^*(x^k, \theta^*) + \epsilon^k,$6, the posterior contracts around the true regression function at a rate

$y^k = m^*(x^k, \theta^*) + \epsilon^k,$7

yielding near-parametric behavior in well-specified cases and near-minimax behavior in misspecified Hölder cases (Roy et al., 24 Sep 2025).

5. Scientific applications and domain-specific extensions

Bayesian symbolic regression has been developed for a wide range of scientific tasks in which interpretability, data scarcity, and uncertainty quantification are all material. In adsorption studies, Bayesian Machine Scientist was adapted with a thermodynamic prior and applied to rediscovering Langmuir-like, dual-site Langmuir, and BET isotherms from historical experimental datasets. Adding the thermodynamic prior increased the share of candidate expressions satisfying physically meaningful constraints, although the paper also stressed that incorrect constraints may hinder the search and that runtime increased by about an order of magnitude because constraints were checked symbolically with SymPy (Fox et al., 2023).

In dynamical systems, Bayesian polynomial neural networks and polynomial neural ODEs convert symbolic neural models from point estimation to posterior inference under Gaussian likelihood and Gaussian parameter priors. A Gaussian process regression preprocessing step smooths noisy trajectories for neural ODE training, and posterior samples over network parameters are propagated to empirical posteriors over recovered polynomial coefficients. The reported conclusion for that model class was that Laplace approximation was the best balance of accuracy, interpretability, and computational efficiency, especially for noisy data (Fronk et al., 2023).

Bayesian symbolic regression has also been coupled to experiment design. In Bayesian optimal experimental design for symbolic discovery, the unknown symbolic law is a discrete random model index $y^k = m^*(x^k, \theta^*) + \epsilon^k,$8 with model-specific continuous parameters $y^k = m^*(x^k, \theta^*) + \epsilon^k,$9, and the next design point is chosen by maximizing mutual information $\sigma^2$0, equivalently minimizing predictive entropy $\sigma^2$1. Hamiltonian Monte Carlo is used to sample from the model-specific parameter priors or posteriors that define the predictive distribution (Clarkson et al., 2022). In the fed-batch bioreactor case of Bayesian inference for missing physics, a well-designed experiment produced a posterior concentrated around Monod-like growth laws, whereas a random design left the posterior diffuse and structurally uncertain (Strouwen, 16 Mar 2026).

For differential-equation discovery from raw time series, integral Bayesian symbolic regression replaces derivative fitting by comparison of numerically integrated trajectories with observations. The posterior is written over the integrated governing function $\sigma^2$2,

$\sigma^2$3

so candidate differential equations are judged by their integrated solutions rather than by estimated derivatives. The reported results indicate that this integral approach had the best exact-recovery performance and robustness among the compared methods on logistic growth and Lotka–Volterra benchmarks, and it was then used to discover bacterial-growth equations from experimental time series (Cabanas-Tirapu et al., 18 Nov 2025).

Materials and chemistry applications show a complementary use of posterior structure. Bayesian SISSO on an $\sigma^2$4 Heusler-alloy magnetic-moment dataset yielded posterior probabilities over descriptor supports, feature-inclusion probabilities, Bayesian-model-averaged predictions, and near-nominal empirical coverage of 95\% predictive intervals, while exposing competing, physically related descriptor families that deterministic SISSO would render artificially unique (Bhattacharjee, 2 Jun 2026). HierBOSSS extended this structural-learning emphasis to simulated problems, canonical Feynman equations, and a single-atom catalysis dataset, combining full posterior inference with a marginal posterior-based Occam’s window selection rule and a graph-edit structural fidelity metric (Roy et al., 24 Sep 2025).

6. Related approaches, misconceptions, and unresolved issues

Not every symbolic-regression method that uses a Bayesian criterion is a full Bayesian symbolic regression method. Some methods are explicitly Bayesian-inspired rather than posterior-based. CADSR, for example, uses a Gaussian likelihood and a BIC reward

$\sigma^2$5

with $\sigma^2$6 defined as the number of nodes in the expression tree plus the number of constant tokens. This yields a Bayesian model-selection surrogate inside reinforcement learning, but it does not maintain a posterior distribution over equations and does not perform posterior inference over expression trees (Bastiani et al., 2024). Similarly, domain-aware Tree-RNN symbolic regression learns symbol priors from domain corpora and regularizes a policy by KL divergence toward those priors, but it is described as Bayesian in spirit rather than a classical Bayesian posterior model over expressions (Huang et al., 12 Mar 2025).

A common misconception is that Bayesian symbolic regression must output a single best formula. Several of the cited formulations argue the opposite. The posterior can remain diffuse because of noise, limited data, or genuine non-identifiability, and the scientifically proper output may then be a sampled ensemble, posterior support probabilities, or a posterior predictive distribution rather than a single canonical equation (Guimera et al., 22 Jul 2025). Another misconception is that priors merely encode simplicity. In the literature they also encode operator frequencies, local symbolic grammar, scientific constraints, domain corpora, screening relevance, and weak physical hypotheses such as the Slater–Pauling-inspired prior in Bayesian SISSO (Bartlett et al., 2023, Bhattacharjee, 2 Jun 2026).

The unresolved issues are largely computational and epistemic rather than conceptual. Exact inference remains difficult because symbolic regression is NP-hard (Virgolin et al., 2022). MCMC-based schemes may suffer from slow exploration of highly multimodal transdimensional spaces, motivating SMC, variational inference, and hybrid deterministic-Bayesian decompositions (Bomarito et al., 11 Dec 2025, Roy et al., 27 Feb 2026). Even in the probabilistic framework, there are regimes in which the true equation may be unlearnable from the available data: one treatment explicitly notes a learnability transition above some noise level, where no algorithm can recover the ground truth from the data alone (Guimera et al., 22 Jul 2025). A plausible implication is that the mature form of Bayesian symbolic regression is not a universal exact solver for equation discovery, but a family of principled approximate inference procedures that combine structural priors, posterior uncertainty, and tractable search over a provably hard expression space.