Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
Authors: Kei Sen Fong
Mehul Motani
Venue: ICML 2025 poster (2025)
Keywords: Symbolic Regression
Benchmark
Primary area: general_machine_learning->supervised_learning

Abstract: Symbolic Regression (SR) algorithms select expressions based on prediction performance while also keeping the expression lengths short to produce explainable white box models. In this context, SR algorithms can be evaluated by measuring the extent to which the expressions discovered are Pareto-optimal, in the sense of having the best R-squared score for a given expression length. This evaluation is most commonly done based on relative performance, in the sense that an SR algorithm is judged on whether it Pareto-dominates other SR algorithms selected in the analysis, without any indication on efficiency or attainable limits. In this paper, we explore absolute Pareto-optimal (APO) solutions instead, which have the optimal tradeoff between the multiple SR objectives, for 34 datasets in the widely-used SR benchmark, SRBench, by performing exhaustive search. Additionally, we include comparisons between eight numerical optimization methods. We extract, for every dataset, an APO front of expressions that can serve as a universal baseline for SR algorithms that informs researchers of the best attainable performance for selected sizes. The APO fronts provided serves as an important benchmark and performance limit for SR algorithms and is made publicly available at: https://github.com/kentridgeai/SRParetoFronts

Pareto-Optimal F Fronts for Benchmarking Symbolic Regression Algorithms
Kei Sen Fong' Mehul Motani 12
Abstract
Symbolic Regression (SR) algorithms select ex- pressions based on prediction performance while also keeping the expression lengths short to pro- duce explainable white box models. this con- text, SR algorithms can be evaluated by measur- ing the extent which the expressions discov- ered are Pareto-optimal, in the sense of having the best R-squared score for given expression length. This evaluation is most commonly done based on relative performance, in the sense that an SR algorithm is judged on whether it Pareto- dominates other SR algorithms selected in the analysis, without any indication on efficiency or attainable limits. In this paper, we explore ab- solute Pareto- Pareto-optimal (APO) solutions instead, which have the optimal tradeoff' between the mul- tiple SR objectives, for 34 datasets in the widely- used SR benchmark, SRBench, by performing ex- haustive search. Additionally, we include compar- isons between eight numerical optimization meth- ods. We extract, for every dataset, APO front of expressions that serve universal base- line for SR algorithms that informs researchers of the best attainable performance for selected sizes. The APO fronts provided serves as an important benchmark and performance limit for SR algorithms and is made publicly available at: https//github.com/eentiigege Frrrrtr
0.7
SBP-GP
AFP
0.8
GP-GOMEA
gploars
AFP
R-Squared Score
0.3
SR Algarithms APO Frant
AlFoynman
Operoe
9
13
15
220
Equation Length
Figure I: Our APO front on a subset of 34 datasets from SRBench (La Cava et al., 2021). The APO front informs us that the benchmarked SR algorithms may not have suf- ficiently explored the search space of short expressions in these datasets. Note: Results here are aggregated over the 34 datasets and we perform per-dataset analysis later.
In contrast to black-box machine learning models, SR al- gorithms produce expressions that are explainable and in- terpretable. Thus, SR has become a first-class algorithm in various fields including physics (Udrescu & Tegmark, 2020), material sciences (Wang al., 2019), engineering (Martinez-Gil & Chaves-Gonzalez, 2020) and healthcare (Christensen et al., 2022: Fong & Motani, 2024).
1. Introduction
Symbolic Regression (SR) is the task of finding closed-form analytical expressions of practical interest that describe the relationship between variables in measurement dataset.
In recent large-scale benchmarking work for SR algorithms, the performance of SR algorithms on black- box regression datasets is measured based on the relative performance of the algorithm (La Cava al., 2021: de Franca et al., 2024). In particular, refer relative Pareto-optimal SR al- gorithm one that produces expressions with the best R- squared (R2) score for expressions up to given size, when compared to other SR algorithms selected in the analysis.
Deepartment of Electrical and Computer Engineering, National University of Singapore, Singapore. "N.1 Institute for Health, Institute for Digital Medicine (WisDM). Institute of Data Sci- ence, National University of Singapore, Singapore. Correspon- dence to: Kei Sen Fong <ffngkkeisen @ www.uu......oommmmmmmme. edu>, Mchul Motani <motani edu.s".
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
In this paper, we propose to move towards absolute Pareto- optimality, where we exhaustively search through the en-

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
tire space of expressions within maximum expression size to find absolute Pareto-optimal (APO) fronts for 34 datasets used in SRBench (La Cava et al., 2021), widely used SR benchmark. In contrast to synthetic datasets that are generated from a known ground-truth equation, black- box datasets do not possess a baseline equation to compare against. Our APO fronts serve as a useful baseline for benchmarking, and they inform SR researchers about the efficiency and attainable limits of SR algorithms.
fronts and compare their performance. This provides SR researchers with large-scale evidence justify their selection of numerical optimization methods.
2. Related Work
2.1. SR Benchmarking
SR algorithms have been evaluated for their relative Pareto- dominance in SRBench (La Cava et al., 2021), an important benchmark an SR with extensive hyperparameter search, However, unlike the APO front that we propose to find, the relative Pareto front does not provide any indication of the efficiency and attainable limits of SR. Furthermore, in the relative analysis, the axes do not reflect the true R² score and true model size. Rather, the axes are the rankings of I²2 score and model size, which means that depending on the SR algorithms selected, the conclusion can vary drastically, In this paper, we discover fixed APO fronts that remain as an unchanging universal reference baseline even when new SR algorithms are developed and propose conventions that would enable SR benchmarking insights to be independent of the selection of competing SR algorithms.
To illustrate, Figure I shows the performance of SR algo- rithms from SRBench (in red) and our APO front (in green), averaged across the subset of 34 datasets. The APO front represents a fundamental limit on the performance of any SR algorithm. Figure shows that many algorithms are far from the APO front, performing much worse on both axes and so there is potential room for improvement.
To find these APO fronts, we use an SR algorithm that ex- haustively searches expressions of a fixed size. To limit the expression lengths, we adopt gene expression program- ming (GEP) (Ferreira, 2002) technique of using a genome- phenome system to manipulate and represent expressions. In GEP, the expressions are encoded in their fixed-length K- K-expression form (analogous to genome) and evaluated in their corresponding decoded variable-length mathematical expression form (analogous to phenome). Hence, GEP allows expressions to be easily modified in their fixed- length K-expression form, while simultancously allowing for varying functional complexity (Ferreira, 2002). Other strengths of K- e-expressions include not requiring a valid- ity check since all K-expressions are guaranteed to decode into a valid expression. We should, however, note that the entire search space of expressions is large and computa- tionally expensive, in the order of Od'), where d is the number of variables in the dataset and is the expression length. Consequently, our data from exhaustive searches on 34 datasets from SRBench establishes a valuable base- line. Sharing this data will prevent redundant computa- tions and enhance SRBench results with more meaning- ful insights. These assets are made publicly available at ParetooFotssssssss....
2.2. Exhaustive Search via K-expression
Gene expression programming combines the strengths from genetic algorithms and genetic programming by utilizing genome-phenome system through the introduction of K- expressions. K-expressions are strings of fixed-length which are subjected to reproduction (genome). These strings are decoded by forming a variable-length expression from a subset of the string (phenome). The decoding process is done by iterating through the K-expression and building an expression tree from top to bottom then from left to right, until no valid symbols can be added. For example, the string ** + -bcdde" is decoded as (a+bb) * (c- with the symbol being in excess and not included in the already full expres- sion. This system enables us to create fixed length strings with expressions of variable size. K-expressions are guar- anteed to decode to form a valid mathematical expression because they have a tail component, in which only terminal symbols are present (Ferreira, 2002). K-expressions have been used to develop DistiISR (Fong & Motani, 2023), an SR algorithm. However, in DistilSR, terminal symbols are replaced with weighted linear combinations of variables, which consequently exclude certain expression structures. In contrast, our algorithm is better suited to exhaustively cover a larger class of expression structures.
The main contributions of this paper are as follows:
1.
We extract an APO front of expressions from the exhaus- tive search for a subset of 34 datasets in the widely-used SR benchmark, SRBench. The APO fronts serve useful baseline for benchmarking and inform SR re- searchers about the efficiency and attainable limits of state-of-the-aa SR algorithms.
2. We propose conventions for analyzing SR benchmark results (SRBench) with the APO front. These can help avoid potentially contradictory conclusions in SRBench.
2.3. Other Exhaustive Search SR
There also exists other SR algorithms which perform ex- haustive search (Kammerer et al., 2020; Bartlett et al., 2023). However, these algorithms tend to make simplifications or
3. We empirically investigate 8 different widely-used nu- merical optimization methods in obtaining the APO
2

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
assumptions that have a different search space. In our work, our primary goal is not develop state-of-the-ar SR al- gorithm, but rather to extract APO front of expressions for black-box datasets and an exhaustive search is simply a means to this end. Thus, our algorithm design choice is targeted at using fewer assumptions to provide more robust assets that encompass a larger class of functions.
the primitive operands set form the primitive_symbolt We select two values of head length. Ideally, this value should be higher to evaluate more expressions. However, the num- ber of expressions in the exhaustive search grows exponen- tially with the head Bead_length, which sets practical constraints. Within our high budget of 1,480,000 core-compute-houve we repeated the search for 10 candom_seed (1 (11284, 1 11964, 15795, 21575, 22118. 23654, 29802. 5390, 6265, 860) for head_length = 3 and did the search for one random_seed (1 1284) for head = 4. The random seeds are the same values as the 10 used in SRBench and were used to generate initial guesses for numerical optimization for the range (-1,1). Eight numerical optimization methods (to optimize the ERCs, R) were selected to be evaluated: (i) L-BFGS-B (Liu & Nocedal, 1989; Zhu et al., 1997), (ii) conjugate gradient (CG) (Hestenes et al., 1952), (iii) Nelder- Mead (Nelder & Mead, 1965), (iv) Powell (Powell, 1977), (v) sequential least squares programming (SLSQP) (Laww son & Hanson, 1995), (vi) truncated Newton constrained (TNC) (Nash, 2000), (vii) trust-region constrained (trust- constr) (Conn et al., 2000), (viii) BFGS (Fletcher, 2000). These methods formed the numerical der these settings, we ensured that the exhaustive search produced exact replicable and reproducible results.
2.4. Numerical Optimization in SR
Numerical optimization has been used in SR to obtain nu- merical constants in expression (Kommenda et al.. 2020; Chen et al., 2015). Instead of using fixed predefined constants in expression trees, the symbol, TR, denoting Ephemeral Random Constant (ERC), is introduced. To Fletcher-
optimize the values of the ERCs, the Broyden-Fletcher- Lodfarb-Shanno (BFGS) algorithm (Fletcher, 2000) is fre- quently used for SR algorithms (Biggio et al., 2021; Petersen et al., 2019). In this paper, we empirically investigate vari- ous numerical optimization methods in obtaining the APO front for SR, beyond the commonly used ones, and compare their performance. Specifically, we consider the following additional 7 methods: (i) L-BFGS-B (Liu & Nocedal, 1989; Zhu et al., 1997), (ii) conjugate gradient (CG) (Hestenes et al., 1952), (iii) Nelder-Mead (Nelder & Mead, 1965), (iv) Powell (Powell, 1977). (v) sequential least squares program- ming (SL.QQP) (Lawson & Hanson, 1995), (vi) truncated Newton constrained (TNC) (Nash, 2000), (vii) trust-region constrained (trust-constr) (Conn et al., 2000).
cumericcl.optimizaaiii optimization list. Un-
At this point, with sufficient details already presented to the reader, we would like to make a short interjection to add an important caveat that the APO fronts are obtained with respect to a fixed primitive function set and a specific local numerical optimizer. Ideally, the fixed primitive function should include all possible function symbols, and the numer- ical optimizer should be true global optimizer that is able to find the best set of numerical constants. In practice, this is not yet possible due to the exponential relationship between the number of primitive functions selected and the lack of a true global optimization algorithm. Starting from this point in the text, it would be more accurate to add the subscript indicating the primitive set used and/or the numerical opti- mization method used (e.g., APO Sub, Mul, Div, Pow), BFGS). ). For readability, we add the subscript sparingly, when we feel it is important to remind the reader of this caveat. Finally, it should be noted that despite this caveat, we are already able to find fronts that have a large performance gap with the equations found via existing SR algorithms.
3. Data Collection and Analysis Conventions
In this section, we describe the procedure used to obtain 34 APO front! for 34 datasets in SRBench. This procedure is supplemented by code implementation, that produces the exact same data if the same random seeds are used. In our data collection, randomness is involved because the numerical optimization methods require an initial guess, which is randomly generated. Since this work uses datasets from in SRBench (La Cava et al., 2021), the value of the random seeds used are selected to match that of SRBench. In the final part of the section, we propose conventions for SR analysis that address some issues with current trends in SR benchmarking analysis.
Algorithm outlines the steps for our exhaustive search SR algorithm. At the start of the algorithm, all possible K-expressions are constructed. The head component of K-expressions is built in Step 3, taking all permutation of primitive symbols (i.e., primitive functions and operands) and each permutation is appended by a tail component (i.c., chain of terminal symbols) in Step 4 that guarantees that all K- K-expressions generated produce valid expressions (Ferreira, 2002). Here, the length of the tail is determined by h x (nmax 1) + 1, where h is the head length and Aaaa is the
3.1. Exhaustive Search Algorithm
In our experiments, we use primitive function set of representing addition, subtrac- tion, multiplication, division and power (the absolute value of the base is taken) respectively. all of which have arity two. The primitive operands comprise of the d features in the dataset, X, and the ERC, R. The primitive function set and
Audd.Sub.Mul..iv....
1 Availlable at http://gghhaa..cmm//eeigag/gaaas
3

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
Algorithm I: Exhaustive Search Pseudo Code
Input:
Primitive-symbole headhh random_seed, numeriial..pptimizaiinn. X, where X is an n x matrix
representing n datapoints with d normalized features (mean O, standard deviation 1)
Output: search_data
max arity Ma MaxOperati wwt.:opprra.tttriiitooo..e. www.iott symbols))
srminal_symbollccommt + Length (max_arity -ttrrtttttttatyyytttraaa
3 population + Permutationn hat.muttttiiiipiiiiiiieeess/asii
heagheng
x Tail component of terminal symbols ensures decodablility
population t -ppenddeerrmnnnlSSymbbl.sse (popullliinniiillliinntnninnn terminal symbol count)
5 search + reull
for Method numerical .opiiiiiaaootiiiiiiiiii Jist do
for K-expression E population do
expression + Decode(K-exprent))
** randoo is required to ensure reproducibility of initial guess of numerical parameters used in numerical optimization expression + Methoodeepressson, X, random seed) expression +
10
X)
11
search data t search data U {(expression, expression score))
12
end
13 end
14 return search data
maximum operand arity of the primitive function set. Then, algorithms. Common practices used, which we discuss below, can help to create visually informative plots. In effort to further improve these practices, particularly in Pareto analysis, we propose two conventions as potential areas of improvement.
iterate through all eight numerical optimization methods (see Step 6) and K-expressions (see Step 7) and perform numerical optimization. The optimized expression and its score on the dataset X are then stored. This forms the raw data which we make available. For head length = 3, we ran Algorithm on 34 datasets from SRBench (see Appendix for dataset details), which had the condition that there were less than 1000 datapoints and less than 10 features, and for headJengee 4, we reduced this 30 datasets by excluding datasets with more than 6 features.
3.3.1. CONVENTION #1: FOR PARETO ANALYSIS SR BENCHMARKING CONSIDER USING THE ACTUAL QUANTITY AS AXES.
In current and recent SR Pareto analysis, the axes do not reflect the true R² score and true model size but instead reflect the rankings of R² score and model size. Depending on the set of SR algorithms selected, the conclusion can vary drastically, such as gplearn's 's performance shown in Figure 2a & 2b. Additionally, ranking removes substantial important information. For example, based on Figure 2a, be very close to DSR in terms of expression length, but in terms of the actual values in Figure 2c, this is not the case, with gplearn having more than three times the length of DSR and is in fact very close in length to BSR, which could not have been inferred from Figure 2a. Further- more, Pareto analysis in other fields more commonly in terms of the actual quantity and not rankings of the quantity.
3.2. Extracting APO Front
For each dataset, extract the APO front for the dataset by taking the expressions with the highest R² score among all the random seeds and all the numerical optimization methods used. One expression for each unique expression length is selected. We set a condition that expressions of larger length must perform better than the expressions of smaller length qualify for inclusion into the APO front, otherwise no expressions are selected for that length. In Table 1, we show the exact expressions on the APO front for 9 datasets. The table for all 34 datasets is included in Appendix B. The expressions in these Tables are the most important main asset this paper.
gplearn
A phenomenon, which term Rank Inversion Paradox (see Appendix M for worked examples and elaboration), can also occur, causing contradictory conclusions when rankings are involved in Pareto analysis.
3.3. Proposed Conventions for SR Benchmarking
Current practices in SR benchmarking have been useful summarizing the relative performances of various SR
It is understood that using ranked axes can provide bet- ter spaced datapoints which can convey messages more

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
Table 1: APO Front Expressions for 9 Datasets (Full table for all 34 datasets in Appendix B).
Dataset APO Front Expressions (in increasing length, seppaaseae by "-+"
1027 ESL
1096 FacultySallriees
T
192_vineyard
30
210.cloud
228 elusage
230 machine cpu
22
485 analcatdata..ehiclle
$19_vinnie
522 pm10
Space
59K4F
0.8
0.5
0.0
FEal
R.Squree Benang
Squuree Score Numting
0.4
0.4
100
600
0.3
200
youp
2120
Length Banking
Langth Sanking
(a) Plot from ranking algorithms used in the work by La Cava et al. (2021). used in the work by Petersen et al. (2019).
(b) Plot from ranking algorithms
(c) Plot from actual values of algorithms used in the work by La Cava et al. (2021).
Figure 2: Three graphs based on the same data (dataset 579_fri_cd www.fro.....o 250_5), but with very different visual conclusions. For transferable results across papers, actual values are preferred over ranks. See Appendix M for inversion paradox with ranks.
efficiently, but the cost of potentially conveying mislead- ing conclusions is high. We hope the trend of using only be reconsidered and
3.3.2. CONVENTION #2: FOR PARETO ANALYSIS SR BENCHMARKING, CONSIDER SUPPLEMENTARY ANALYSIS WITH PER-DATASET RESULTS TO CONFIRM THE TREND.
rankings in SR Pareto analysis supplemented with analysis using actual quantities.
Finally, presenting results in actual values also has the addi- tional benefit of making the plots and conclusions transfer- able to other SR works.
In typical SR analysis, results often aggregated across datasets, while trends on individual datasets are omitted in the work. By aggregating results, too much simplifica-
5

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
0.7
0.9
10
ISPO5
0.9
(122))
0.5
Scom
102
Sooee
0.5
(A3R
0.5
14500
0.3
0.4
0.4
40
40
40
Equual Langth
Equaition Langt-
Eqquation Langth
(a) A Aggregated all datascts
(h) Dataset: 594_fri_c2.100.5
(c) Dataset: 366_fri_c1_10005
Figure 3: By aggregating results, information is lost. For example, Operon's good performance individual datasets is not communicated to researchers via the aggregated plot.
tion may have occurred, and information is lost. In Figure 3, it appears that DSR, AFP, EPLEX and SBP-GP are on the relative Pareto-optimal front from the aggregated results. However, it turns out that for the methods DSR, gplearn GP- GOMEA, Operon, AFP FE, AFP, FEAT, EPLEX. SBP-GP, FFX, AlFeynman, their respectively frequency of being the relative Pareto-optimal front across individual datasets are 100%, 67%, 63%, 50%, 37%, 33%, 30%, 27%, 17%, 10%, 3%, respectively. Note that contrary to aggregate re- sults, AFP, EPLEX and SBP-GP are not on the APO front for most datasets. To better represent results, we recom- mend providing supplementary Pareto analysis conducted on individual datasets.
Backpropagation Genetic Programming (SBP-G) (Virgolin et al., 2019), gplearn (Stephens, 2016).
4.1. Finding #1: Are current methods benchmarked in SRBench close to the APO front?
Based on Figure 1, it seems as though none of the meth- ods are close to the APO Sub, Mul. Div, Pow) front, and only SBP-GP has an R82 score that exceeds that of the APO front, but at the cost of much higher expression length. Looking at individual datasets, in Figure 4, we show examples of individual APO fronts and classify each into one of 4 types (the characteristic of each type is described in the captions). Among the 34 APO fronts, 15 fall under Type I, 5 fall un- der Type II, 5 fall under Type III and 9 fall under Type IV (see Appendix C for exact details). The performance varies largely across different datasets. Notably, for datasets that fall under Type IV, all SR algorithms fail to sufficiently search the space of small expressions, otherwise, they would have performance close to the APO front. It is, however, promising to see DSR consistently perform relatively close to the APO front consistently, as shown in the examples in Figures 4a, 4b, 4c.
4. Results and Discussion
In this section, we analyze the results obtained from our col- lected data and extracted APO fronts. We compare against the 14 SR methods benchmarked in SRBench (La Cava et al., 2021), under GPL-3.0 license: Age-Fitness Pareto Optimization (AFP) (Schmidt& & Lipson, 2010), Age-Fitness Pareto Optimization with Co-evolved Fitness Predictors (AFP FE) (Schmidt & Lipson, 2010), AlFeynman (Udrescu et al., 2020), Bayesian Symbolic Regression (BSR) (Jin et al., 2019), Deep Symbolic Regression (DSR) (Petersen al., 2019), Epsilon-Lexicase Selection (EPLEX) (La Cava et al., 2016), Feature Engincering Automation Tool (FEAT) (La Cava et al., 2018), Fast Function Extraction (FFX) (Mc- Conaghy, 2011), Genetic Programming- based Gene-pool Optimal Mixing Evolutionary Algorithm (GP-GOMEA) (Virgolin et al., 2017), Interaction-Traun Transformation Evolu- tionary Algorithm (ITEA) (de Franca & Aldeia, 2021), Multiple Regression Genetic Programming (MRGP) (Ar- naldo al., 2014), Operon (Burlacu al., 2020), Semantic
4.2. Finding #2: How rare is it to obtain an expression with performance close the APO front?
In most datasets, it is rare to obtain an expression with performance close to the APO Suh, Mul. Div, Pow) front. In Figure 5, we show selected examples of the distribution of n²2 score in the search space. These three examples (Figures 5a, 5b, 5c) were selected as their top-bin in the histogram had the minimum, median and maximum value among all other histograms. Additionally, among all the histograms, 85.3% have less than 1 within the top- in
6

---

Pareto- Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
0.9
SPP-GP FEAT
SAP-GP
0.8
GP GOMEA
GONEA
TAF AP PE FEAT
0.7
R-Squared Score
R-Squared Score
0.6
EPLEX
0.5
AFP
gpleam
0.4
[SR
0.2
MRGP
0.1
Algorithms APO Front
Algorithms APO Front
AlFeynman
=0
NRGP
11
13
17
200
9
11
15
17
220
Equation Length
Equation Length
(a) Type (Example Dataset: 611_fri, [1__fr_____10____01
(b) Type II (Example Dataset: 522 522.pm10)
0.5
0.9
0.8
0.4
0.7
DSR
GP-GOMEA
SBP-GP
R-Squared Score
0.6
AlFoynman
AFP
0.3
DSR
0.5
R-Squared
AFF.FE
AFP FE Operon
PEAT
0.4
0.2
0.3
0.2
0.1
GP-GOMEA
AlFeynman
0.1
EPLEX
SR Algorithms
SR Algorithms
6BPGF
=00
=0
NRGB
Equation Length
Equation Length
(c) Type III (Example Dataset: v78.visua))
(d) Type IV (Example Dataset: 228.elusage)
Figure 4: Examples of APO fronts for the 34 datasets. Type refers to cases where there are SR algorithms of larger length that exceed the R² score of the APO front. Type the same as Type I, but the SR algorithms do not exceed by more than 0.1 R² score difference. Type III refers cases where there no SR algorithms of larger length that exceed the R² score of the APO front, but there is least one SR algorithm that is less than 0.1 R score away. Type IV is the same as Type III, but with all SR algorithms more than 0.1 R² score away. The dotted lines are extensions of the longest APO front expression we could attain due to realistic computation constraints in exhaustively searching larger expressions.
the histogram, 66.9% have less than 0. within the top-bin in the histogram, 35.2% have less than 0.01% % within the top-bin the histogram. This tells us that simple random
sampling of expressions from the exhaustive search space of short expressions is unlikely to provide performance close to the APO front.

---

Pareto- Optimal Fronts for Benchmarking Symbolic Regression Algorithms
wabe: 0:224
70
80
Percentage of Search Space (%6
40
Percentage of Search Space (%6)
70
50
50
40
40
20
Percentage
30
20
1032N
20
susogh
Aree 525 cm 200 des
50
0-81
014-2
0.24.3
0.3-44
0.4.55
8.5..66
06:0.7
0.7-8.8
0.-4.9
091100
0-8.1
0402
02.83
0.344
0 0-0.6
0.60.77
0.7-8.8
0.8-9
0.100
6-0.1
0.102
02:03
0.304
0440.
0.5-0.6
0.600
0.7-888
0.09
039-1.
Squared Value
Squared Value
R-Squared va le
(a) Dataset: 522.pm10
(b) Dataset: 230 mpchinee.ccu
(c) Dataset: 112.chhaase.geyaee.
Figure 5: Histograms of 3 examples of the distribution of It score among the exhaustive search space of expressions. The left graph's 's top-bin contains 0.0002% (the lowest among all experiments), the middle graph's 's top bin contains 0.06% (the median) and the right graph's top-bin contains 6.61% (the most).
4.3. Finding #3: How is the loss landscape, in terms of expression structure around expressions in the APO front, like?
On average across the datasets, 23.7% of expressions that are two-mutations away from the APO front are able to greedily mutate back to the APO front. Dataset 659 "e__sseutheex177 had the lowest rate of 0% and dataset 690 visualizing galaxy had the highest rate of 68.8%. Our collected data provides us with powerful database that allows us to analyze the loss landscape in terms of expres- sion structure, which has hardly been analyzed in SR. Fur- ther metrics can be created and analyzed in depth (such as changing two-mutations to k-mutations), but this requires a thorough separate analysis and far extends beyond the scope of this paper, leave extensions of this idea future works in SR.
Unlike the loss landscape on numerical parameters, visu- alizing the loss landscape in terms of expression structure is challenging because there is no 'nice' continuous axis to represent changes in structure. To this end, we describe the loss landscape by considering if expressions that are two- mutations away from the APO front are likely to mutate back into the APO front.
To measure this, we take all expressions that are two- mutations away from the APO front and apply two greedy one-step mutations. To illustrate this, let consider the expression 3.5MM
(on the APO We generate expressions two-
4.4. Finding #4: How does the choice of numerical optimization affect results?
front for dataset mutations away structurally from this expression. One example is which then considers all one-step mutations and greed- ily picks the best to mutate into Div(Mul(xdatta Mulxxa). 1,1),, which in turn considers all one-step mutations and greedily picks the best to mutate into (ub(MMu((,,0055),MM. (0.531)), In this case, the initial expression that is two-mutations away from the APO front is able to greedily mutate back to the APO front.
To answer this, we compare the various APO fronts (i.e., comparing APO Sub, Mol, Div, Pow), Powell, APO Sub, Mul, Div. Pow), BFGS, APO Sub, Mul, Div, Pow), TNC and 5 other APO fronts from the other numerical opti- mization methods). Empirically, there is not a significant difference the results across different numerical optimiza- display
tion methods evidenced in Figure 6, where the histograms with the largest difference. In terms of KL divergence, Figures 6a & 6b showed the largest difference, followed by Figures 6a & 6c. Yet, the distribution of R² score in the search space appears largely similar despite being the ones with the largest difference. We also consider only the expressions on the APO front, and evaluate the standard deviation of the R² score of these expressions across varying numerical optimization methods. The mean standard deviation of R score ranged from a minimum of 0.0 to a maximum of 0.043, with a mean and median
Another example is which then considers all one-step mutations and greedily picks the best to mutate into Mul(Div(z3.- 240),Sub(rr ,1 180),, which in turn considers all one-step mutations and greedily picks the best to mutate into 330-://2gg....73).Subbb..
0.0416)), (.teiiicrra--668...... 0.446)))
In this case, the initial expression that is two-mutation away from the APO front unable to greedily mutate back to the APO front.
8

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
10 190
Vakue 0.59
3224
Waker
70
60
Percentage of Search Space (%6
Percentage of Searcn Space (%6
Percentage of Search Space
50
40
30
30
20
20
10
10.%%
100%
5-77
007% 1Ak
6-1
030-3
0.2..33
4304
0.4.55
03..
06607
0.708
0..0
0.910
00.11
0.102.2
0.20333
0.30.
0.40.5
05.8.66
0.6077
0,7.8.8
0.3-3.9
0.9-1.0
04,11
0.102
020.
033.4
044.5
0.8.6
0.0-7
0.7.8
0.00.
0.9-1.
R Squared Value
A R-Suurred Value
Squared Va la
(a) Optimization: Powell
(b) Optimization: BFGS
(c) Optimization: TNC
Figure 6: Histograms of the distribution of R² score across different numerical optimization methods on the dataset 192.vineyard. We chose to display this dataset and these three optimization methods because the difference in probability distributions (measured via KL divergence) is the greatest here (among all other numerical optimization and datasets).
of 0.00043 and 3.42e-19 respectively. Only 1.47% of the standard deviation measured exceeded 0.01. These observations indicate the expressions are stable with respect to the numerical optimization method selected and in terms of R² score, there is no strong reason to favor one more than another. Though the differences between the numerical optimization methods are small, for completeness of results, the percentage of times in which each method produces the expression on the APO front are 32.1%, 25.0%, 20.0%, 10.9%, 9.41 %, 2.06%, 0.588% and 0.00% for the methods Powell, CG, trust-constr, L-BFGS-B, BFGS, TNC, SLSQP and Nelder Mead respectively. The commonly used BFGS does not have clear advantage in obtaining expressions the APO front compared to the other methods.
the prediction performance of APO fronts is similar across different random seeds (the standard deviation for expres- sions across the random seeds is less than 0.1 for 86.2% of the expressions), which gives some confidence that the parameters are not often stuck at local optima, otherwise the prediction performance would have been different. We also could not include all possible operators. Finally, though the computation resources spent on the project is high, we hope that by making these data publicly available, it reduces the computational burden other researchers and allows researchers to focus new discoveries without the need for redundant computational efforts.
5. Conclusion
In this work, we move towards absolute Pareto optimal fronts, which enables SR researchers to determine the abso- lute difference between SR algorithms and the best attain- able expressions. This complements current SR benchmark- ing efforts in which the relative performance with other SR algorithms is made instead. We also propose conventions to improve upon common practices in SR benchmarking efforts to reduce the tendency for misleading conclusions to be drawn. Then, we compare among different numeri- cal optimization methods in obtaining the APO fronts to investigate the choice of numerical optimization method to use in SR. Finally, potentially help SR algorithm design. A key finding is that the search space of short, simple equations is sufficiently expressive and should be explored more before expanding the search space to longer equations, a mechanism that is related to increasing explainability, which is primary rea- for practitioners to pick SR over alternative machine learning algorithms in the first place.
4.5. Finding #5: Which datasets should SR focus on?
Based on the discussion earlier in Finding #1 and Figure 4, for datasets with results under Type III and Type IV, current SR algorithms are not sufficiently searching the space of short expressions, which SR researchers should focus on. Even for datasets with results under Type II, the performance of short expressions are comparable with SR algorithms performance, which demands investigation.
APO Sub, Mul, Div, Pow) front
4.6. Limitations and Societal Impact
report several findings that can
In the most precise and definite sense of word 'absolute", the numerical optimization method needs to provide the global minimum. However, there are no numerical optimization methods that have this property. To help mitigate this, we ran experiments over different initial guesses (set by the random seed) for numerical optimization. We found that

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
Acknowledgements
Cranmer, M., Sanchez Gonzalez, A., Battaglia, P., Xu, R., Cranmer, K., Spergel, D., and Ho, S. Discovering sym- bolic models from deep learning with inductive biases. Advances in Neural Information Processing Systems, 33: 17429-17442, 2020.
This research/project supported by the National Re- search Foundation, Singapore under its AI Singapore Pro- gramme (AISG Award No: JISG3-PDD-2023-08-05222
and A*STAR, CISCO Systems (USA) Pte. Ltd and Na- tional University of Singapore under its Cisco-NUS Ac- celerated Digital Economy Corporate Laboratory (Award 121001E0002). The computational work for this article was partially performed on resources of the National Supercom- puting Centre, Singapore (httss//www.nscc.sg))
de Franca, F. O. and Aldeia, G. Interaction- transformation evolutionary algorithm for symbolic re- gression. Evolutionary Computation, 29(3):367-390, 2021.
de Franca, F. O.. Virgolin, M., Kommenda, M., Majumder, M., Cranmer, M., Espada, G., Ingelse, L.,, Fonseca, A., Landajuela, M., Petersen, B., et al. SRBench++; Princi- pled benchmarking of symbolic regression with domain- expert interpretation. IEEE Transactions on Evolutionary Computation, 2024.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Ferreira, C. Gene expression programming in prob- lem solving, pp. 635-653. Springer London, Lon- doi: 10.1007/ 978-1-4471-0122-99 54. URL htttps://don //dgi.org/
References
don, 2002. ISBN 978-1-4471-0123-9,
Arnaddo, I., Krawiee, K., and O' D'Reilly, U.-M. Multiple regression genetic programming. In Proceedings of the 2014 Annual Conference on Genetic and Evolutionary Computation, pp. 879-886, 2014.
Fletcher, R. Practical methods of optimization. John Wiley
& Sons, 2000.
Bartlett, J., Desmond, H., and Ferreira, P. G. Exhaustive symbolic regression. IEEE Transactions on Evolutionary Computation, 2023.
Fong, K. and Motani, M. DistiISR: distilled version of gene expression programming symbolic regression. In Proceedings of the Companion Conference Genetic and Evolutionary Computation, pp. 567-570, 2023.
Biggio, L., Bendinelli, T., Neitz, A., Lucchi, A., and Paras- candolo, G. Neural symbolic regression that scales. In International Conference on Machine Learning, pp. 936- 945, PMLR, 2021.
Fong, S. and Motani, M. Explainable and privacy- preserving machine learning via domain-aware symbolic regression. In Conference on Health, Inference, and Learning, pp. 198-216. PMLR, 2024.
Burlacu, B., Kronberger, G., and Kommenda, M. Operon C++: An efficient genetic programming framework for symbolic regression. In Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion, pp. 1562-1570, 2020.
Friedman, H. Greedy function approximation: A gradient boosting machine. Annals of Statistics, pp. 1189-1232. 2001.
Chen, Q.. Xue, B., and Zhang, M. Generalisation and do- main adaptation in GP with gradient descent for symbolic regression. In 2015 IEEE Congress on Evolutionary Com-
Hestenes, M. R., Stiefel, E., et al. Methods of conjugate gradients for solving linear systems, volume 49, NBS Washington, DC, 1952.
putation (CEC), pp. 1137-1144. IEEE, 2015.
Chèze, G. and Fieux, E. The inversion paradox and ranking methods in tournaments. arXiv preprint arXiv:2503.02429 2025.
Jin, Y., Fu, W., Kang, J., Guo, J., and Guo, J. Bayesian symbolic regression. arXiv preprint arXiv: 1eXiiv/910008892 2019.
Christensen, J., Demharter, S., Machado, M., et al. Iden- tifying interactions omies data for clinical biomarker discovery using symbolic regression. Bioinformatics, 38 (15):3749-3758. 2022.
Kammerer, L., Kronberger, G., Burlacu, B., Winkler, S. M., Kommenda, M., and Affenzeller, M. Symbolic regres- sion by exhaustive search: Reducing the search space using syntactical constraints and efficient semantic struc- ture deduplication. Genetic Programming Theory and Practice XVII, pp. 79-99. Springer, 2020.
Conn, A. R., Gould, N. I., and Toint, P. L. Trust region methods. SIAM, 2000.
10

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
Kommenda, M., Beham, A., Affenzeller, M., and Kron- berger, G. Complexity measures for multi-objective sym- bolic regression. In Computer Aided Systems Theory- EUROCAST 2015: 15th International Conference, Las Palmas de Gran Canaria, Spain, February 8-13, 2015, Revised Selected Papers 15, pp. 409-416. Springer, 2015.
Petersen, B. K., Larma, M. L., Mundhenk, T. N., Santi- ago, C. P., Kim, S. K., and Kim, J. T. Deep symbolic regression: Recovering mathematical expressions from data via risk- policy gradients. The International Conference on Learning Representations, 2019.
Powell, M. J. D. Restart procedures for the conjugate gra- dient method. Mathematical Programming, 12:241-254, 1977.
Kommenda, M., Burlacu, B., Kronberger, G., and Affen- zeller, M. Parameter identification for symbolic regres- sion using nonlinear least squares. Genetic Programming and Evolvable Machines, 21(3):471-500 2020.
Schmidt, M. D. and Lipson, H. Age-fitness Pareto opti- mization. Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation, pp. 543,544, 2010.
La Cava, W., Spector, L.,, and Danai, Epsilon-lexicase selection for regression. In Proceedings of the Genetic and Evolutionary Computation Conference 2016, pp. 741- 748, 2016.
Stephens, T. Genetic programming in Python, with scikit- learn inspired API: gplearn. 2016.
La Cava, W., Singh, T. R., Taggart, J., Suri, S., and Moore, J. H. Learning concise representations for regression by evolving networks of trees. The International Conference on Learning Representations, 2018.
Udrescu, S.-M. and Tegmark, M. Al Feynman: A physics- inspired method for symbolic regression. Science Ad- vances, 6 6(16):eaay2631, 2020.
La Cava, W., Orzechowski, P., Burlacu, B., de França, F. O., Virgolin, M., Jin, Y., Kommenda, M., and Moore, J. H. Contemporary symbolic regression methods and their relative performance. Neurips Track on Datasets and Benchmarks., 2021.
Udreseu, S.-M., Tan, A., Feng, J., Neto, O., Wu, T., and Tegmark, M. AI Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Neural Information Processing Systems, 33:4860-4871, 2020.
Lawson, C. L. and Hanson, R. J. Solving least squares problems. SIAM, 1995.
Virgolin, M., Alderliesten, T., Witteveen, C., and Bosman, P. A. Scalable genetic programming by gene-pool opti- mal mixing and input-space entropy- based building-block learning. In Proceedings of the Genetic and Evolutionary Computation Conference, pp. 1041-1048, 2017.
Liu, D. C. and Nocedal, J. On the limited memory BFGS method for large scale optimization. Mathematical Pro- gramming, 45(1):503-528, 1989.
Virgolin, M., Alderliesten, T., and Bosman, P. A. Linear scaling with and within semantic backpropagation-base genetic programming for symbolic regression. In Pro- ceedings of the Genetic and Evolutionary Computation Conference, pp. 1084-1092, 2019.
Martinez-Gil, J. and Shaves-Gonzaller J. M. novel method based on symbolic regression for interpretable semantic similarity measurement. Expert Systems with Applications, 160: 13663, 2020.
Matsubara, Y., Chiba, N., Igarashi, R., and Ushiku, Y. Re- thinking symbolic regression datasets and benchmarks for scientific discovery. Journal of Data-centric Machine Learning Research, 2024.
Virgolin, M., De Lorenzo, A., Medvet, E., and Randone, F. Learning a formula of interpretability to learn inter- pretable formulas.
Parallel Problem Solving from Nature-PPSN XVI: 16th International Conference, PPSN 2020, Leiden, The Netherlands, September 5-9, 2020, Proceedings, Part II 16, pp. 79-93. Springer, 2020.
McConaghy, T. FFX: Fast, scalable, deterministic symbolic regression technology. Genetic Programming Theory and Practice IX, pp. 235-260, 2011.
Vladislavleva, E. Smits, G. F,, and Den Hertog, D. Order of nonlinearity complexity measure for models gener- ated by symbolic regression via Pareto genetic program- ming. JEEE Transactions on Evolutionary Computation, 13(2):333334, 2008
Nash, S. G. A survey of truncated-Newton methods. Journal of Computational and Applied Mathematics, 124(1-2): 45-59, 2000.
Nelder, J. A. and Mead, R. A Simplex method for func- tion minimization. The Computer Journal, 7(4):308-313 1965.
Wang, Y., Wagner, N., and Rondinelli, J. M. Symbolic regression in materials science. MRS Communications, (3):793-805, 805, 2019.

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
Zahir, S. Normalisation and rank reversals in the additive analytic hierarchy process: A new analysis. International Journal of Operational Research, 4(4):446-467, 2009.
Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound- constrained optimization. ACM Transactions on Mathe- matical Software (TOMS), 23(4):550-560 1997.
12

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
A. More Details the 34 Selected Datasets
Table Details of the subset of the 34 selected datasets from SRBench.
Dataset
Number of Datapoints
Number of Features
1027 ESL
488
Saso_aaeultySllariss
50
192_vineyard
52
2
210 210_cloud
108
5
228 elusage
55
2
230_machine_cyy
209
6
485 analcatdata_
48
519_vinnia
380
522 pm10
500
523 aneaaot
100
547_no2
500
556_analcatdataaappea22
475
3
557 analcatdata apneal
475
3
_61_cpu
209
579 779_fri.00.250.3
250
594 fri c2 100
100
596 196.ri.co wwr.fr.c...... frl.222500..
250
597 fri c2 500 5
500
5
601.fi.cc_2500
250
5
611.fri.c3 c3.00.5
100
5
613 1.3.fri_33.00 250.5
250
617 fri c3 500 5
500
024_fric0_100.55
100
631_fri.c15500.
500
649 fri e0 500 5
500
656_fri.c11100.5
100
5
659 1191
47
7
663_rabe_266
120
665 sleuth case2002
147
6
678 378.visualizing.emvie
111
687 687_sleuth.ex
62
690 1300visualizingggabbee
323
4
706 wleuth.caseeeeeeeeeeeeeeaaaaa 1202
93
6
712. chscaseec .eesser
222
2
Note: The Friedman datasets (Friedman, 2001) are synthetically generated by summing collection of smooth, bell-shaped components, where each component is given a randomly chosen scale and depends only on a small, randomly selected subset the inputs. Each component is centered at a random point in the input space and its shape is determined by picking a random orientation and applying random stretch factors along each direction. Inputs themselves are drawn from a standard normal distribution. For finer details of the process, please refer the original work by Friedman (2001).
13

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
B. APO Front Expressions for All Datasets
Table 3: APO Front Expressions for all 34 Datasets.
Dataset
APO Front Expressions (in increasing length, separated by (-))
1027.E8L
SSabbulllAddd J.0 Mull),,l.
1096 SacultySalaries
33
210.cloud
228_elusage
- 1.666-09 (0aa:0~~MMMaap000000-3-A//Addd--
+
230 machine cpu
Sy o-t.huu..............o...//a//aa
485 vesie
1,55e 3555155555lllrr000770005-55/11i0
+
519 vinmie
Divi D-
522_pm10
8.23c 10
523 janaalattaaaaneevoo
2.13e Mall- (.969))-....99-9.-.--
2D0.2000SSSSiiitllll111110000000 133,z)
547.102
$56.amalcutdaat..pmaa2
357 analcallaaaa apncal
561
579 c0
594_fn_c2_100.5
1.8e
596_fri.e..250.3
397 197.0...2.500.5
M0tp::/A.M.....SSSS....-AAASSaS0 -1.05-09
601. wwl.fn.cc..... 25..5
611_fm.c3.100
613 fri. c3. 250.5
617_fn_c3_500.3
1555 -De
- (ask
624_fri.c0_100
31
$5005
649 fri c0 00.00.0
2g
656_fn.cl_100.
4.72
659 sr9.sleuth.ex7111
2g
- DDicMMalle,...500.5u66...
663_ahe_266
30
665. 665_sleuth.cass200
178.visualizingeevvvvoo enviroomental
6. 6.88--09
687 sleuth ex1605
+ Sttpaaasee..aaaaaia.../....00000
690. 090..iismliiing galaxy
706 sleuth caselluss22222
11__chcuse_geyeerl
14

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
C. Classification of Dataset Type
Table 4: Classification of dataset type based on performance of SR algorithms with respect to the APO front.
Type
Datasets
Type I
547 no2, 556 analcatdata 2 557_ htt.aaaacattttaappppa
579 wtt_fri.c0222..011 594_fri.c2_100 596_fri.c2.2500
597 fri c2 500 5, 601 fri cl 250 5, 611 fri c3 100 5,
413_fri.c3.250.3 617_fri_c3_500. 624_fri_c0_100) 631 fri cl 500 5, 649 fri c0 500 5, 656 fri el 100
Type II
1027 ESL, 522 pm10, 561 opu, 663 rabe 266, 690 visualizing galaxy
Type III
210 ..loud $19_vinnie 523 analcatdata neavote,
678 visualizing visualizing.eevioonneeen environmentall, 112_chscase cesccsee.eesserr
Type IV
1096 FacultySallrries 192 vineyard, 228 elusage,
230 230._achine.cpu, 485 a.t.aaalaaatt..echil.... 659 359_sleuth.ex177 1714, 665 sgleuthccase2002, 687 3le_the ex1605, 706 sle.theehhcc case 1202
D. Validity of K-expressions
All K-expressions can be decoded into valid mathematical expressions. In Algorithm terminol.symboo .count is the length of the tail given by a formula, headJength is hyperparameter that is the number of operators and operands in the head of K-expression. For example, the K-expression ** + f is decoded as (a + b) (c d). where the length of the tail is 5 and the number of operators and operands in the head is 4.
In the main text, we claim the length of the tail determined by X (Amarr - 1) 1, where h is the head length and Amma the maximum operand arity of the primitive function set. The length of the tail is ensure that there are no non-terminal symbols with empty arguments. It is computed from assuming the worst-case scenario where each symbol in the head is an operator, there will be at most Amaz arguments. We have h symbols at the head, each with most maax arguments to fill, but since each of the h symbols fill up an empty spot, each symbol only creates at most (Nmar 1) empty spots to fill. This leads to h X (Amax 1) potential spots for the tail symbols to fill. Including the initial empty spot to be filled, we get x (Amar - 1) 1.
E. Extension to SRSD Metrics
SRSD (Matsubara al., 2024) introduced standard for metrics to in SR analysis. For their accuracy metric, this is simply R²2 > 0.999 and can be obtained from the Ra column. Of greater interest are the solution rate and NED introduced by SRSD. However, these are only computable on datasets with closed-for foom ground-truth, so they are not directly applicable to this work. Nonetheless, meaningful way to incorporate these metrics is by taking the APO equations we found (see Appendix B) and treating them as 'proxy ground-truth', enabling a new metric to assess SR algorithms performance on box datasets.
F. Verification on Physics Experiments
We also ran experiments on data from the Newtonian dynamics experiments by Cranmer et al. (2020) to demonstrate its applicability to recovering already verified equations for verification. Specifically, we applied the APO front on the data to recover the known force laws via the representations learned by the internal 'message function'. We used the 1-D version of orbital force equation and charged particle force equation in Cranmer et al. (2020). The ground truth equation structures appear on the APO front we find using our approach, and these equations are the true underlying physical law.
15

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
G. Expanded Function Set
To perform a similar analysis that is extended to a much larger function set (i.e., {Add, Sub, Mul, Div, Pow, Sin, Cos, Aresin, Arcos, Exp, Log, Max, Min)), our sampling of run-time estimates a 12.7x x increase in compute resources required. However, we recognize the interest in having some indication of expansion of the function set and have done so on single random seed 860, and using only BFGS optimization in the supplementary materials.
H. Regularization in SR
An alternative approach that we could have used is to add regularization and train on the test set, but it is not clear which regularization scheme or regularization parameters should be used to obtain the best upper bound. Thus, we opted to exclude regularization (e.g.. length, complexity, parameter magnitudes) in obtaining the Pareto front.
From the supplementary materials given in the publicly available GitHub link, the relation between train and test R is generally strong with little trade-off. However, the gap (e.g., train R² test R²) is more interesting and the supplementary materials data show some trade-off between the gap and model length, which warrants a more detailed analysis in future work.
I. More Details on Two-mutation Analysis
In finding #3, are interested in studying properties of the loss landscape that tell the tendency of getting stuck local optima and the difficulties in assessing the global optima. For typical machine learning, for loss landscapes similar to the Rastrigin function, in regions close to the global optimum, the function has many small depressions or "basins of attraction' (valleys) that can trap an optimization algorithm in local minima. Although the global minimum exists, the surrounding parameter space is filled with numerous local minima that resemble shallow valleys. For SR, the loss landscape of greater interest is on the function structure rather than on the numerical parameters and 'considering if expressions that are two mutations away from the APO front are likely to mutate back into the APO front' is our way of assessing if the equations on the APO front are surrounded by many 'basins of attraction' that makes tougher for SR algorithms to discover them.
J. Operator Frequency Analysis
We can also analyze the frequency of operators appearing on the APO front to gain insights which can possibly inform hyperparameter settings in SR. For example, in evolutionary SR, a operator that appears with higher frequency on the APO front may benefit from having higher chance of being selected with parts of an expression has be filled at random (e.g.. I-mutation). From our data in the supplementary materials, for equations of length 9, the average frequency of Add, Sub, Mul, Div, Pow is 0.7000, 0.8625, 0.9500, 0.4833, 1.0042, respectively. Note that the frequency can exceed because there can be more than 1 of the same operators in each equation.
K. Suggestions for New APO Front-based Metrics
The APO plots allow users to come up with a large range of interesting new measures. A non-exhaustive list of suggested measures include: i). R² closeness, by measuring the percentage of datasets where the R² value obtained by the SR algorithm is within 0.1 (i.e., 0.1 vertical distance on the Pareto plot) of the extended APO front (plateau after max length searched), ii). Euclidian distance from single algorithm with the front, though the relative magnitudes of axes need to be decided, iii). compare the entire field of SR with respect to the front with existing measures like hypervolume (HV). For the algorithms AFP, AFP FE, AlFeynman, BSR, DSR, EPLEX, FEAT, FFX, GOMEA, ITEA, MRGP, Operon, SBP-GP, gplearn, the values for R² closeness are 62% 59%, 12% 15% 62%, 68%, 68%, 35% 65%, 47% 32% 62%, 68%, 59%, respectively.
16

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
L. Extended Plots
0.7
6.3
0.6
0.
Sa--r)
0.5
29776
0.7
TIEA
0.1
0.0
-02
0.1
0.1
0.3
02
033
Spare
0.
4.9
0.
-10
7020
-2.
-10
300
570 055
600 Equattin Longes
900
300
600
Equation Length
Equation Length
(a) Aggregated on all datasets
(b) Dataset: 594_fri.c2.100.5
(c) Dataset: 650
Figure 7: Plots from Figure 3 but with reduced truncation. Unlike rankings, using the actual metric values may be tricky to use in visualizing the results since the range of values can be large with largely varying distribution.
M. Worked Examples of Rank Inversion Paradox
Below are worked examples with exact values to illustrate, with concrete example with real values, examples of when using rankings in Pareto analysis can lead to potentially contradictory results. We show an example where the aggregating method is to take the median of the ranks and another where the aggregating method is to take the average of the ranks.
M.1. Using Median of Ranks
Problem Setup:
Consider 4 SR algorithms: A, B, C D and 3 datasets, Ds1, Ds2, Ds3.
On Ds1, the raw performances are:
R2 (higher is better)- C: 0.9, B: 0.8, A: 0.7, D: 0.6
Model size (lower is better) B: 3, A: 5, D:
On Ds2:
R C: 0.9, D: 0.8, A: 0.7, B: 0.6
Model size- A: B: 5, C: D:9
On Ds3:
R B: 0.9, A: 0.8, C: 0.7, 0.6
Model size- D: 3, A: 5, B: C:
Now, using the procedure in SRBench to determine which algorithms are relatively Pareto optimal (using SRBench github and consistent with Fig. 2 of the NeurlPS paper): first take the rank per dataset, then take the median ranks.
Case 1: Using only Algorithm A, B, C in Pareto analysis
Here, we use the notation (3, 2) to represent that on Ds1, Algorithm A is Rank 3 in R2, and Rank 2 in model size.
On Dsl:
A: (Rank 3 in R Rank 2 in model size), B: (2, 1), C: (1,
On Ds2:
17

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
A: (2, 1), (3. C: (1.3)
On Ds3:
A: (2, I), B: (1, 2), C: (3, 3)
Median rank across Ds1, Ds2, Ds3:
A: (2, D) (2, 2). C: (1,3)
Thus, one would conclude Algorithm & C relatively Pareto optimal.
Case 2: Using Algorithm A, B, C & D in Pareto analysis
On Dsl:
A: (3, 2), B: (2, 1), C: (1, (4, (4,44
On Ds2:
A: (3, B: 2), C: (1, 3), D: 42.
On Ds3:
A: (2, 2), B: (1, 3), C: (3, 4), D: (4, I)
Median rank across Ds1, Ds2, Ds3:
A: (3, 2), (2, 2), C: (1, 3), (4.
Thus, would conclude Algorithm B & C are relatively Pareto optimal. Note that although only Algorithm D was added, Algorithm B is now suddenly optimal, and Algorithm A is suddenly not optimal. We can call this the 'Rank Inversion Paradox', inspired by (Chèze & Fieux, 2025; Zahir, 2009). Note that in taking actual quantities, any Algorithm that is not Pareto optimal will never be Pareto optimal with the addition of new Algorithms(s).
M.2. Using Average of Ranks
Problem Setup:
Consider 4 SR algorithms: A, B, C D and 3 datasets, Ds1, Ds2, Ds3.
On Ds1, the raw performances are:
R2 (higher is better)- A: 0.9, 0.8, C: 0.7, 0.6
Model size (lower is better) C: 3, A: 5, D: 7, B:
On Ds2:
R² A: 0.9, B: 0.8, C: 0.7, D: 0.6
Model size B: 3, C: 5, A: D:9
On Ds3:
R2 C: 0.9, D: 0.8, B: 0.7, A: 0.6
Model size B: 3, C: 5, 7, 9
Now, we first take the rank per dataset, then take the average ranks.
Case I: Using only Algorithm A, B, Cn Pareto analysis
Here, use the notation (1,2) represent that Ds1, Algorithm A is Rank 1 in R², and Rank 2 in model size.
On Dsl:
A: (Rank in R², Rank 2 in model size), B: (2, 3), C: (3, 1)
On Ds2:
18

---

Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms
A: (I.3), B: 1), (3.2)
On Ds3:
A: (3, 3), B: (2, 1), C: (1, 2)
Average rank across Ds1, Ds2, Ds3:
A: (1.67, 2.67), (2, 1.67), C: (2.33, 1.67)
Thus, one would conclude Algorithm & are relatively Pareto optimal.
Case 2: Using Algorithm A, B, C & D in Pareto analysis
On Dsl:
A: (1, 2), B: (2, 4), (3, (4.
On Ds2:
A: (I, 1), (3,
On Ds3:
A: (4, 3), B: (3, 1), C: (1, 2), D: (2,
Average rank across Ds1, Ds2, Ds3:
A: (2, 2.67), (2.33, 2), C: (2.33, 1.67), 33.33,3.677
Thus, would conclude Algorithm A & C are relatively Pareto optimal. Note that although only Algorithm D was added, Algorithm C is now suddenly optimal, and Algorithm B is suddenly not optimal. We can call this the 'Rank Inversion Paradox', inspired by (Chèze & Fieux, 2025; Zahir, 2009). Note that in taking actual quantities, any algorithm that is not Pareto optimal will never be Pareto optimal with the addition of new algorithms(s).
N. Other Complexity Measures
In the supplementary materials, also include other complexity measures such the count of operators, the count of numerical constants, Kommenda's complexity (Kommenda et al., 2015), Virgolin's trained linear elastic net interpretability estimator (we use their trained rescaled coefficients) (Virgolin et al., 2020) and Vladislavleva's 's order of non-linearity (using =1c-6 as done in their work) (Vladislavieva et al., 2008). We have also added Peterson's 's complexity that is used in DSR (Petersen et al., 2019).
19


Reviews:
[{"review_id": "CwAzCJu5Dz", "reviewer": "Reviewer_NEMV", "summary": "This paper proposes the absolute Pareto optimal (APO) front as a new benchmarking methodology for evaluating symbolic regression (SR) algorithms. Conventional SR evaluation is based on relative Pareto dominance with respect to other algorithms, but this does not provide any measure of efficiency or achievable limits. The authors established the theoretical limits of the trade-off between expression length and prediction performance (R-squared) through exhaustive search on 34 real-world datasets in SRBench. Specifically, they generated all possible expressions (within a specific length limit) using the K-expression format, and applied eight different numerical optimization methods to find the optimal parameters. As a result, the difference in performance between the current SR algorithm and the APO front was clarified, and it was shown that many algorithms were unable to find potential optimal solutions, especially in the search for short expressions. The paper also proposes new conventions for SR benchmark analysis (using actual values, and performing individual analysis rather than just aggregating results across datasets). The results of this research, which required a huge amount of computing resources (approximately 1.48 million core compute hours), are published as a valuable baseline for future SR algorithm development.\n\n## Update After Rebuttal\n\nAfter careful consideration of the authors' rebuttal, I am changing my evaluation from \"weak reject (2/5)\" to \"weak accept (3/5)\". The authors have sufficiently addressed my main concerns and promised the following improvements:\n\n1. Clear explanation of the limitations of \"absolute Pareto optimality\" and improved notation (A.1)\n1. Addition of complexity metrics beyond expression length (Kommenda, Virgolin, Vladislavleva, etc.) (A.3, C.2)\n1. Provision of performance evaluations on both training and test data (A.3)\n1. Self-contained explanation including proofs of K-expression properties (A.4)\n1. Validation of the practical value of the APO front in modeling physical phenomena (A.11)\n1. Provision of additional analyses such as function frequency analysis (A.14)\n\nWhat particularly influenced my evaluation is that, as other Reviewers (especially Gtam) have emphasized, this research represents an important contribution to the evaluation of SR algorithms and provides a valuable resource to the research community despite its high computational cost. The improvements indicated by the authors address the limitations of the paper and enhance the practical utility of the results.\n\nI also understand the concerns raised by other reviewers (particularly Reviewer GBr3) regarding the analysis conventions, but the authors' rebuttal (especially D.13-D.19) provides compelling examples of why traditional ranking-based analysis is problematic. On this point, while recognizing that there are various perspectives, I believe there is value in the practices proposed by the authors.\n\nOverall, I am confident that this paper makes a valuable contribution to SR research, and with the implementation of the improvements promised by the authors, its value will be further enhanced. Therefore, I support the acceptance of this paper.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> General\n\nPlease refer to point (G.1) in the response to Reviewer Gtam.\n\n> Claims And Evidence\n\n(A.1) We thank the reviewers for advice on our discussion on the subtleties with the term “absolute” in the limitations section that includes our mitigation strategies. To further improve the clarity, we will denote APO with a subscript that states the primitive function set used and the numerical optimization method (e.g., ${APO}_{(+,-,*,/,**,sin,cos,...),SLSQP}$) as a caveat. Despite not having a true global numerical optimizer, we are already able to find fronts that have a large performance gap with the equations found via existing SR algorithms.\n\n> Methods And Evaluation\n\n(A.2) For complexity indicators, please refer to point (C.2) in the response to Reviewer Gtam.\n\n(A.3) We implement the reviewers’ suggestions and include more data for SR research. In the $\\texttt{Extracted\\\\_APO\\\\_Fronts}$ folder, instead of having only one column for $\\texttt{EquationParameters}$, we now have 3 columns, each consisting of a vector of numerical constants optimized based on all data, train data only and test data only, respectively. For each of the 3 vectors of numerical constants, we also include their performance on all data, train data only and test data only. These will be made publicly available in the next upload opportunity. Below is an example (numbers truncated to 4 s.f.) of the improvement (e.g., $\\texttt{1027\\\\_ESL\\\\_BFGS\\\\_3\\\\_860\\\\_summary.csv}$):\n\n|…|EquationStructure|EquationParametersAll|R2FitAllEvalAll|R2FitAllEvalTrain|R2FitAllEvalTest|…|\n|-|-|-|-|-|-|-|\n|...|\n|…|$\\texttt{Sub(Mul(xdata[3],x[0]),Mul(xdata[1],x[1]))}$|[0.5503,-0.5307]|0.8056|0.8009|0.8172|…|[0.5713,-0.5118]|0.8052|0.8013|0.8136| …|[0.4931,-0.5881]|0.8016|0.7941|0.8216|...|\n\n|…|EquationParametersTrain|R2FitTrainEvalAll|R2FitTrainEvalTrain|R2FitTrainEvalTest|… |EquationParametersTest|R2FitTestEvalAll|R2FitTestEvalTrain|R2FitTestEvalTest|...|\n|-|-|-|-|-|-|-|-|-|-|-|\n|...|\n|…|[0.5713,-0.5118]|0.8052|0.8013|0.8136| …|[0.4931,-0.5881]|0.8016|0.7941|0.8216|...|\n\n> Theoretical Claims\n\n(A.4) We will make the paper more self-contained by including proofs of the properties from Ferreira in Appendix D.\n\n> Experimental Designs, Supplementary, Broader Literature\n\n(A.5) We thank the reviewer for the validation.\n\n> Essential References\n\n(A.6) We thank the reviewer and will add all SRSD metrics. For accuracy metric, this is simply $R^2>0.999$. For solution rate and NED, these are only computable on datasets with closed-form ground-truth, so not available. However, we thought that a meaningful way to incorporate these metrics is by taking the APO equations we found and treating them as “proxy ground-truth”, enabling a new metric to assess SR algorithms performance on black-box datasets. We thank the reviewer for inspiring this additional use of the APO equations and will add this and the discussion of other less commonly-used but high quality datasets in literature (e.g., SRSD) in Appendix E.\n\n> Other Strengths And Weaknesses\n\n(A.7) For “Limited versatility”, please refer to point (C.3) for more results and discussion.\n\n(A.8) For “Computational cost issues”, to the best of our knowledge, there is no empirical or theoretical alternative which could create a similarly universal baseline, so we felt the benefits far outweighed the cost. Our cost reducing strategy is to open-source these files so that it is a one-off cost instead.\n\n(A.9) For “Lack of practical guidance”, we will add guidance in Section 4 similar to point (G.1) para 3.\n\n(A.10) For “Interpretability definition\", please refer to point (C.2) where we add other complexity metrics.\n\n> Other Comments\n\n(A.11) We thank the reviewer for the creative and impactful suggestion. In the revision, we will add Appendix F which shows that in [1] in the Newtonian dynamics experiments where SR is applied on the internal functions, the equation structures discovered by Cranmer et al. appear on the APO front we find using our approach, and these equations are the true underlying physical law.\n\n[1] Cranmer, M., et al. (2020). Discovering symbolic models from deep learning with inductive biases. NeurIPS.\n\n(A.12) For APO front on a wider range of functions, please refer to point (C.3).\n\n(A.13) We thank the reviewer for the idea. Using the data in point (A.3), we will include results of train vs test $R^2$ in Appendix H. In most cases, there is no trade-off, so we also include results of train minus test $R^2$ against length.\n\n(A.14) We thank the reviewer again and add function frequency analysis in Appendix J. For example, for equations of length 9, the average frequency of Add, Sub, Mul, Div, Pow is 0.7000, 0.8625, 0.9500, 0.4833, 1.0042, respectively, which can possibly inform the settings when generating the candidate solutions in SR algorithms.\n\n(A.15) We plan to explore using this approach for various permutations of common activation functions in small neural network architectures as a baseline for NAS algorithms.", "claims_and_evidence": "The arguments in this paper are well supported by the evidence presented. First, the argument regarding the construction of the APO front is supported by data obtained as a result of exhaustive search using a vast amount of computing resources. The authors explain the search method in detail and present it as a reproducible algorithm.\n\nThe authors' claims regarding the comparison of the current SR algorithm and APO fronts are supported by detailed analysis of 34 datasets, which are classified into four types to clearly explain the trends. In particular, Figure 4 visually shows representative examples of each type to support the claims.\n\nThe claims regarding the comparison of numerical optimization methods are also supported by specific data, such as the measurement of differences in distribution using KL divergence and statistics on the ratio of equations generated on the APO front.\nHowever, there are some limitations to the claim of “absolute” Pareto optimality. As the authors themselves acknowledge, there is no guarantee that the numerical optimization methods used can find the true global optimum solution, and in this respect, it cannot be said to be truly “absolute”. However, efforts have been made to mitigate this problem by using multiple random number seeds and eight different optimization methods.", "methods_and_evaluation_criteria": "The proposed method and evaluation criteria are appropriate for the SR benchmarking problem. Exhaustive search using fixed-length expressions with K-expression provides a reasonable trade-off between limiting the search space and guaranteeing the generation of valid formulas.\n\nThe Pareto front of expression length and R-squared values as evaluation criteria directly correspond to the essential goal of SR (balancing predictive performance and interpretability). In addition, the evaluation criteria proposed by the authors (using actual values rather than rankings, and performing individual analysis as well as dataset aggregation) appropriately point out the problems with current SR benchmarking.\n\nHowever, the fact that they do not take into account model complexity indicators other than formula length (e.g. the type of mathematical operations or the number of numerical constants) and that they only consider the R-squared value of the training data and do not evaluate generalization performance are methodological limitations. However, the authors are aware of these limitations, and given the constraints on computing resources, this is a reasonable choice.", "theoretical_claims": "This paper is mainly an empirical study, and there are few theoretical assertions that require formal mathematical proof. The assertions regarding the properties of K-expressions (such as the fact that all K-expressions can be decoded into valid mathematical expressions) are based on the cited literature (Ferreira, 2002) and are not directly proven within the paper.\n\nThe claim about the complexity of the exhaustive search space (O(d^l), where d is the number of variables in the dataset and l is the length of the expression) is based on basic combinatorics and is considered to be correct.\n\nThe core claim that the APO front represents the performance limit of any SR algorithm is logically derived from the nature of exhaustive search, but is subject to the aforementioned constraint of the limits of numerical optimization.", "experimental_designs_or_analyses": "The experimental design is sound overall, and the results are valid. The authors describe the experimental settings in detail, including the selection of random number seeds, the set of primitive functions used, and the numerical optimization method, in order to ensure reproducibility.\n\nThe method for extracting APO fronts is clearly defined, and the procedure of selecting the formula with the highest R-squared score for each formula length for each dataset is logical. The condition that longer formulas must perform better than shorter formulas is also reasonable.\n\nOf particular note is that the authors compared eight different numerical optimization methods and showed that there was little difference in their performance. This provides valuable insight into the selection of numerical optimization methods in SR research.\n\nThe limitations of the experiment include the restriction on the length of the formula (head length = 3 and 4 only), the restriction on the function set used (only 5 binary operations), and the selection of a specific data set (less than 1000 data points, less than 10 features), but these are reasonable choices given the constraints on computing resources.", "relation_to_broader_scientific_literature": "In terms of the evaluation method for symbolic regression, this study extends SRBench (La Cava et al., 2021). While SRBench was based on relative Pareto dominance, this study introduces the concept of absolute Pareto optimality.\n\nThe K-expression-based approach is based on Ferreira (2002)'s genome-phenome system as an approach that combines the advantages of genetic algorithms and genetic programming. It has also been shown to be related to other K-expression-based algorithms, such as DistilSR (Fong & Motani, 2023).\n\nThe comparison of numerical optimization methods contributes to the research on the optimization of numerical constants in SR (Kommenda et al., 2020; Chen et al., 2015). In particular, the fact that it compares methods other than the BFGS algorithm (frequently used in Biggio et al., 2021; Petersen et al., 2019) on a large scale is new.\n\nThis study also mentions SR applications in various fields, such as physics (Udrescu & Tegmark, 2020), materials science (Wang et al., 2019), engineering (Martinez-Gil & Chaves-Gonzalez, 2020), and healthcare (Christensen et al., 2022), showing its relevance to a wide range of scientific literature.", "essential_references_not_discussed": "The paper cites and discusses a certain amount of relevant major literature. The major previous studies are appropriately cited in each category of symbolic regression benchmarking (SRBench), utilization of K-expression, exhaustive search SR, and numerical optimization in SR.\n\nAs a recent development in SR benchmarking, the paper also mentions SRBench++ (de Franca et al., 2024), but on the other hand, it does not follow benchmarks that propose new datasets and evaluation metrics, such as SRSD (Matsubara et al., 2024).", "other_strengths_and_weaknesses": "**Strengths:**\n\n・Originality: The introduction of the concept of absolute Pareto optimality into the evaluation of SR algorithms can be seen as an innovation in the evaluation paradigm.\n\n・Practicality: By providing a baseline in the form of the APO front, the authors have given concrete targets for future SR algorithm development. This is particularly valuable for algorithm development that focuses on the search for short formulas.\n\n・Contribution to the research community: By publishing results that require a huge amount of computing resources, the authors have reduced the computational burden on other researchers.\n\nThorough analysis: The SR is examined from multiple angles, including a comparison of eight numerical optimization methods and an analysis of the loss landscape around the APO front.\n\n**Weaknesses:**\n\n- Limited versatility: The applicability to problems with more complex relationships is limited due to the set of functions used and the restriction on the length of the expression.\n\n- Computational cost issues: The proposed method requires a large amount of computational resources, making it difficult to extend to new datasets and function sets.\n\n・Lack of practical guidance: Although the gap between the APO front and the current algorithm is shown, there are few specific algorithm design guidelines to fill the gap.\n\n・Interpretability definition: Only the formula length is used as an interpretability indicator, and other aspects (such as conceptual simplicity) are not considered.", "other_comments_or_suggestions": "- By verifying the formula used to construct the APO front in actual application examples (e.g. modeling physical phenomena), it may be possible to further demonstrate its practical value.\n\n- It may be beneficial to conduct research on an extended version of the APO front that includes a wider range of functions (e.g. trigonometric functions, exponential functions, etc.).\n\n- By reformulating the problem as a multi-objective optimization that includes R-squared values for test data as well as training data, it may be possible to construct an APO front that also takes into account generalization performance.\n\n- By performing feature analysis (such as what structures and functions are frequently used) on the equations on the APO front, it is possible to gain further insights for designing efficient SR algorithms.\n\n- The methods used in this research have the potential to be applied to benchmarking interpretable machine learning methods other than symbolic regression, and it is also worth exploring this direction.", "questions_for_authors": "I look forward to your response to the concerns I have raised above, but I have no further questions.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nVe5LT06Tb", "reviewer": "Reviewer_YJ4h", "summary": "This paper proposes an absolute evaluation criterion for symbolic regression, namely the Absolute Pareto Optimality (APO). At the same time, the effects of eight different optimization algorithms are analyzed. The establishment of this criterion is of great significance, as it provides a fair and reliable benchmark for researchers to evaluate the performance of symbolic regression algorithms, and solves the problem that in previous comparisons of algorithms, only relative indicators such as the recovery rate and R-squared value could be used for evaluation. The exhaustive search that consumed 1,480,000 supercomputer core-compute-hours demonstrates the author's determination. In this way, future comparisons of various methods will be more fair. It is expected to address the problem where one paper claims that its own method is excellent, but the results are different in another paper.\n\n## update after rebuttal\n\nI agree with the assertion from other reviewers about the attractiveness of the proposed method, and I also consider that there are some problems in this paper. The responses from authors address most of my concerns, even though they are not quite satisfactory. Since I acknowledge the value of the idea proposed by this paper, I am inclined to accept it and will keep my rating.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> General\n\nPlease refer to point (G.1) in the response to Reviewer Gtam.\n> Claims, Methods and Evaluation, Experimental Designs, Supplementary, Broader Literature, References\n\n(B.1) We thank the reviewer for the validation.\n\n> Other Strengths And Weaknesses\n\n(B.2) We thank the reviewer for the reference to problems in using ranking in other fields such as multiagent learning evaluations. We will add this to the revision to link it to broader literature.\n\n(B.3) For “absoluteness”, please refer to point (A.1) in the response to Reviewer NEMV.\n\n> Other Comments\n\n(B.4) We thank the reviewer and will refine the sentence.\n\n> Questions For Authors\n\n(B.5) We used multiple random seeds and also multiple numerical optimization methods. We can guarantee all possible structures (within a certain complexity) are searched, but to the best of our knowledge, there are no theoretical or empirical guarantees in literature for the numerical optimization methods for all the different structures, so these are the best results we can produce with current technology and knowledge. We will be clearer on this in the revision as outlined in (A.1).\n\n(B.6) We thank the reviewer for the question. Because the main reasons to select SR in the first place is its explainable and interpretable models via short concise equations, when handling datasets with higher dimensions, feature selection can be performed first. Feature selection is also already a common preprocessing step in SR algorithms since SR algorithms do not scale well with dimensions. However, the best feature selection method to use for SR is highly debatable, hence we focused on datasets with low-dimensions. SRBench is already a compilation of datasets from multiple environments and can be said to be the de facto benchmark for SR, with extensively-tuned performance of SR algorithms, we did not use other datasets.\n\n(B.7) For improving existing algorithms, please see point (G.1) para 3. Additionally, using loss landscape analysis (“Finding #3”) but using k-mutations instead of 2-mutations, and finding the best k, can be incorporated into SR algorithms, where greedy k-mutations are performed for every candidate (i.e., for every candidate structure, greedily search the best k-mutations). Note that the second suggestion is an ongoing future work which needs more thorough analysis, but is only enabled by having the APO fronts.\n\n(B.8) max_arity is the number of arguments of the function with the most arguments, terminal_symbol_count is the length of the tail given by a formula, head_length is a hyperparameter that is the number of operators and operands in the head of a K-expression. We will include a longer explanation with simple examples of K-expression as well as proofs of the properties of K-expressions by Ferreira (2002) in Appendix D.\n\n(B.9) The length of the tail is to ensure that there are no non-terminal symbols with empty arguments. It is computed from assuming the worst-case scenario where each symbol in the head is an operator, there will be at most $n_{max}$ arguments. Since the symbol itself would fill an empty argument as well (with the exception of the first symbol), we have $h \\times (n_{max}-1)+1$. We will include this in Appendix D.\n\n(B.10) We thank the reviewer for the question about Figure 5 in “Finding #2”. We did this because we wanted to show a variety of rarities of equations in the top-bin, so for each unique combination of i). random seed, ii). optimization method and iii). dataset, we made a histogram. The 3 in Figure 5 “were selected as their top-bin in the histogram had the minimum, median and maximum value among all other histograms”. As we learn later in “Finding #4”, no optimization method provides a clear prediction performance advantage. To communicate our findings more effectively in the revision, we will fix the histograms in Figure 5 to BFGS only, so as not the conflate the message of “Finding #2” with “Finding #5”. \n\n(B.11) For “Finding #3”, we are interested in studying properties of the loss landscape that tell us the tendency of getting stuck at local optima and the difficulties in assessing the global optima. For typical machine learning, for loss landscapes similar to the Rastrigin function, in regions close to the global optimum, the function has many small depressions or \"basins of attraction\" (valleys) that can trap an optimization algorithm in local minima. Although the global minimum exists, the surrounding parameter space is filled with numerous local minima that resemble shallow valleys. For SR, the loss landscape of greater interest is on the function structure rather than on the numerical parameters and “considering if expressions that are two mutations away from the APO front are likely to mutate back into the APO front”, is our way of assessing if the equations on the APO front are surrounded by many \"basins of attraction\" that makes it tougher for SR algorithms to discover them. We will include this discussion in Appendix I.", "claims_and_evidence": "The authors claim three main contributions in this paper, and there are related evidences to support these claims:\n\n1.They establish absolute Pareto-optimal fronts for 34 real-world datasets by exhaustively searching expressions up to fixed sizes using gene expression programming and K-expressions. This is rigorously supported by their computational effort and empirical comparisons showing gaps in current SR algorithms. The publicly released APO expressions and performance metrics provide a concrete, reproducible baseline.\n\n2. The authors propose standardized conventions for SR benchmarking, advocating the use of actual metric values instead of rankings. This is validated through examples shown in Figure 2 and Figure 3, demonstrating how ranked axes distort conclusions. The proposal's logic aligns with the need for transferable and interpretable benchmarks.\n\n3. The authors conduct an empirical comparison of numerical optimization methods for SR, showing minimal impact on APO front quality. Quantitative results and low $R^2$ variance across methods robustly support this claim.", "methods_and_evaluation_criteria": "The method proposed in this paper is important to the field of symbolic regression and its benchmarking practices.\n\n1.By constructing absolute Pareto-optimal fronts through exhaustive search, the authors provide an objective, domain-agnostic baseline for SR algorithms. Unlike relative Pareto fronts, APO fronts define the theoretical performance ceiling for any SR method, enabling researchers to quantify how close their algorithms are to the \"best possible\" expressions for real-world datasets. This addresses a critical gap in SR benchmarking, where prior evaluations lacked universal reference points.\n\n2.The APO fronts reveal that state-of-the-art SR algorithms systematically underperform on short expressions. This insight directs the SR community to prioritize improving compactness-aware search strategies, balancing accuracy and interpretability.\n\n3.The proposal to use actual metrics ($R^2$, expression length) instead of rankings mitigates biases introduced by algorithm selection and enhances reproducibility. For example, Figure 2 demonstrates how ranked axes can misleadingly compress performance differences, while actual values expose true gaps. This standardization fosters fairer comparisons and accelerates progress by aligning the community on shared evaluation criteria.", "theoretical_claims": "I do not find theoretical claims in this paper.", "experimental_designs_or_analyses": "The experimental design and result analysis in this paper are methodologically sound and rigorously structured. The authors construct absolute Pareto-optimal fronts through exhaustive search using gene expression programming and K-expressions across 34 real-world datasets. This approach aligns with multi-objective optimization principles and leverages prior validated techniques. Testing eight numerical optimization methods further ensures robustness, mirroring engineering practices where Pareto solutions are evaluated across diverse algorithms.\n\nThe experimental results are also analyzed rigorously. (1) Datasets are classified into four types, with some exposing insufficient exploration of compact expressions, which is a finding consistent with Pareto front \"inaccessibility\" observed in optimization studies. (2) Stability tests reveal minimal impact of numerical methods on APO front quality, demonstrating experimental control over confounding variables.", "relation_to_broader_scientific_literature": "The paper’s contributions are related to broader scientific literature:\n\n1.The absolute Pareto-optimal front concept extends classical Pareto optimality principles, widely applied in fields like chip design and biomedical modeling. By exhaustively generating APO fronts for symbolic regression, the work mirrors deterministic Pareto methods while adapting them to data-driven modeling challenges, thereby bridging optimization theory and interpretable machine learning.\n\n2.The empirical comparison of numerical optimization methods aligns with studies evaluating multi-objective algorithms, such as orthogonal evolutionary strategies or learning automata. Findings that numerical methods minimally impact APO quality resonate with literature emphasizing Pareto front stability under parameter variations, reinforcing the validity of SR algorithm evaluations.", "essential_references_not_discussed": "I do not recognize missing essential references that are not discussed.", "other_strengths_and_weaknesses": "Strengths\n\n1.This paper builds on Pareto-optimality principles, proposing a framework for generating absolute Pareto fronts through exhaustive search. This aligns with multi-objective optimization methodologies. The use of gene expression programming and K-expressions ensures valid mathematical structures, reflecting prior work in symbolic regression and genetic algorithms.\n\n2. The public release of all expressions, optimization parameters, and performance data reduces redundant computational efforts. The proposal to use actual metrics ($R^2$, expression length) over rankings addresses biases seen in multiagent learning evaluations.\n\nWeakness:\n\nWhile multi-seed experiments mitigate local optima risks, the “absoluteness” of APO fronts remains contingent on numerical optimization convergence. It is a limitation also noted in game-theoretic learning, where Pareto outcomes depend on assumptions like non-deceptive opponent behavior.", "other_comments_or_suggestions": "I have a minor suggestion to refine the expression. In the abstract, there is a sentence \"serves as an important benchmark that serves as a performance limit\" --> “serves as an important benchmark and performance limit\".", "questions_for_authors": "1.How do you address potential limitations in achieving global optimality due to local minima in numerical optimization, despite using multiple random seeds? Are there theoretical or empirical guarantees that the APO front represents the global Pareto-optimal set for the expression spaces explored?\n\n2.The APO fronts are derived from 34 SRBench datasets. How do you ensure these results generalize to domains with higher dimensionality, dynamic environments, or datasets outside SRBench?\n\n3.How should SR researchers leverage the APO front to improve existing algorithms? For instance, would integrating APO-based heuristics enhance search efficiency?\n\n4.Please explain the following terms: max_arity, terminal_symbol_count, head_length.\n\n5.Why the length of the tail is determined by $h\\times (n_{max}-1)+1$？\n\n6.In Figure 5, why do different optimization methods need to be used to obtain optimal results for different data sets?\n\n7.In Finding #3, why do you say you can describe the loss landscape by considering if expressions that are two mutations away from the APO front are likely to mutate back into the APO front?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "W0eAyQeH6e", "reviewer": "Reviewer_Gtam", "summary": "One common way to evaluate symbolic regression (SR) algorithms is to judge whether one method Pareto-dominates other SR algorithms. That means it has better performance for a given expression length. This paper proposes to evaluate SR methods against absolute Pareto-optimal solutions instead. It finds an absolute Pareto-optimal front of expressions for 34 real-world datasets from SRBench, a widely used SR benchmark. This is achieved by exhaustive search over all possible expressions of a given length and over eight different numerical optimization methods. The main contribution is a new baseline for benchmarking SR that informs SR researchers about the achievable limits of SR algorithms. Additionally, the paper proposes conventions for analyzing SR benchmark results and discusses several findings from the experiments.\n\n## update after rebuttal\nThe current paper has certain limitations:\n- not enough operators/functions included in the search space\n- relatively small maximum length of the equations\n- not all datasets in SRBench considered\n\nHowever, I also do understand that some choices needed to be made to limit the already large computational resources needed.\n\nEven in light of these limitations, I support the acceptance of this paper.\n\nIt is the first paper to provide some APO fronts, and I believe these can already be useful to researchers. From what I understand, the authors will release all evaluated equations, their fitness and complexity measured with respect to different metrics. These would be very interesting to analyze, especially for the datasets where there is a clear gap between the current methods and the found APO front. Very often, equations are only useful if they are \"simple\" enough to interpret. This paper shows that in some cases it is possible to find shorter equations with much better performance than the current methods (Figures 4c, 4d).", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> General\n\n(G.1) We thank all 4 reviewers for their comprehensive reviews, identifying strengths of the work (e.g., “the authors have given concrete targets for future SR algorithm development”, “addresses a critical gap in SR benchmarking, where prior evaluations lacked universal reference point”) while providing actionable recommendations.\n\nOur personal motivation for this work began after receiving questions from SR researchers such as “is it still worth it to develop better SR algorithms?” and “is there still meaningful room for improvements in SR algorithms or has performance saturated?” in the context of new SR algorithms’ performance on SRBench being somewhat incremental (on average ~0.0077 $R^2$ improvements). Previously, one could only answer those questions subjectively based on trends, in the same spirit that Moore's Law does with speculating transistor density. With this work, we now confidently know that there are many cases in SRBench where there is still a large potential improvement. This is in spite of some computational limits we had to make given our already-large resource budget. Thus, this work is a data-driven justification that continuing research into developing new SR algorithms is still worth it.\n\nThough the main contribution is the Pareto front equations (i.e., given explicitly as equations in Appendix B and available in supplementary materials), the results also provide insights on how SR algorithms can achieve these performances. Particularly, it shows that search space of short, simple equations is sufficiently expressive and should be explored more in SR algorithms before expanding the search space to longer equations – a mechanism that is also related to increasing explainability and interpretability, which are the primary reasons for practitioners to pick SR over alternative machine learning algorithms in the first place.\n\nAmong the various strengths, if you i). find that the Pareto fronts contributed in this paper are highly relevant and applicable to many, if not all, SR algorithms’ evaluation and ii). think that our work introduces and makes it convenient to add an informative and important baseline that you would like new SR algorithms and research papers to adopt, we hope that you could consider helping us advocate for acceptance of this work. This work makes publicly accessible, to communities of various funding levels, an important computationally-expensive baseline.\n\nWe address the main concerns in the individual replies and hope that in light of these, the reviewers would consider increasing their recommendation.\n\n> Claims And Evidence, Methods And Evaluation, Experimental Designs, Supplementary\n\n(C.1) We thank the reviewer for the validation.\n\n> Broader Literature\n\n(C.2) We thank the reviewers for the suggestion to include more complexity measures to increase utility. Our provided data allows the computation of these with negligible costs as they are derived quantities from the obtained equation. In the csv files, we have extracted, for each row in the $\\texttt{EquationStructure}$ column i). the count of operators, ii). the count of numerical constants, iii). Kommenda’s complexity, iv). Virgolin’s trained linear elastic net interpretability estimator (we use their trained rescaled coefficients), v). Vladislavleva’s order of non-linearity (using $\\epsilon$=1e-6 as done in their work). We have also added Peterson’s complexity that is used in DSR. The files will be updated and the code for this processing will be made publicly available in the next upload opportunity. We are also committed to adding new metrics, and will update the files when new works such as [1] are publicly available. Below is an example (numbers truncated to 4 s.f.) of the improvement (e.g., $\\texttt{1027\\\\_ESL\\\\_BFGS\\\\_3\\\\_860\\\\_summary.csv}$):\n\n|EquationLength|EquationStructure|...|OperatorCount|ConstantsCount|Kommenda|Virgolin|Vladislavleva|Peterson|...|\n|-|-|-|-|-|-|-|-|-|-|\n|...|\n|7|$\\texttt{Sub(Mul(xdata[3],x[0]),Mul(xdata[1],x[1]))}$|...|3|2|6|0.7620|2|7|...|\n\n[1] Kacprzyk, K., & van der Schaar, M. (2025). Beyond Size-Based Metrics: Measuring Task-Specific Complexity in Symbolic Regression. In AISTATS.\n\n> Other Strengths And Weaknesses\n\n(C.3) We thank the reviewers for the suggestions for expansion of the function set. To perform a similar analysis, our sampling of run-time estimates a 12.7X increase in compute resources required. However, we recognize the interest in having some indication of expansion of the function set and have done so on a single random seed 860, and using only BFGS optimization, which will be included in supplementary materials available to readers and discussed in Appendix G. However, for the expansion of length, even a small increase in length led to an estimated 84.0X increase cost, so we are unable to provide that even with a single random seed and only BFGS. We will make these limitations clearer in the limitations section.", "claims_and_evidence": "Claims in the papers are supported by clear and convincing evidence. The observations made about the current practices in SR literature are justified and their limitations discussed.", "methods_and_evaluation_criteria": "The paper proposes a new baseline for a well-known and appropriate benchmarking dataset - SRBench. In particular, it focuses on the real-world datasets in SRBench which do not have a known ground truth. Thus establishing Pareto-optimal equations for these datasets is important.", "theoretical_claims": "N/A", "experimental_designs_or_analyses": "The experimental designs and analyses seem valid. The paper describes the search space of equations and exhaustively searches through it. It also checks eight different numerical optimization algorithms and runs them for different seeds to ensure that they do not tend to be stuck in some suboptimal local minima. The random seeds are described and chosen to match the ones used in SRBench.", "relation_to_broader_scientific_literature": "The paper addresses an important problem in SR benchmarking. Although the number of terms is currently the main measure used to judge the interpretability of an expression, it would be beneficial to briefly discuss other approaches to measuring equation's complexity (Vladislavleva et al., 2009; Vanneschi et al., 2010; Kommenda et al., 2015; Virgolin et al., 2020; Virgolin at al., 2021).\n\nKommenda, M., Beham, A., Affenzeller, M., and Kronberger, G. (2015). Complexity Measures for Multiobjective Symbolic Regression. \n\nVanneschi, L., Castelli, M., and Silva, S. (2010). Measuring bloat, overfitting and functional complexity in genetic programming. \n\nVirgolin, M., De Lorenzo, A., Medvet, E., and Randone, F. (2020). Learning a Formula of Interpretability to Learn Interpretable Formulas. \n\nVirgolin, M., De Lorenzo, A., Randone, F., Medvet,  E., and Wahde, M. (2021). Model learning with per\nsonalized interpretability estimation (ML-PIE).\n\nVladislavleva, E. J., Smits, G. F., and den Hertog, D. (2009). Order of Nonlinearity as a Complexity Measure for Models Generated by Symbolic Regression via Pareto Genetic Programming.", "essential_references_not_discussed": "All essential references are discussed.", "other_strengths_and_weaknesses": "I think the paper contributes a very valuable asset to the community. Not only does it allow us to better judge the performance of SR algorithms in absolute terms, but it also gives us insights into the complexity of different datasets. I also believe that the dataset of all evaluated equations may provide additional insights, similar to the ones already present in the paper.\n\nTable 1 (and Table 3) shows the APO front expressions for multiple datasets from SRBench. These are very interesting, especially in conjunction with Figure 4, which shows that for certain datasets, none of the tested methods is on the absolute Pareto front. This demonstrates a large gap in the capabilities of the current algorithms and hopefully will stimulate further research. Having \"ground truth\" equations for real-world datasets will allow for a better evaluation of how far we are from closing this gap. \n\nThe main weakness of the paper lies in a still relatively limited search space. I am fully aware that the current search space already requires extensive computational resources. However, it would be beneficial to understand what the Pareto front looks like for bigger equation lengths and for a richer function set (e.g., containing trigonometric functions or exponentials). It is possible that some of the methods are, in fact, on the absolute Pareto front for some datasets - we just do not know what this front looks like for larger equations. For instance, GP-GOMEA could, in theory, be optimal for datasets in Figures 4a and 4b. This is a significant limitation for datasets where the performance of the current algorithms is much higher than the maximum performance on the found absolute Pareto front (e.g., Figure 4a). I believe the limited search space (and thus unknown parts of the absolute Pareto front) should be emphasized in the limitations section.", "other_comments_or_suggestions": "None", "questions_for_authors": "None", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1eZgKe41zN", "reviewer": "Reviewer_GBr3", "summary": "The authors aim to generate a set of \"absolute pareto optimal\" model results for 34 of the datasets in SRBench, as a way of having an upper bound on the performance of symbolic models on those datasets up to a given equation length. They exhaustively search for equations for those datasets up to a given length and use different numerical optimizers to tune the constants. They then compare to SRBench results and argue for a number of analysis conventions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> General\n\nPlease refer to (G.1).\n> Claims\n\n(D.1) We thank the reviewer for the suggestions for “Convention #1”. We will improve the clarity of our discussion by explicitly stating that our comments refer only to Pareto fronts/optimality and not statistical tests in (Demšar, 2006).\n\nThe reviewer noted that “One shouldn’t compare rankings for different experiments or to assess magnitude differences”, and this is precisely the issue with the SRBench paper and accompanying code (in SRBench github master/postprocessing/blackbox_results.ipynb and Figure 2 of their NeurIPS paper) which is the motivation for our convention. We felt it was urgent to highlight this point, given its widespread adoption in SR papers. We appreciate the reviewer’s agreement on this matter and will incorporate the quoted phrasing to enhance clarity.\n\nFurthermore, in other fields where Pareto fronts/optimality is used, e.g., economics, engineering, operations, multi-objective optimization, raw performance metrics serves as the axes. To our knowledge, SRBench’s analysis is one of, if not the only, to deviate from this practice by using rankings. Other reasons to use raw performance for Pareto fronts/optimality is to ensure the results are “transferable” and comparable across different studies, which other reviewers highlighted. We apologize for the previous wording that suggested a disagreement with the reviewer. We will add these clarifications to ensure that our position–fully aligned with the reviewer–is accurately conveyed.\n\n(D.2) For “Convention #2”, the aggregation affects both raw quantities and ranking, both convention #1 and #2 needs to be applied simultaneously. We will include median plots in Appendix L for completeness. Per the convention, we actually don’t recommend either the mean or median. Finally, Figure 1 is meant for illustration purposes, which we will address in convention #2 that unaggregated results are required for a holistic assessment.\n> Methods And Evaluation\n\n(D.3) Actually, SRBench does not use a consistent operator set for SR algorithms according to their paper, code and results file (i.e., $\\texttt{.feather}$ file). We found {Add,Sub,Mul,Div,Pow} to be commonly included in most SR algorithms. Thus, we used it so that distance away from the APO front can be attributed to search inefficiency rather than a lack of certain functions. We will make this clearer in the paper. However, we still create some new assets for the expanded function set (please see (C.3))\n\n(D.4) We thank the reviewer and will add the full untruncated plots in Appendix L, the ordering of the SR algorithms still remains the same under truncation, unlike taking average of the ranks as the axis for SRBench (which we recommend against), which can change the ordering depending on the algorithms selected.\n\n(D.5) We will include the data sampling process for the Friedman datasets in Appendix A.\n\n(D.6) In order to obtain the best upper bound on the performance with the best possible equation in the search space, we train on the full datasets. The alternative approach of using only the train set is problematic because we have to add regularization to the objective and it is not clear which to use to obtain the best upper bound. However, to provide more measurements to SR researchers we have included data on all combinations of optimization procedures (please refer to point (A.3) in the response to Reviewer NEMV). A meaningful future work enabled with the data is to explore which regularization techniques would create the closest performance to this upper bound on performance. We will include this discussion with plots in Appendix H.\n> Supplementary, Broader Literature, Essential References\n\n(D.7) We thank the reviewer for the validation.\n> Other Strengths And Weaknesses:\n\n(D.8) For the 4 numbered weaknesses please refer to points 1. (D.3), 2. (D.6), 3. (D.1) & (D.2), 4. (B.6), respectively.\n> Comments\n\n(D.9) For conventions, please refer to (D.1) & (D.2).\n\n(D.10) For our improvements in limitations, please see (A.1) & (C.3).\n\n(D.11) The APO plots allow users to come up with a large range of measures. A non-exhaustive list of suggested measures include: i). R-squared closeness, by measuring the percentage of datasets where the R-squared value obtained by the SR algorithm is within 0.1 (i.e, 0.1 vertical distance on the Pareto plot) of the extended APO front (plateau after max length searched), ii). Euclidian distance from a single algorithm with the front, though the relative magnitudes of axes need to be decided, iii). compare the entire field of SR with respect to the front with existing measures like IGD and HV. For the algorithms AFP, AFP_FE, AIFeynman, BSR, DSR, EPLEX, FEAT, FFX, GP-GOMEA, ITEA, MRGP, Operon, SBP-GP, gplearn, the values for i). are 62%, 59%, 12%, 15%, 62%, 68%, 68%, 35%, 65%, 47%, 32%, 62%, 68%, 59%, respectively. We will include the full results in Appendix K.\n> Questions\n\n(D.12) For both questions, please refer to point (D.11).", "claims_and_evidence": "> Convention #1: For SR benchmarking, axes in performance trade-offs plots should be in terms of the actual quantity and not in terms of ranking.\n> \n\nRaw performance comparisons and performance rankings achieve different goals and answer different questions. \n\nRaw performance metrics answer the question “how much different is the performance of alg A from alg B on dataset X?” and work well when comparing methods on a single dataset like in Figure 2. (Even there, though, the collapse of equation lengths in 2c) hides the fact that SBGP, ITEA, FFX, MRGP, and AI-Feynman produce models of different complexity, which is made apparent when comparing rankings in 2a). However, they bias the comparison of algs across datasets difficult because the distribution of values varies among problems.\n\nIn contrast, performance rankings are the more robust way to answer “how often does alg A outperform alg B over (many) datasets?” and that is primary question multi-dataset benchmarks like SRBench try to answer. It is also why rankings are the accepted standard for statistical comparisons of algorithms over multiple datasets. (e.g., Demšar, J. Statistical Comparisons of Classifiers over Multiple Data Sets. *Journal of Machine Learning Research* **2006**, *7* (Jan), 1–30.)\n\nThe problems the authors point out are that 1) rankings change when the set of algorithms change; and 2) rankings lose the magnitude of the changes. These points are true but they are, IMO, well-known. One shouldn’t compare rankings for different experiments or to assess magnitude differences; they are used because they are useful for benchmarking over many datasets and answering whether algorithm performances differ or not.\n\n> Convention #2: Aggregating results across datasets can have contradictory conclusions, it should be supplemented with analysis on individual datasets to confirm the trend.\n> \n\nIronically, the problem illustrated in Figure 3a)-3c) is caused by aggregating the *actual quantities* rather than the rankings, which contradicts convention #1 recommendation. \n\nIt also appears (judging by operon R2) that the authors are using the mean R2 over datasets, when they should use the median. R2 is unbounded on the negative end, so median scores would  give a more robust estimate of performance across datasets. An even more robust way to estimate relative performances of algorithms across multiple datasets would be to use rankings, which are not dataset-dependent like the R2 scores and equation lengths are. \n\nThis problem (mean aggregation of raw performance scores) also biases the results in Figure 1.", "methods_and_evaluation_criteria": "> In our experiments, we use a primitive function set of {Add,Sub,Mul,Div,Pow}\n> \n\nThe operator set the authors used doesn’t match the set used in developing/benchmarking SR algorithms on SRBench. La Cava et al 2021 / SRBench allowed the operator set {+, −, ∗, /, sin, cos, arcsin, arccos, exp, log, pow, max, min}. So, strictly speaking, the absolute pareto optimal (APO) fronts generated by the authors are not pareto-optimal with respect to the task definition in the SRBench results. I would expect the APO front to be heavily influenced by the chosen operator set, and the exhaustive space of equations to be much larger than what the authors have generated. So it is not clear to me that the APO results are directly comparable to SRBench without some large caveats.\n\nIn section 3, the authors argue against aggregating results, but then in all results figures, results with equation lengths greater than 20 (sometimes 100) are collapsed together. In some cases this presentation skews the Pareto front interpretation since most of the benchmarked methods produce results that are $\\geq$ 20. Isn’t this expressly against the recommendation?\n\n\npretty much all of the Type 1 datasets are Friedman datasets, which are synthetically generated from ground-truth analytical equations. (ref: Jerome H Friedman. Greedy function approximation: A gradient boosting machine. Annals of statistics, pages 1189–1232, 2001.) if the goal is to produce absolute pareto optimal front for those equations, presumably you would want to include the ground truth equation.\n\nthe authors don’t make it clear how they are handling train/test splits to compare to SRBench. if they are just training (model parameters) on the full datasets, the APO models will be biased relative to ones that could be possibly obtained if they followed train/test procedure of the benchmarked methods. i.e., even if a “optimal” model in some sense exists, it is not a given that one could possibly find it via finite sample training.", "theoretical_claims": "n/a", "experimental_designs_or_analyses": "See methods & evaluation criteria.", "relation_to_broader_scientific_literature": "The contribution of the paper builds on SRBench and subsequent papers that use this benchmarking resource. Its main contribution is to specify an upper front of performance for a subset of the datasets in that benchmark. to my knowledge this hasn't been suggested previously.", "essential_references_not_discussed": "Prior work is sufficiently discussed.", "other_strengths_and_weaknesses": "Strengths: Having a set of pareto-optimal models for SRBench is, in general, a good idea, as it sets a measurable goal for these datasets. \n\nWeaknesses: I think the authors work has some serious caveats that cause it to fall a bit short of delivering on that goal, including: \n\n1) restriction of operator space in generating equations\n2) apparent optimization of model constants on test data, which likely overestimates performance potential of APO \n2) proposed analysis conventions which are misguided IMO \n3) coverage of only a fraction of the datasets in SRBench", "other_comments_or_suggestions": "I would suggest the authors totally remove the analysis conventions suggestions. I think they are misguided and lack nuance. \n\nI would also suggest the authors improve their discussion of limitations to incorporate the points above. \n\nFinally I encourage the authors to discuss what types of measures an APO could provide within SRBench.", "questions_for_authors": "- Is there a measure like distance to the optimal Pareto front that could viably be used? \n\n- what is the likelihood that some of the APO models actually appear in the pareto fronts of population-based SR methods, and aren't selected because of finite sample limitations?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This submission received four reviews: 3 weak accept and 1 accept. The four reviewers all joined the internal post-rebuttal discussion.

While some of them are on the fence about the decision, they are inclined to accept this work. They value its analysis and attractive proposal of the absolute pareto-optimal fronts for the SRBench datasets. I sometimes find recent SR papers with ranking-based evaluations that a little bit misleading, specifically for top-ranked methods for each category, and I would like to welcome the proposed conventions.

Based on the internal discussion and reviews, I recommend that ICML 2025 accepts this work.
At the same time, I share non-trivial concerns from the reviewers that I hope the authors address in their camera-ready or their future studies:

> inability of the method to scale much larger than it currently does

> I would prefer to see more operators/functions included in the search space, a larger maximum length of the equations, and more datasets. 

> it is not easy for other users to build upon the work. The authors answer my question about how to use their work to improve existing algorithms in (B.7), which is not helpful enough.

> The paper has several limitations including constraints on "absolute" Pareto optimality due to restricted operator sets, limitations on equation length, and coverage of only a subset of SRBench datasets. These limitations may restrict the applicability of the findings to broader symbolic regression contexts.

> coverage of only a fraction of the datasets in SRBench

> optimization of model constants on test data, which likely overestimates performance potential of APO

Lastly, I want the authors to make sure that the camera-ready has the following changes for future studies in the research community that follow this work.

- Publish the data used in this study (including analysis) and add the URL to the camera-ready as promised in the paper since some of the reviewers and I expect the assets to be published
- Publish code with clearer, self-contained instructions to reproduce the reported results. Based on the current materials, a few of the reviewers find it difficult for other users to build on this work, which should be resolved prior to the publication as this work proposes new conventions
- Make claims a little bit more specific, especially when referring to SRBench. Since this work uses a fraction of SRBench datasets, it should be clarified when the authors make claims involving the datasets. E.g., Figure 1 should explicitly say 'using 34 of 122 black-box datasets from SRBench" instead of "Our APO front with SR benchmark (SRBench (La Cava et al., 2021)) results". This is not limited to Figure 1 caption, and the authors should carefully check and modify the claims regarding SRBench accordingly.
Decision: Accept (poster)