PySR: Symbolic Regression Toolkit

• PySR is a genetic-programming-based symbolic regression toolkit that derives closed-form expressions through evolutionary search.
• It utilizes parallel, multi-population strategies with mutation and crossover over user-defined operators to balance error and complexity.
• Validated across physics, astronomy, and dynamical systems, PySR delivers state-of-the-art accuracy with customizable constraints and efficient deployments.

PySR is a high‐performance genetic‐programming‐based symbolic regression library, written in Julia with a Python interface, designed to discover interpretable, closed-form mathematical expressions that accurately map input data to target outputs. PySR implements a parallel, multi-population evolutionary strategy, applying mutation and crossover operations over expression trees built from user-defined operators, variables, and constants. Its objective function incorporates both data-fit loss (typically MSE or a domain-specific variant) and an explicit complexity penalty, enabling Pareto optimization over accuracy and compactness. Beyond generic regression, PySR offers extensive customizability in operator selection, complexity regularization, constraint enforcement, and post-processing; it has rapidly become a tool of choice in physics, astronomy, dynamical systems, space weather, high-energy analysis, and beyond.

1. Algorithmic Principles and Evolutionary Framework

PySR frames symbolic regression as an evolutionary search over the space of mathematical expression trees. Each candidate model is a tree composed of operators (e.g. $+$, $-$, $\times$, $/$, $\sin$, $\exp$), input features, and trainable constants. Several populations (“islands”) are evolved in parallel, where each is subjected to rounds of selection, mutation (subtree replacement, constant perturbation), and crossover (subtree exchange between parent trees) (Brum et al., 27 Aug 2025). Fitness evaluation comprises both a numerical error metric (e.g., MSE, MAE, cross-entropy) and a parsimony or complexity penalty proportional to the node count of the expression tree. This penalty is tunable (parameter $\lambda$), and enables PySR to construct a Pareto frontier of candidates trading off error against symbolic complexity.

The model selection step extracts optimal candidates according to multi-objective criteria, including minimal loss, steepest improvement per complexity increment (score), or “hall-of-fame” bests within a loss threshold (Brum et al., 27 Aug 2025).

2. Hyperparameterization and Operator Customization

Users configure PySR via key hyperparameters, including: population size, number of evolutionary generations, choice of binary and unary operators, maximum allowed complexity, and complexity penalty weight. The operator set can be tailored precisely to the problem domain, ranging from purely polynomial or rational forms, to transcendental or piecewise functional blocks (e.g., $\sin$, $\exp$, $\log$, $\max$, $\min$) (Tsoi et al., 2023, Markidis et al., 25 Apr 2025). Nested constraints (e.g., forbidding $\sin(\sin(x))$) and domain-specific custom functions (e.g., physically-motivated kernels, threshold or saturation operators) can also be imposed.

Tuning guidance from benchmarking indicates that smaller populations and limited operator sets yield stable recoveries of known dynamical systems (Lorenz, Lotka–Volterra, SIR/SEIR), while more expressive libraries and weaker parsimony allow high-fidelity surrogates in multi-dimensional regression (Brum et al., 27 Aug 2025, Bendavid et al., 1 Aug 2025).

3. Loss Functions, Regularization, and Constraint Enforcement

PySR supports arbitrary loss functions defined on the fit residuals, including MSE, MAE, cross-entropy, and robust domain-specific variants. Custom constraints—e.g., soft penalties implementing physical laws (zero loading, monotonicity, finite slope), symmetry enforcement, or hypothesis testing (factorization, Regge behavior)—may be integrated directly into the fitness function (Fox et al., 2023, Dotson et al., 17 Apr 2025). In problems where background knowledge is essential, soft constraints multiply the loss by fixed penalty factors for each violated condition, whereas hard constraints prune the search space at every generation. Constraint checking leverages symbolic algebra (e.g., via SymPy), with a resulting $\sim10\times$ increase in per-generation runtime (Fox et al., 2023).

Bayesian symbolic regression frameworks (“Machine Scientist”) related to PySR can also incorporate constraint priors for sharper adherence to theoretical expectations, albeit at additional computational cost (slow MCMC mixing) (Fox et al., 2023).

4. Benchmark Performance in Scientific and Engineering Applications

PySR has demonstrated state-of-the-art accuracy and interpretability in diverse scientific domains:

• Dynamical systems: Exact recovery of analytical forms for nine benchmark models, including chaotic and epidemiological processes, with $R^2 > 0.99$ and low structural error (Brum et al., 27 Aug 2025).
• PDE discovery: In conjunction with denoising attention neural networks (ANN-PySR), robust extraction of governing equations under up to 200% noise and sparse sampling, consistently beating classical SINDy and genetic algorithms in coefficient accuracy and symbolic fidelity (Zhang et al., 22 Jun 2025).
• Astrophysics: Learning analytic neutron-star radius estimators $R(M,\Lambda)$ from GW-only data with $\lesssim 0.1$ km error across broad EOS variation, outperforming conventional parametrized-EOS approaches in speed and transparency (Bejger, 28 Apr 2025).
• LHC physics: Surrogate models for angular observables and CP-sensitive detector-level quantities, yielding compact, interpretable formulas and ultrafast evaluation in high-statistics analyses (Bendavid et al., 1 Aug 2025, Bahl et al., 8 Jul 2025).
• FPGA/Edge deployments: Symbolic regression surrogates, translated via hls4ml, slash inference latency and resource use, achieving $>90\%$ of neural network accuracy at $13\times$ lower latency and order-of-magnitude reductions in DSP/LUT footprints (Tsoi et al., 2023).

Representative accuracy metrics

Domain	Structural recovery	$R^2$	Typical Error
Lorenz/SIR/SEIR	Correct form	$>0.99$	$10^{-4}\,-\,10^{-3}$ trajectory MAE
PDEs (ANN-PySR)	Correct coefficients	N/A	$<0.05$ with 200% noise
GW NS radii	$R$ vs. ground-truth	N/A	$<0.1$ km mean, $<0.5$ km worst-case
LHC A$_i$	Surrogate/MC agree	N/A	$1$–$6\times 10^{-3}$ fit loss

5. Interpretability and Model Selection

A defining advantage of PySR is the explicit algebraic form of its output models. Each expression can be analyzed for limiting behavior, symmetries, and domain-specific structure, facilitating physical insight, rapid forward evaluation, and direct incorporation in inference pipelines (Bendavid et al., 1 Aug 2025, Dotson et al., 17 Apr 2025). Pareto-front model selection gives practitioners a choice of candidate formulas across a spectrum of complexity and accuracy, enabling “dialing up” of resource budget (e.g., in FPGA deployment) or prototyping for further scientific investigation (Tsoi et al., 2023).

In multi-replica studies, statistical clustering of expansion coefficients (ECC) can be used to quantify convergence and identify distinct families of local optima in discovered forms (Dotson et al., 17 Apr 2025).

6. Practical Guidelines for Application and Tuning

Recommendations from published studies include:

• Begin with minimal operator sets and moderate population sizes; introduce complexity only as needed (Brum et al., 27 Aug 2025).
• Employ strong complexity penalties during initial search to favor parsimonious forms; relax penalties for coefficient refinement.
• Leverage domain knowledge via custom operators, constraints, or input scaling (e.g., $\ln\Lambda$ for NS radii).
• For robust modeling under high noise, couple PySR with advanced pre-processing (filtering, denoising neural networks), batchwise training, and feature engineering (Zhang et al., 22 Jun 2025).
• Always inspect the Pareto frontier for physically interpretable models, validating candidates via forward simulation or domain-specific statistical tests.
• For FPGA deployments, select models targeting the latency–accuracy “knee-point” given resource constraints, and employ symbolic-to-HLS translation tools for synthesis (Tsoi et al., 2023).

7. Extensions, Limitations, and Ongoing Directions

Possible extensions include:

• Augmentation of the operator library (e.g., composite functions, delay kernels, domain-specific nonlinearities).
• Incorporation of physical invariances or asymptotics as hard/soft constraints (Bejger, 28 Apr 2025).
• Multiobjective optimizers to decouple loss and complexity (Markidis et al., 25 Apr 2025).

Limitations observed:

• Computational cost scales with population size, operator set, and constraint overhead.
• Hyperparameter tuning can be nontrivial, especially in highly coupled or high-dimensional settings (Brum et al., 27 Aug 2025).
• For extreme or pathological ground-truth models, constraint choices or operator limitations may exclude recovery of the true form (Fox et al., 2023).

Continued development addresses rapid parallelization, GPU integration, constraint-aware search, and hybrid symbolic–deep modeling.

---

PySR is now established as a versatile and robust symbolic regression toolkit for physical, dynamical, and inference-driven modeling, combining the interpretability of closed-form mathematics with the flexibility of evolutionary computation, and underpinning advances in scientific discovery across theory, data analysis, and hardware acceleration (Brum et al., 27 Aug 2025, Markidis et al., 25 Apr 2025, Tsoi et al., 2023, Dotson et al., 17 Apr 2025).