Symbolic Regression Automates Particle Physics Background Modeling With Single-Run Uncertainty Estimates

A research team has introduced SymbolFit, an open-source computational framework that automates parametric modeling by using symbolic regression to search for mathematical functions that fit binned experimental data while simultaneously providing uncertainty estimates within a single execution run. Published in Computing and Software for Big Science, the method addresses a longstanding operational bottleneck in high-energy physics data analysis at the CERN Large Hadron Collider, where researchers routinely construct parametric models to describe complex distributions for downstream statistical inference. Historically, fitting parametric curves to binned experimental observations has demanded an exhaustive and manual iterative process, requiring physicists to guess an adequate functional form before any regression fit could take place. The newly developed framework eliminates the necessity of predefining rigid closed-form equations by treating the underlying functional structure itself as a trainable parameter, allowing evolutionary machine learning algorithms to explore vast mathematical function spaces and discover optimal descriptive formulations without relying on manual trial and error.

In experimental collider studies searching for new physical phenomena, particle interactions are recorded across millions of collision events and compiled into binned histograms representing kinematic observables, most notably the invariant mass distributions of reconstructed final-state particles. Potential signatures of novel physics frequently manifest as subtle, narrow resonance peaks situated atop smoothly declining background event rates. Determining whether an observed fluctuation signifies a genuine discovery or a random statistical variation requires rigorous hypothesis testing that compares the recorded data against highly accurate parametric descriptions of both signal and background components. However, because collision events undergo intricate, multi-stage selection cuts, triggering criteria, and reconstruction filters, the resulting distributions incorporate arbitrary geometric and kinematic shape distortions. Consequently, the observed distributions rarely adhere to underlying closed-form theoretical laws derived strictly from first principles, leaving empirical parameterizations as the only viable mechanism for continuous modeling in standard particle search workflows.

Traditional parametric approaches force analysts to propose empirical function families manually, assessing variations through iterative trial fits and separate model selection criteria such as statistical F-tests. When experimenters analyzed high-mass trijet and dijet events at the Compact Muon Solenoid experiment, they crafted specialized empirical formulations incorporating fractional powers, logarithmic terms, and tunable hyperparameters to capture the falling background spectra. Although these bespoke equations satisfactorily modeled specific historical datasets, their inherent rigidity renders them susceptible to breakdown whenever detector configurations, operating energies, or event selection strategies shift. Alternative fitting techniques like Gaussian process regression offer flexible non-parametric distributions, but they require the numerical inversion of an n-by-n covariance matrix that scales with a cubic time complexity of order n cubed. This scaling creates severe computational bottlenecks for multi-variable analyses, while the resulting probabilistic predictions require complex, non-standard mathematical treatments to interface with standard high-energy physics statistical frameworks.

To overcome these limitations, SymbolFit establishes an automated, multi-tiered computational pipeline linking symbolic regression directly to rigorous statistical re-optimization. The architecture ingests binned data records containing bin centers, observed counts, and asymmetric one-standard-deviation error bounds that capture measurement and modeling uncertainties. Central to the framework is PySR, a high-performance Python package executing genetic programming to machine-search for functional forms without requiring a predefined template. Candidate functions are represented as mathematical expression trees whose internal and terminal nodes correspond to variables, constants, and algebraic operators such as addition, multiplication, division, exponentiation, hyperbolic tangents, and custom Gaussian functions. Through successive generations, candidate expression trees evolve via node mutations and subtree crossovers across multiple independent sub-populations. Guided by a customized chi-squared objective function that directly incorporates data uncertainties, this evolutionary search discovers diverse functional structures that closely conform to arbitrary data shapes while constraining overall tree complexity.

Because genetic programming focuses strictly on finding algebraic topologies and does not inherently compute parameter uncertainties, SymbolFit implements a crucial parameterization and re-optimization fit stage using the LMFIT minimization library. Once symbolic regression yields a batch of promising functional equations, the framework freezes their symbolic expressions and automatically parameterizes all embedded numerical constants into variable parameters. LMFIT then conducts a secondary nonlinear least-squares fit against the data points to minimize the chi-squared loss and refine parameter magnitudes. If complex formulations trigger numerical convergence failures during minimization, an automated fallback loop progressively reduces the degrees of freedom by fixing selected parameters to their initial values until the algorithm converges with all relative parameter errors remaining beneath an established threshold. This secondary fitting protocol reliably yields best-fit parameter values, symmetric parameter uncertainties, and a complete correlation matrix, supplying the rigorous error modeling essential for downstream statistical evaluation.

Unlike many symbolic regression studies where the primary scientific objective is symbolic interpretability or the extraction of exact physical laws, SymbolFit focuses purely on constructing valid, highly adaptive parametric curves for statistical inference. Because the generated outputs are explicit closed-form mathematical equations, they integrate seamlessly with widely adopted CERN statistical tools including Combine and pyhf without requiring structural code alterations. This direct architectural compatibility bypasses the intricate likelihood modifications and non-standard mathematical treatments required when deploying Bayesian machine-learning predictors such as Gaussian processes. Furthermore, by generating an entire batch of viable candidate equations during a single fit configuration, the framework enables particle physicists to treat the functional model choice itself as a systematic uncertainty through established discrete profiling methods. Analysts can evaluate objective goodness-of-fit metrics, inspect automatically exported diagnostic visual plots, and quantify fit quality across differing mathematical functional complexities, providing complete control over the balance between expressive modeling fidelity and structural simplicity.

The team validated the framework on synthetic data lacking an underlying theoretical functional form, designated as Toy Dataset 1. This one-dimensional binned spectrum featured an engineered sharp peak coupled to an extended high-energy tail, paired with synthetic vertical uncertainty bars. Using a configuration constrained to a maximum complexity of sixty nodes with binary operators, unary functions, and a custom Gaussian operator, a single SymbolFit execution produced forty-six candidate functions spanning complexity scores from one to sixty. The re-optimization fit substantially improved the reduced chi-squared scores across candidate formulas compared to initial symbolic regression outputs. For instance, candidate function twenty-seven possessed six adjustable parameters whose best-fit values and correlation matrix defined a robust uncertainty model. Generating an ensemble of functions by sampling parameter values from a multidimensional normal distribution allowed researchers to map the sixty-eight percent quantile uncertainty band, verifying that the automated pipeline successfully tracks complex asymmetrical structures while properly estimating statistical uncertainties.

To establish practical efficacy in authentic collision physics, the authors benchmarked SymbolFit against real proton-proton collision datasets recorded at a center-of-mass energy of thirteen tera-electronvolts by the CMS experiment during Run 2. The primary demonstration evaluated background modeling in a dijet invariant mass spectrum obtained from public data repositories, alongside validations across diphoton, trijet, paired-dijet, and dimuon resonance channels where CMS originally reported no evidence of new physics. In their dijet analysis, CMS had traditionally relied on an empirical four-parameter mathematical equation to model the continuum background. To evaluate whether machine-searched equations could reliably isolate prospective physics discoveries without absorbing real signals into the background estimate, the researchers synthesized pseudodata by injecting an artificial, narrow Gaussian resonance into the empirical dijet spectrum at an invariant mass of 3.1 tera-electronvolts, with a peak width of 0.2 tera-electronvolts, an amplitude strength of ten, and ten percent random bin perturbation noise.

The investigation applied a standard blinded analysis protocol, masking the invariant mass range between 2.7 and 3.5 tera-electronvolts surrounding the injected signal peak while fitting the continuum background. Three distinct SymbolFit runs were initiated using the exact same operational configuration but varying random initialization seeds, producing three separate background models whose reduced chi-squared values converged close to unity after re-optimization. When unblinded background-only fits were executed on the complete pseudodata, all three symbolic regression models exhibited a distinct excess of events at the resonance location, confirming the background was not artificially overfitted. Subsequent signal-plus-background fits incorporated a parameterized Gaussian template, successfully absorbing the local event excess. In extensive signal extraction trials testing injected mass values between 2980 and 3150 gigaelectronvolts and signal strengths ranging from 2 to 38, all three symbolic regression models accurately recovered true injected signal parameters within expected statistical uncertainties, demonstrating performance comparable to or exceeding the handcrafted CMS empirical formulation.

Beyond one-dimensional background modeling, the researchers demonstrated the framework’s scalability to multidimensional data distributions, addressing experimental needs such as deriving smooth data-to-simulation scale factors and transfer factors for data-driven background estimation. In collider analyses, corrections for jet energy scales, heavy-flavor tagging, and tau lepton identification are routinely parameterized across multiple kinematic variables like transverse momentum and pseudorapidity, traditionally necessitating coarse binned weights or cumbersome products of Bernstein polynomials. Using Toy Dataset 3, composed of three distinct two-dimensional binned distributions manually generated without underlying mathematical functions, single SymbolFit executions generated batches of 22, 41, and 46 candidate functions that smoothly described the surfaces with well-behaved uncertainties. By replacing laborious manual trial-and-error procedures with an automated symbolic regression search, the framework establishes a versatile, computationally efficient methodology for parametric modeling across high-energy physics analyses and other data-intensive scientific disciplines where exact theoretical distributions remain fundamentally unknown.

Subject of Research: Automated parametric modeling and uncertainty estimation using symbolic regression in high-energy physics

Article Title: SymbolFit: Automatic Parametric Modeling with Symbolic Regression

Article References: Tsoi, H. F., Rankin, D., Caillol, C., Cranmer, M., Dasu, S., Duarte, J., Harris, P., Lipeles, E., & Loncar, V. (2025). SymbolFit: Automatic Parametric Modeling with Symbolic Regression. Computing and Software for Big Science, 9(1), Article 12. https://doi.org/10.1007/s41781-025-00140-9

Image Credits: AI Generated

DOI: 10.1007/s41781-025-00140-9

Keywords: Symbolic Regression, SymbolFit, High-Energy Physics, Large Hadron Collider, Parametric Modeling, Uncertainty Estimation, CMS Experiment, Machine Learning, Automatic, Parametric, Modeling, Symbolic