SymbolicArena: A Unified Infrastructure for Benchmark Distillation and Dynamic Evaluation in Symbolic Regression

Ziwen Zhang; Xiju Wu; Yuheng Jing; Runxiang Wang; Boxiao Wang; Yifan Zang; Yifan Zhang; Yang Wang

cs.LGSep 28, 2026

SymbolicArena: A Unified Infrastructure for Benchmark Distillation and Dynamic Evaluation in Symbolic Regression

Authors: Ziwen Zhang, Xiju Wu, Yuheng Jing, Runxiang Wang, Boxiao Wang, Yifan Zang, Yifan Zhang, Yang Wang, +4 more

Organizations: School of Artiﬁcial Intelligence, University of Chinese Academy of Sciences · C2DL, Institute of Automation, Chinese Academy of Sciences · Nanjing University of Science and Technology · School of Future Technology, University of Chinese Academy of Sciences

Abstract

Symbolic regression (SR) seeks concise and interpretable mathematical expressions from data for scientific equation discovery. Existing SR benchmarks face a tradeoff between evaluation cost and benchmark validity. Repeated evaluation of large task pools is expensive, and compact benchmarks lack systematic evidence of preserved task diversity and algorithm discriminability. SymbolicArena provides a unified infrastructure for benchmark distillation and dynamic evaluation. The framework standardizes 664 heterogeneous tasks with executable ground truth expressions and distills the Full Task Set into Core50, a validated benchmark of 50 tasks. The distillation process preserves task coverage and algorithm discrimination under explicit balance constraints. SymbolicArena applies a unified execution protocol to heterogeneous SR algorithms and produces comparable outputs and search trajectories. Multi Axis Evaluation characterizes numerical quality, symbolic quality, and search behavior. Core50 reduces evaluation workload by 92.5% and maintains agreement with Full Task Set evaluations. Experiments show that SymbolicArena achieves 72.6% to 86.7% lower approximation error than alternative selectors, further supporting its fidelity to the Full Task Set. Evaluation reveals a substantial gap between numerical fitting and symbolic recovery across current SR methods, suggesting that reliable equation recovery remains an open challenge.

Figures & tables

Benchmark	Tasks	Focus	HP	CP	VC	ST
SRBench ( La Cava et al., 2021 )	252	Algorithm comparison	$\checkmark$	$\checkmark$	$\times$	$\times$
SRBench 2025 ( Imai Aldeia et al., 2025 )	24	Accuracy and complexity	$\checkmark$	$\checkmark$	$\times$	$\times$
SRSD ( Matsubara et al., 2024 )	240	Scientific equation rediscovery	$\times$	$\times$	$\times$	$\times$
LLM-SRBench ( Shojaee et al., 2025b )	240	LLM based equation discovery	$\times$	$\times$	$\times$	$\times$
Nguyen ( Uy et al., 2011 ) Keijzer ( 2003 ) Korns ( 2011 ) Vladislavleva ( 2009 )	12 15 15 8	Equation recovery and extrapolation	$\times$	$\times$	$\times$	$\times$
SymbolicArena	664 $\rightarrow$ 50	Validated benchmark distillation	$\checkmark$	$\checkmark$	$\checkmark$	$\checkmark$

Table 1: Comparison with representative symbolic regression benchmarks.

Figure 1: Overview of SymbolicArena . a. Benchmark distillation retains 664 tasks from 800 candidates, forming the Full Task Set . The Calibration Set tasks support selection of the four probes . Task properties and probe responses across Full Task Set guide selection of 50 tasks in Core50 . b. SymbolicArena runs 15 SR algorithms on Core50 using a common protocol and records outputs and search trajectories. Multi Axis Evaluation reports numerical quality (ID and OOD), symbolic quality (SYM and MIN), and search behavior (EFF and STAB).

Panel of Figure 2 — Figure 2: Aggregate fidelity and task representation. a. PCA projection of structural and probe response descriptors for Full Task Set , Calibration Set , and Core50 . b. Mean task information and aggregate score MAE for Core50 and six alternative selectors. c. Source proportions for Full Task Set and the evaluated 50 task subsets.

Panel of Figure 3 — Figure 3: Numerical score and rank agreement. ID and OOD scores for nine algorithms on Full Task Set and Core50 . FePySR, SymbolFit, and JAXSR are excluded from benchmark construction. Scores use identical evaluation records and a one hour budget.

		Numerical		Symbolic		Search
Paradigm	Algorithm	ID $\uparrow$	OOD $\uparrow$	SYM $\uparrow$	MIN $\uparrow$	EFF $\uparrow$	STAB $\uparrow$
Structured Search	iMCTS	77.91	73.35	42.15	80.37	91.66	23.31
	QLattice	34.79	25.37	29.82	55.12	87.81	31.12
	JAXSR	37.99	30.21	32.00	75.03	99.99	97.43
Evolutionary Search	PySR	70.34	66.28	44.60	82.62	91.94	31.03
	PyOperon	38.89	30.64	25.06	44.88	97.28	15.34

Table 2: Six-axis profiles of 15 symbolic regression methods on Core50 (clean). Methods are grouped by their primary search paradigm. Higher values denote better performance on each axis.

Figure 4: Per algorithm six axis profiles on Core50 . Each panel displays the six axis scores of one algorithm. Normalized scores are multiplied by 100 for a common 0 to 100 display scale. Gray profiles show the remaining algorithms for reference. The panels reveal differences beyond a single numerical ranking.

Appendix figures & tables22 assets

Supplementary material from the paper’s appendix.

Appendix

Field group	Stored information
Identity and provenance	dataset_id , source family, and subgroup. Original basename and source benchmark
Variables and target	Ordered feature names and target name. Number of features and dummy variable indicator
Ground truth expression	Raw source formula and executable expression. Normalized symbolic forms and expression tree
Data splits	Fixed train, ID test, and OOD test samples. Split sizes and parameters required for deterministic regeneration
Structural metadata	Formula complexity and operator types. Variable count bin and OOD type
Duplicate metadata	Semantic duplicate group identifier and provenance links used for duplicate auditing

Appendix

Table 3: Unified schema of Full Task Set .

Dimension	Categories and task counts
Task sources	LLM-SRBench (240), SRSD (232), SRBench 1.0 (133), Keijzer (15), SRBench 2025 first-principles subset (12), Nguyen (12), Korns (12), Vladislavleva (8)
Formula types	Rational (282), Trigonometric (137), Polynomial (124), Exponential / logarithmic (86), Mixed (34), Unknown (1)
Formula complexity	Simple (284), Moderate (188), Complex (192)

Appendix

Table 4: Composition of Full Task Set (664 tasks).

Figure 5: Dual probe mining diagnostic. Each point represents one task in the PySR and LLM-SR response space. Off diagonal points indicate large performance differences. Selected Calibration Set tasks cover distinct ID and OOD response regions.

Figure 6: Calibration usability of candidate algorithms. Points show finite result rate $q(a)$ and explosion rate $b_{100}(a)$ on Calibration Set . The reference region marks high availability and low explosion rates. All 12 algorithms remain eligible for probe selection.

Rank	Candidate combination	Mean usability
1	DSO, PyOperon, iMCTS, uDSR	0.971250
2	DSO, PyOperon, QLattice, uDSR	0.969375
3	PyOperon, iMCTS, QLattice, uDSR	0.968750

Appendix

Table 5: Highest scoring feasible four probe combinations. Mean calibration usability uses

\alpha=0.5

and

\tau=100

. Each reported subset covers three primary paradigms.

Selector	MAE	Empirical 95% range
Uniform random	$0.5298\pm 0.2550$	$[0.1605,\,1.1819]$
Smoothed family stratified	$0.5438\pm 0.2072$	$[0.2257,\,1.0102]$
Metadata diverse	$0.5060$	—
Response K-medoids	$0.6044$	—
Top information	$1.0339$	—
Difficulty balanced	$1.0430$	—

Appendix

Table 6: Aggregate score fidelity of 50 task selectors. Random selectors report mean MAE

\pm

standard deviation and empirical 95% ranges across 100 draws from Full Task Set . Other selectors report point estimates. Core50 is fixed, and the final row gives a proportional family quota check.

Idx	Task	Variables	Ground Truth
SRSD 15 tasks
1	Feynman II.27.18	$x_{0},x_{1},x_{2}$	$8.854\times 10^{-12}x_{0}^{2}$
2	Feynman I.43.31	$x_{0},x_{1}$	$1.3806\times 10^{-23}x_{0}x_{1}$
3	Feynman II.34.11	$x_{0},\allowbreak\dots,\allowbreak x_{6}$	$\frac{x_{0}x_{2}x_{4}}{2x_{5}}$
4	Feynman I.18.16	$x_{0},\dots,x_{5}$	$x_{0}x_{1}x_{3}\sin x_{4}$
5	Feynman II.34.2a	$x_{0},\dots,x_{4}$	$\frac{x_{1}x_{3}}{2\pi x_{4}}$

Appendix

Table 7: Core50 ground truth formula manifest. Tasks are grouped by source family. Variables follow the recorded order, and constants use at most four decimal places. Exact expressions and task identifiers are provided in the released metadata.

Figure 7: SymbolicArena execution architecture. A shared task contract connects isolated algorithm wrappers to standardized evaluation artifacts.

Algorithm	Probe4	$H$	Inclusion (%)	Any optimum (%)
uDSR	Yes	0.9875	100.0	100.0
PyOperon	Yes	0.9600	98.8	99.3
DSO	Yes	0.9700	88.1	91.5
iMCTS	Yes	0.9675	81.4	85.4
QLattice		0.9600	30.5	35.0
gplearn		0.9300	1.2	1.9

Appendix

Table 8: Calibration usability and bootstrap inclusion of the 12 eligible probe candidates. Inclusion frequency averages tied optima uniformly.

Removed	Replacement	Best $U$	$\Delta U$	Paradigms
iMCTS	QLattice	0.969375	0.001875	3
DSO	QLattice	0.968750	0.002500	3
uDSR	QLattice	0.964375	0.006875	3
PyOperon	gplearn	0.963750	0.007500	3

Appendix

Table 9: Probe4 replacement analysis.

U

denotes panel usability, and

\Delta U

measures the loss from the default Probe4 score

U=0.971250

Perturbed weight	Weight range	Median Jaccard vs. default	Overlapping tasks	Median within- group Jaccard	Full-ranking rate
Coverage	$[0.40,0.50]$	0.724	42	0.667	$100\%$
MeanInfo	$[0.30,0.40]$	0.710	41.5	0.667	$100\%$
Balance	$[0.15,0.25]$	0.710	41.5	0.667	$100\%$

Appendix

Table 10: Single-parameter weight perturbation results. Jaccard similarity and overlapping tasks are measured against the default-weight selection under 30 random restarts per setting. Within-group Jaccard similarity compares the 30 restarts of one setting.

Optimal set	Selections	Share	Jaccard vs. default
S0 ( Core50 )	2012	$40.24\%$	1.000
S1	1803	$36.06\%$	0.724
S2	1173	$23.46\%$	0.786
S3	12	$0.24\%$	0.754

Appendix

Table 11: Distribution of optimal sets under joint weight sampling. Selections report the number of the 5,000 sampled weight vectors that produce each optimal set. S0 denotes the default-weight selection.

Statistic	Result
Mean Jaccard vs. default	0.850
Minimum overlapping tasks	42 / 50
Full-ranking rate	$100\%$
Tasks selected in all runs	37
Tasks in $\geq 90\%$ of runs	38
Tasks in $\geq 75\%$ of runs	42

Appendix

Table 12: Joint random weight analysis results. Statistics summarize the 5,000 sampled weight vectors.

Figure 8: Spearman correlations among the six evaluation axes across 15 algorithms on clean Core50 .

Figure 9: Mean scores and 95% task bootstrap intervals on clean Core50 . Points denote mean scores and bars denote percentile intervals from 1,000 resamples of 50 task blocks. Each sampled block retains the three recorded runs. Scores use the 0 to 100 display scale.

		ID				OOD
Paradigm	Algorithm	Clean $\uparrow$	1% $\uparrow$	5% $\uparrow$	$\boldsymbol{\Delta_{5\%}}$ $\uparrow$	Clean $\uparrow$	1% $\uparrow$	5% $\uparrow$	$\boldsymbol{\Delta_{5\%}}$ $\uparrow$
Structured Search	iMCTS	77.91	55.08	46.69	-31.22	73.35	49.40	40.73	$-32.62$
	QLattice	34.79	32.36	29.18	$-5.61$	25.37	23.09	20.27	$-5.10$
	JAXSR	37.99	32.62	31.00	$-6.99$	30.21	25.91	24.56	$-5.65$
Evolutionary Search	PySR	70.34	48.42	41.05	$-29.29$	66.28	41.50	35.53	$-30.75$
	PyOperon	38.89	36.29	32.97	$-5.92$	30.64	27.85	24.42	$-6.22$

Appendix

Table 13: Numerical robustness of 15 symbolic regression methods under training label noise. ID and OOD scores are reported under clean, 1%, and 5% noise conditions.

\Delta_{5\%}

denotes the score change from clean to 5% noise.

		Numerical		Symbolic		Search
Paradigm	Algorithm	ID $\uparrow$	OOD $\uparrow$	SYM $\uparrow$	MIN $\uparrow$	EFF $\uparrow$	STAB $\uparrow$
Structured Search	iMCTS	55.08	49.40	37.22	72.88	89.34	12.82
	QLattice	32.36	23.09	29.59	54.60	90.87	38.32
	JAXSR	32.62	25.91	30.88	76.16	99.95	91.33
Evolutionary Search	PySR	48.42	41.50	30.53	60.26	86.08	6.02
	PyOperon	36.29	27.85	25.14	43.76	97.30	15.92

Appendix

Table 14: Six-axis profiles of 15 symbolic regression methods on Core50 at

1\%

training-label noise (

\sigma=0.01

). Methods are grouped by their primary search paradigm. Higher values denote better performance on each axis.

		Numerical		Symbolic		Search
Paradigm	Algorithm	ID $\uparrow$	OOD $\uparrow$	SYM $\uparrow$	MIN $\uparrow$	EFF $\uparrow$	STAB $\uparrow$
Structured Search	iMCTS	46.69	40.73	36.43	66.26	86.11	5.24
	QLattice	29.18	20.27	29.19	54.99	85.95	29.66
	JAXSR	31.00	24.56	31.29	77.47	99.95	91.95
Evolutionary Search	PySR	41.05	35.53	30.66	59.16	85.28	6.73
	PyOperon	32.97	24.42	23.60	39.83	97.27	15.33

Appendix

Table 15: Six-axis profiles of 15 symbolic regression methods on Core50 at

5\%

training-label noise (

\sigma=0.05

). Methods are grouped by their primary search paradigm. Higher values denote better performance on each axis.

Paradigm	Method	Description
Structured Search	iMCTS ( Huang et al., 2025 )	Monte Carlo tree search for guided expression exploration.
	QLattice ( Broløs et al., 2021 )	Probabilistic graph search for compact symbolic models.
	JAXSR ( Kitchin, 2024 )	Structured symbolic model construction with optimization in JAX.
Neural Policy Search	uDSR ( Landajuela et al., 2022 )	Neural symbolic search with evolutionary refinement.
	DSO ( Petersen et al., 2021 )	Policy search using risk seeking optimization.
Evolutionary Search	PySR ( Cranmer, 2023 )	Population search with complexity aware model selection.

Appendix

Table 16: Symbolic regression algorithms evaluated in SymbolicArena . Bold marks validation algorithms.

Algorithm	Search setting	Operators and basis	Fixed settings
gplearn	Population 1000.	add, sub, mul, div, sqrt, log, sin, cos	1 job, tournament 20, initial depth 2 to 6, MAE metric, parsimony 0.001, crossover 0.9, subtree, hoist, and point mutation 0.01.
PyOperon	Population 500, pool 500.	add, sub, mul, div, aq, exp, log, sin, cos, tanh, sqrt, square, constants, variables	4 threads, max length 50, max depth 10, tournament 5, keep best reinsertion, LM optimizer, local search probability 1.0.
PySR	8 populations of size 64, 500 cycles per iteration.	Binary $+,-,\times,/$ . Unary square, cube, exp, log, sin, cos.	Serial execution with 1 process, max size 30, max depth 10, parsimony 0.001, deterministic mode, precision 32, best model selection.
QLattice	BIC criterion.	QLattice/Feyn internal symbolic basis	4 threads, regression mode, BIC model selection, 4 significant digits.
DSO	Batch size 1000.	add, sub, mul, div, sin, cos, exp, log	4 batch cores, reward inv NRMSE, $\epsilon=0.05$ , learning rate $5\times 10^{-4}$ , entropy weight 0.03, expression length 4 to 64.
uDSR	Batch size 1000, GP meld population 100 for 20 generations.	add, sub, mul, div, sin, cos, exp, log, sqrt, 1.0, constants, poly token	1 batch core, reward inv NRMSE, $\epsilon=0.05$ , learning rate $5\times 10^{-4}$ , entropy weight 0.03, expression length 4 to 100, GP meld and linear poly enabled.

Appendix

Table 17: Algorithm specific configurations for the formal Core50 evaluation. Gray rows mark validation algorithms.

Model	Stage	Primary role	Affects search
Opus 5	Postprocessing	Symbolic processing and judgment	No
Llama	LLM-SR search	Equation structure generation	Yes
Llama	DrSR search	Generation and diagnostic feedback	Yes

Appendix

Table 18: Language model roles in SymbolicArena . Opus 5 operates only during evaluator side postprocessing.

Axis	Role of Opus 5
ID	None
OOD	None
SYM	Simplification and equivalence judgment
MIN	Simplification before deterministic tree size computation
EFF	None
STAB	Pairwise structural judgment of final expressions

Appendix

Table 19: Role of Opus 5 in Multi Axis Evaluation.

Explore similar work

May 12, 2026cs.SC

FePySR: A Neural Feature Extraction Framework for Efficient and Scalable Symbolic Regression

A fundamental challenge in symbolic regression (SR) is efficiently recovering complex mathematical expressions from observational data. Although this problem is NP-hard, many expressions of practical interest decompose naturally into combinations of nonlinear feature modules, concentrating structural complexity into a small number of reusable components. Here, we introduce FePySR, a two-stage framework that reduces the SR search space by extracting valid features prior to equation search. FePySR first employs a heterogeneous neural network to constrain observational data to a set of candidate expressions, then performs structural optimization within this refined expression space using PySR. Across five standard benchmarks, FePySR outperforms state-of-the-art methods by achieving higher equation recovery rates. On a set of 75 highly complex synthesized equations, FePySR recovers 36 equations, while producing substantially smaller mean squared errors on the remaining unrecovered cases, with reduced computation time compared to PySR. FePySR's first stage also maintains consistent performance under varying numbers of selected top features and increasing levels of noise in the observational data. Applied to ordinary differential equations governing biological systems, FePySR successfully identifies governing equations in 24 out of 100 tests where PySR recovers none. Taken together, FePySR is a generalizable framework that can enhance the SR solvers, enabling the efficient and reliable recovery of symbolic expressions across scientific domains.

Zhiming Yu, Wangtao Lu, Xin Lai

School of Mathematical Sciences, Zhejiang University, Hangzhou, Zhejiang, China · Faculty of Medicine and Health Technology, Tampere University, Tampere, Finland

Symbolic Regression Unified Path-Integral Update Scientific Discovery

Ask PDF

Sep 1, 2026cs.LG

Neural Symbollic Regression Using Deep Learning and Sparse Modelling

Symbolic Regression (SR) seeks to find succinct mathematical expressions that represent the fundamental relationships within data, providing interpretability and scientific understanding that exceeds that of black-box models. Nevertheless, traditional methods like Genetic Programming face challenges with scalability and are highly sensitive to noise, while sparse regression techniques such as SINDy rely significantly on predetermined feature libraries. In this work, we present a Neural Symbolic Regression (NSR) framework that treats neural networks as functional preconditioners for symbolic discovery. Our approach uses a decoupled pipeline: a neural network first learns a smooth, noise-robust approximation of the target function in an interaction- aware nonlinear feature space. LASSO is then applied to extract sparse, interpretable closed-form expressions. To improve predictive accuracy and symbolic fidelity by integrating distributed hyperparameter optimization with Ray Tune and ASHA scheduling. Experiments on the Nguyen benchmark suite show that our approach consistently outperforms SINDy and non-tuned neural baselines in RMSE, noise robustness, and out-of-distribution generalization. Ablation studies confirm the significance of feature interactions, neural depth, and tuning strategies. In general, this study presents a scalable and understandable neural-symbolic framework, creating a solid link between neural approximation and the discovery of sparse equations for scientific machine learning.

Ravi Kumar U, Sumitra S

Department of Mathematics Indian Institute of Space Science and Technology Thiruvananthapuram, India

Symbolic Regression Neuro-Symbolic Framework Deep Learning Spectral Preconditioning

Ask PDF

May 22, 2026cs.LG

When Good Equations Get Bad Scores: Improving Symbolic Regression Through Better Parameter Optimization

Symbolic Regression (SR) plays a central role in scientific knowledge discovery by distilling mathematical equations from observational data. Most existing SR methods function within a bi-level optimization framework: an outer loop that searches for the discrete equation structure, and an inner loop that optimizes the continuous parameters of that structure. Crucially, parameter-fitting quality directly determines a structure's score and thus the outer-loop search. However, nonlinear operators make the inner loop highly non-convex, and budget-driven reliance on fast local solvers (e.g., BFGS) often yields poor local minima and underestimated scores for correct structures. This ``Good Structure, Bad Score'' phenomenon becomes a key bottleneck, degrading efficiency and misguiding the search away from the true equation. To resolve this, we propose SAGE-Fit (Structure-Aware and Semantics-Guided Evaluator for Symbolic Regression), an SR-native fitting framework that exploits the dual native priors of symbolic expressions. By capitalizing on the structural and semantic priors unique to SR, we design tailored modules for each property, thereby effectively mitigating this optimization bottleneck. Extensive experiments demonstrate that our approach, as a plug-and-play module, significantly enhances evaluation fidelity and universally improves the performance of various SR systems.

Boxiao Wang, Kai Li, Zhiwei Chen +5

Institute of Automation, Chinese Academy of Sciences

Symbolic Regression Unified Path-Integral Update Structure-Aware Bad Discrete

Ask PDF