Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large $p$
Authors: Shengbin Ye
Meng Li
Venue: ICML 2025 poster (2025)
Keywords: Symbolic regression
Nonparametric variable selection
Extreme-scale datasets
Primary area: general_machine_learning->supervised_learning
TL;DR: PAN+SR is a scalable symbolic regression framework that uses nonparametric variable selection to efficiently handle high-dimensional datasets, improving performance across multiple SR methods.

Abstract: Symbolic regression (SR) is a powerful technique for discovering symbolic expressions that characterize nonlinear relationships in data, gaining increasing attention for its interpretability, compactness, and robustness. However, existing SR methods do not scale to datasets with a large number of input variables (referred to as extreme-scale SR), which is common in modern scientific applications. This "large $p$'' setting, often accompanied by measurement error, leads to slow performance of SR methods and overly complex expressions that are difficult to interpret. To address this scalability challenge, we propose a method called PAN+SR, which combines a key idea of ab initio nonparametric variable selection with SR to efficiently pre-screen large input spaces and reduce search complexity while maintaining accuracy. The use of nonparametric methods eliminates model misspecification, supporting a strategy called parametric-assisted nonparametric (PAN). We also extend SRBench, an open-source benchmarking platform, by incorporating high-dimensional regression problems with various signal-to-noise ratios. Our results demonstrate that PAN+SR consistently enhances the performance of 19 contemporary SR methods, enabling several to achieve state-of-the-art performance on these challenging datasets.

Ab Initio Nonparametrie Variable Selection for
Scalable Symbolic Regression with Large p
Shengbin Ye" 12 Meng Li
Abstract
the 1600s when Johannes Kepler used astronomical data to discover that Mars' orbit was elliptical. This discovery, along with Kepler's other parsimonious and analytically tractable laws of planctary motion, helped launch a scientific revolution.
Symbolic regression (SR) is a powerful technique for discovering symbolic expressions that char- acterize nonlinear relationships in data, gaining increasing attention for its interpretabilityy, com- pactness, and robustness. However, existing SR methods do not scale to datasets with a large num- ber of imput variables (referred to as extreme-scale SR), which is common in modern scientific ap- plications. This "large p" setting, often accom- panied by measurement error, leads to slow per- formance of SR methods and overly complex ex- pressions that are difficult to interpret. To address this scalability challenge, we propose a method called PAN+SR, which combines a key idea of ab initio nonparametric variable selection with SR to efficiently pre-screen large input spaces and reduce search complexity while maintain- ing accuracy. The use of nonparametric methods eliminates model misspecification, supporting a strategy called parametric-ase assistee nonparametric (PAN). We also extend SRBench, an open-source benchmarking platform, by incorporating high- dimensional regression problems with various signal-to- signal-to-noise ratios. Our results demonstrate that PAN+SR consistently enhances the perfor- mance of 19 contemporary SR methods, enabling several to achieve state-of-the-ar performance on these challenging datasets.
With the recent progress in theoretical modeling and ex- perimental instrumentation, researchers have entered new era of big data. The development of SR models is partic- ularly important, as they have emerged as a powerful tool for developing machine learning models that are intelligible, interpretable, and compact. Unlike large numerical models, the mathematical expressions used SR models enable an easy understanding of their behavior, making them valuable in fields such as physics, where they can connect newly dis- covered physical laws with theory to facilitate subsequent theoretical developments (Wu & Tegmark, 2019). Moreover, SR models offer safe and responsible option for machine learning applications with high societal stakes, such as those related to human lives, as they are well-suited for human interpretability and in-depth analysis. As such, SR models have found successful applications across a range of fields, including astrophysics (Lemos et al., 2023), chemistry and materials science (Hernandez et al., 2019; Liu et al., 2020; 2022), control (Derner et al., 2020), economics (Verstyuk & Douglas, 2022), mechanical engineering (Kronberger et al., 2018), medicine (Virgolin et al., 2020), and space explo- ration (Märtens Izzo, 2022), among others (Matsubara et al., 2024).
SR literature has traditionally focused datasets with low- dimensional inputs, often with I 10, and primarily con- sidered only relevant variables- -those used in the ground truth (La Cava et al., 2021; Kamienny et al., 2022; Sho- jace et al., 2023; Tenachi et al., 2023; Li et al., 2024). In these settings, variable selection has not been critical, as SR has largely been viewed as an optimization problem under low-noise conditions. However, modern scientific applica- tions increasingly involve datasets with far larger numbers of variables (p- = 102 to 459 in this work), often including irrelevant variables, rendering variable selection a critical yet underexplored concept in SR pipelines.
1. Introduction
Symbolic regression (SR) is mathematical technique for finding symbolic expression that matches data from an unknown function. An early example of SR dates back to
'Department of Statistics, Rice University, Houston, TX, USA 'Department of Statisties and Data Science, Northwest- em University, Evanston, IL, USA. Correspondence to: Meng Li <meng rice.edu>.
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
While variable selection a well-established topic in statis- tics, its adoption in SR has been limited and its effective-

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
ness in SR remains unclear. Existing approaches, such as random forest (RF)-based pre-selection in PySR (Cran- mer, 2023), have demonstrated limited utility. Indeed, the PySR documentation explicitly notes that options like se l ect.____ feauurss are rarely used, suggesting that cur- rent methods are not well-suited to SR tasks. This observa- tion further supported by our analysis in Appendix D.2, where RF is shown to perform unsatisfactorily. The limited performance of off-the-shelf methods like RF highlights the unique challenges of variable selection in the context of SR. Unlike typical variable selection tasks, SR variable selection demands a near-zero false negative rate (FNR), as excluding even a single relevant variable from the search space pre- vents the recovery of the true underlying function. While false positives (FPs) primarily increase computational bur- den, they do not fundamentally impede the discovery of the underlying model. This asymmetry in performance require- ments explains why standard methods often fall short and underscores the importance of designing variable selection methods specifically tailored to SR.
features, X (m1, r) in this case, and set of math- ematical operations, e.g., O {+; ,,, +;
exp, e 2). Even for this low-dimensional problem, it has been shown that exploring all expressions F(XX, O), induced by X and O is NP-hard (Virgolin & Pissis, 2022). Hence, typical SR algorithms only traverse through a small subset of the full search space, such as limiting the complexity of the can- didate SR models, total runtime, number of mathematical operations, etc.
In realistic scientific applications. particularly in the era of big data, scientists often include as many intuitively reason- able features as possible, many of which may be irrelevant to the target y. This practice causes the search space F(X, O) to expand exponneexponentiallly quick (Ye et al., 2024), mak- ing it extremely challenging-if not impossible-to recover fo(-) using algorithmic approaches alone. To this end, we propose the PAN+SR framework, which integrates the non- parametric module of PAN as model- pre-screening step. This framework excludes irrelevant features prior to applying SR methods, thereby mitigating the explosion of the search space in high-dimensional problems. Here, we assume that a high-dimensional SR problem in (I) can be reduced to
In this paper, we introduce a versatile framework, PAN+SR,
for improving SR methods at extreme scales. PAN+SR leverages the Parametric Assisted by Nonparametries (PAN) strategy (Ye et al., 2024) for an ab initio screening of large influx of input variables before expression synthesis, en- abling SR tasks at extreme scales. In light of the unique challenge of SR pre-screening, we propose a novel non- parametric variable selection method designed to minimize FN; we refer to this method as PAN throughout this pa- per. Furthermore, to evaluate PAN+SR extreme scales, we extend the open-source SR benchmarking database, SRBench (La Cava et al., 2021), with high-dimensional problems containing white noise at various signal- -oiiie ratios. In Section 6, we showcase the performance up- lift of 19 contemporary SR methods under PAN+SR. The PAN+SR framework is available as an open-source project at
y/ = (o( + Ex, for i i= 1,.,..,,0 n,
(2)
where only small subset So of po =Sol < p of fea- tures exert influence on y. Then the oracle search space F F(Xs., O) is a significantly smaller subspace of the full search space F(X, O). Thus, the successful identification of Sp. or at least a superset of So. is critical for reduc- ing high-dimensional SR problems into manageable low- dimensional ones. With this reduction, the dataset (y. Xs.) becomes sufficient for discovering fo(-), enabling SR meth- ods to handle high-dimensional problems without requiring any modifications to their algorithms.
3. Related Work
2. Background and Motivation
SRBench (La Cava et al., 2021) is a reproducible and open- source henchmarking platform for SR that has made signifi- cant strides in the field through its curation of 122 real-world datasets and 130 ground-truth problems and its comprehen- sive evaluations of 14 contemporary SR methods. SRBench has quickly gained adaptions with numerous studies leverag- ing it to evaluate accuracy, exact solution rate, and solution complexity (Kamienny et al., 2022: Landajuela et al., 2022: Kamienny et al., 2023; Keren et al., 2023; Shojace et al., 2023; Makke & Chawla, 2024). Despite its widespread use, SRBench primarily focuses on low-dimensional prob- lems, which limits its applicability in the context of high- dimensional problems, a hallmark of the era of big data. In particular, the 130 ground-truth problems from the Feynman Symbolic Regression Database (Udrescu & Tegmark, 2020)
Given a dataset (y, X) with target y E R² and features X = (#1..... xp)p E RTXP, SR assumes the existence of an analytical data geteraennnating function that links X to y:
y: = fo( fo(aa....,, Tip) + t for i= n. (1)
in the presence of observation noise EG. The goal of SR is to recover the unknown regression function fo(-) symbolically. For example, consider regressing the gravitational force be- tween two objects, F, on their masses (m) m2) and the dis- tance between their centers (r). An SR algorithm would ide- ally discovr the Newton's 's Law of Universal Gravitation, F = 6.6743 x 10 10-11, 11, mymp/r?. This is typically done by randomly constructing mathematical expressions using the
2

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
and the ODE-Strogatz repository (Strogatz, 2015) contain only the oracle features Xs. with at most p 9 features. This low and narrow dimensional scope leaves SRBench less suited for analyzing SR at extreme scales, underscoring the need for a high-dimensional SR database.
tion fo(.) by sum of regression trees,
M
Σ Z() (Pp)+ tE, EN N,,,,0,00P,,)) (3)
where each regression tree (.......... Fpp) partitions the feature space based on the values of 21,.... Zp. For each posterior sample, we calculate the proportion of splits in the ensemble (3) that use 2 as the splitting variable, for P.
4. Method
Inspired by PAN, the PAN+SR framework utilizes a one- step nonparametric variable selection strategy to pre-screen high-dimensional dataset (y, X) and parse the reduced dataset (y, Xg) to SR methods for subsequent expression synthesis and selection. Unlike traditional variable selection literature, where the primary focus is controlling the false discovery rates, the PAN criterion calls for minimizing the false negative rate (FNR) while controlling the false positive rate (FPR) is secondary. In other words, the selected set of features S should be a superset of So and as small as possible. When S fails to be the superset of So (i.e., there is at least one FN), the reduced search space F(Xg: g: O) no longer contains f.).). rendering any subsequent discovery based on Xg to be false.
j = ,.....p. The variable inclusion proportion (VIP) qj of æt is then estimated as the posterior mean of these pro- portions across all posterior samples (Chipman et al., 2010), Intuitively, 91.............................. Up encode the relative importance of each feature, where large VIP Tj suggests x, being an impor- tant driver of the response y. However, deciding on how large a VIP value must be to indicate relevance remains a challenge. For instance, BART-G.SE addresses this by using a permutation test on Q1 ........ Op to identify significant features, whereby controlling the family-wise error rate.
Here, we propose alternative approach that utilizes the rankings of VIPs instead of their raw values. Specifically, let r3 denote the ranking of the VIP qj. Relevant features Xs, are expected to occupy top-ranking positions, namely [l.....Po} due to their strong associations with y. In con- trast, irrelevant features Xs. S1 [p] \ So are expected to appear in lower-ranking positions, namely {Po + I,....p). since they are only selected sporadically by chance (Chip- man et al., 2010; Bleich et al., 2014). Consequently, a natural decision rule is to select feature æ, if ry falls within Po,....Pp}.
Nonparametric model-free variable selection has been ex- tensively studied in the literature. Lafferty and Wasserman (2008) propose the RODEO method for nonparametric vari- able selection through regularization of the derivative expec- tation operator. Candès et al. (2018) propose a model-free knockoff procedure controlling FDR with no assumptions on the conditional distribution of the response. Fan et al. (2011) propose sure independence screening method for B-spline additive model. In the Bayesian literature, Bleich et al. (2014) design permutation tests for variable inclusion proportion of Bayesian Additive Regression Tree (BART); Liu et al. (2021) deploy spike-and-sla priors directly on the nodes of Bayesian forests.
However, this decision rule is impractical in real-world ap- plications since the sparsity Po is unknown. To address this limitation, we propose method that leverages multiple in- dependent runs of BART to estimate the feature rankings more robustly. Let Tj.k denote the VIP ranking of æ; in the kth run. Assume that the rankings of æy are randomly distributed over the K independent runs (see Appendix D. for empirical justification):
Despite this diverse array of methods, few meet the unique proposition of the PAN criterion. Among the few recent methods investigated in Ye et al. (2024), they found BART- G.SE (Bleich et al., 2014), a BART-based permutation variable selection method, to be particularly suitable for PAN. However, our comprehensive simulation study in Ap- pendix D.2 reveals that BART-G.SE, along with three other methods, exhibit insufficient TPR, particularly under noisy or low-sample-size conditions. This deficiency renders these methods unsuitable for the PAN+SR framework.
iid ~
if E So Unif((1o + +1, (...,,,)), if $ So
Tj. T.
Then the average ranking rj. = =K Tj,/ K of æ, across K independent runs forms two distinct clusters, Co for Xs. and C₁ for Xs,- Specifically, r,. for Xs, are expected to cluster in Co with mean (1 + po)/2. while those for Xs, tend to cluster in C1 with mean (Po 1 + p)/2. Although both cluster means are unknown due to the unknown spar- sity Po, their separation can be identified using clustering techniques.
In this paper, we introduce a novel BART-based variable selection method and demonstrate its PAN criterion consis- tency through an extensive simulation study in Section 6.2. The key idea behind BART is to model the regression func-
To illustrate, consider the extended Feynman 1-38-12 12 dataset (defined in Section 5.2) with p = 204 features, of which
3

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
Po 4 are relevant. Without loss of generality, we as- sume that the relevant features Xs, are 11; 22, 33; az. 2, 4) and Si = { 20.).. When K = 20 independent BART models are trained on the dataset, the rankings TI,A: "2,,:733,A T frequently fall within 2,, 3, 43 across all = I,...,,3 20 runs. This is 2,
signal-to-noise ratio. Datasets were split 75%/25% in train- ing and testing. For black-box problems, hyperparameters were either set to the optimal values published by SRBench or to values recommended by the original authors of the respective methods. The best hyperparameter settings in black-box regression problems were used in ground-truth problems. Instructions for reproducing the experiment is available in Appendix A, and detailed experimental settings are described in Appendix C.
i.e., So
because the relevant features are frequently selected for tree splits due their strong associations with the response vari- able y. leading to high VIPs and consistently top rankings. In contrast, irrelevant features 25, ......, 2004 are included sporadically in BART, with is.A: ..........2... 2204.6 distributed ran- domly across {5...., 2043. As evident in Figure 5 in Ap- pendix D.I, the average VIP rankings Fy. of the relevant features form a low-mean cluster Co with cluster mean of (1 + po)/2
5.1. Symbolic Regression Methods
Here we summarize the SR methods evaluated in this pa- per. A long strand of SR methods is based on genetie pro- gramming (GP), a technique for evolving executable data structures, such as expression trees. The most vanilla ver- sion we test is gplearn (Stephens, 2020), which performs random expression proposal and iterates through the steps of tournament selection, mutation, and crossover. Advanced GP-based methods utilize different evolutionary strategies and optimization objectives, ranging from Pareto optimiza- tion for efficient trade-offs between accuracy and model complexity to program semantics optimization for increas- ing coherence in expression. Here we test an array of ad- vanced GP-based SR algorithms, including Age-Fitness Pareto optimization (AFP) (Schmidt & Lipson, 2010), AFP with co-evoloved fitness estimate (AFP FE) (Schmidt & Lip- son, 2010), Epigenetic Hill Climber (EHC) (La Cava et al., selection (EPLEX) (La Cava et al., 2019a), Feature Engineering Automation Tool (FEAT) (L. Cava et al., 2019b), Fast Function Extraction (FFX) (McConaghy, 2011). GP version of Gene-pool Optimal Mixing Evolu- tionary Algorithm (GP-GOMEA) (Virgolin et al., 2021), 1), Interaction-Transse Transformatinn Evolutionary Algorithm (ITEA) (de Franca & Aldeia, 2021), Multiple Regression Genetic Programming (MRGP) (Arnaldo et al., 2014), Operon (Burlacu et al., 2020), PySR (Cranmer, 2023), and Semantic Back pack-propagation Genetic Programming (SBP-GP) (Vir- golin et al., 2019).
2.5, while those of the irrelevant fea- tures form a high-mean cluster C₁. concentrating around (Po + +1+ /2 = 104.5.
However, the sparse regression setting naturally leads to a class imbalance problem as |Col Po is much smaller than |C1| = - Pu- To this end, we propose to apply agglomera- tive hierarchical clustering (AHC) with Fuclidean distance and average linkage to (f1....., Fp.)) and cut the dendro- gram to form two clusters: Cg and C.. Then, features in Co are retained, while those in C₁ are discarded. Notably, the proposed data-driven selection criterion does not require any knowledge about the sparsity level po or a tunable selection threshold. An ablation study evaluating the effect of differ- ent clustering algorithms on selection accuracy is available in Appendix D.3. We herein refer to this variable selection method for SR pre-screening as PAN; see Appendix C.2 for implementation details.
5. Experiment Design
Using an open-source benchmarking platform, SRBench, we evaluate the PAN+SR framework on two separate tasks. First, we assess its ability to make accurate predictions on "black-box" regression problems in which the underlying regression function remains unknown. Second, we test PAN+SR's ability to find the correct generatinn func- tion In on synthetic datasets with known data-generating functions originating from Feynman Lectures on Physics (Feynman et al., 2010; Udrescu & Tegmark, 2020).
Additional methods include Bayesian Symbolic Regression (BSR) (Jin et al., 2020), which places a prior on the ex- pression tree; Deep Symbolic Regression (DSR) (Petersen et al., 2021), Unified Deep Symbolic Regression (uDSR) (Landajuela et al., 2022), and Dynamic Symbolic Network (DySymNet) (Li et al., 2024) utilize recurrent neural net- works to propose symbolic expressions; Transformer-based Planning for Symbolic Regression (TPSR) (Shojaee et al., 2023) leverages pretrained transformer models; AlFeynman 2.0 (Udrescu al., 2020) which uses a divide-and-conquer technique to recursively decomposing complex problems into lower-dimensional sub-problems.
The experiment settings are summarized in Table 1. All experiments were run on a heterogeneous cluster. Each al- gorithm was trained on each dataset in 10 repeated trials with a different random state to control both the train/test split and the seed of the algorithm. Each run was performed until 24-hour time limit was reached or up to 500,000 expression evaluations for black-box problems or 1,000,000 for ground-truth problems. For ground-trum -trth problems, we chose a few representative algorithms in the black-box prob- lems and investigated additional settings of sample size and
4

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
Table I: Settings used in the experiments.
SETTING
BLACK-BOX PROBLEMS
GROUND-TRUTH PROBLEMS
# OF DATASETS
35
100
# OF ALGORITIMS
19
19
OF TRIALS PER DATASET
10
10
TRAIN/TEST SPLIT
.75/.25
.75/.25
TERMINATION CRITERIA
500K EVALUATIONS OR 24 HOURS
IM EVALUATIONS OR 24 HOURS
SAMPLE SIZE
ALL
500, 1000, 1500. 2000
SIGNAL- TO-NOSE RATIO
NONE
0.5, 2, 5, 10, 15, 20, NONE
TOTAL COMPARISONS
12250
142000
COMPUTATION COST
34K CORE HOURS
104K CORE HOURS
MEMORY ALLOCATION
16 GB
16 GB
5.2. Datasets
X = (ar). jam........................ 23,.n E < is the irreverent feature matrix induced by the jth relevant feature for j = 1, ........... Po- totaling p = p1(++)) features. In Section 6.2, we fix = 50
We curated a database of high- high-dimensional regression prob-
lems for testing the capability of PAN+SR. We selected 35 black-box regression problems available in PMLB v1.0 (Ro- mano al., 2021) using the following criteria: n < 200 and p 2 10 or n 2 200 and p 2 20. These problems were used in SRBench and overlap with various open-source reposito- ries, including OpenMI. (Vanschoren et al., 2014) and the UCI Machine Learning Repository (Kelly et al., 2013).
so the total number of features is p = 51pb- Additional dataset information and sampling process available in Appendix B.
Besides the 3,200 distinct simulation settings described in Table 1 (100 datasets, 8 SNRs, and 4 sample sizes), we include additional simulation settings in Appendix D.4 to further assess PAN+SR's 's behavior under alternative feature structures. These include (1) additive noise in features, (2) duplicated features, and (3) correlated features.
We also curated 100 high-dimensional ground ttutt regres- sion problems by modifying the Feynman Symbolic Re- gression Database (Udrescu & Tegmark, 2020) include irrelevant features and white noise. For each equation f-() in the Feynman Lectures on Physics, we generated the rele- vant features Xs, following Udrescu and Tegmark (2020):
5.3. Metrics
Predictive Accuracy We assessed predictive accuracy using the coefficient of determination, defined
(I (Figi In.)) Uni((a;, by), for 1 S j S po, (4)
Ra 11 i= (3: 0.)
where po = |So| is the number of relevant features, is the sample size, and aj and bj the lower and upper bounds for feature æ described in Udrescu and Tegmark (2020). To study the effect of noise on PAN+SR, we tuned the signal- to-noise ratio (SNR) by adding a Gaussian error term when generating the response variable:
M1 g)
Model Complexity In line with SRBench, define model complexity as the total number of mathematical operators, features, and constants in the model. To avoid redundancy, symbolic models are first simplified using SymPy (Meurer et al., 2017), a Python library for symbolic mathematics.
y: = fo Ii.po) ) + te, for 1 s S n, (5)
Solution Criteria For ground- -ruth regression problems, we follow SRBench's definition of symbolic solution. A model f(X) is considered a solution to the SR problem of = fo(X) e if x(x) does not reduce to a constant and (1) f f fo = for some a E R or (2) f/Jo b for some + 0. That is, the predicted model f only differs from the true model fo by either additive or multiplicative constant.
where E; sid NO. oF SNR When of 0
or SNR = xo, (4) and (5) generate the original Feynman Symbolic Regression Database.
In addition to the relevant features Xs. = (x) ) we included an array of irrelevant features XIx, represent- ing the era of big data where all reasonable features are in- cluded in the dataset. Specifically, for each relevant feature
Zz. E Sp. we generate (x)
iid (z)),nr jim ... ))m) ~ Unif(a,, by)) representing s copies of independent and irrelevant features coming from the same distribution as xy. Then, the final fea- ture matrix is X [Xs, X!
While predictive accuracy can be influenced by the simula- tion design, the symbolic solution criterion offers a more reliable metric for assessing whether an SR method can uncover the true -ana-gnnng process. However, since
P E Raxp, where in
im
5

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
Opronn
PySR
SBP-GP
FEAT
EPLEX
UDSR
TPSR
GP-GOMEA
AFP FE
DySymNet
EHC
AFP
FFX
DSR
MRGP
oplear
ITEA
AIF eynman
BSR
0.00
0.25
0.50
0.75
1.00
2000
4000
6000
20000
40000
60000
80000
R² Test
Model Size
Training Time (s)
PAN+SR
SR
Figure 1: Results on the black-box regression problems. Points indicate the mean test set performance and bars represent the 95% confidence intervals. Training time for PAN+SR includes the runtime of PAN, which averages only 74. seconds.
SymPy's simplification process is not always optimal, it is possible that some symbolic solutions not identified in the process.
In addition accuracy gain, PAN+SR significantly reduces training times for several SR algorithms, including SBP-GP, uDSR, AFP_FE, AlFeynman, and BSR. Notably, AlFeyn- man, the 2nd slowest running SR algorithm. achieves 5-fold speedup (from 71250 seconds to 1 13997 seconds), while uDSR benefits from nearly 3-fold speedup (from 7628 seconds 2612 seconds) with PAN pre-screening. The computational overhead introduced by PAN is mini- mal, averaging only 74.14 seconds single core. As PAN relies on independent MCMC chains, this overhead can be further reduced through parallel processing, making PAN+SR both efficient and scalable.
Feature Usage Accuracy The irrelevant features present a unique challenge for SR methods to identify the correct Kata-generating model fo. When the predictive model f includes irrelevant features (FPs), cannot be considered excludes some
a symbolic solution to fo- Conversely, relevant features (FNs), it also fails meet the symbolic so- lution criteria. Although neither FPR nor FNR corresponds directly to symbolic solution rate, they provide insights into why does not qualify as a symbolic solution.
6.2. Ground-truth Datasets
6. Results
Figure 2 summarizes performance on the ground-truth re- gression problems with n = 1000, SNR = 0o, and x = 50. Methods are sorted by their standalone R²² on the test set. PAN+SR consistently improves both I22 and solution rate across all 19 SR methods. Due the high dimensionality of the ground-truth problems, the standalone AlFeynman encountered out-of-memory errors and failed to complete any of the 1000 runs. However, PAN significantly improves AlFeynman's performance lifting it from last place to 2nd overall in symbolic solution rate. Furthermore, PAN consis- tently outperforms all other nonparametric variable selection methods tested, achieving the highest TPR among four other methods and delivering the best R when paired with SR, detailed in Appendix D.2. This underscores the effec- tiveness and necessity for nonparametric pre-screening in
6.1. Blackbox Datasets
Figure 1 shows that PAN+SR consistently improves test set R2 across 18 out of 19 SR algorithms, with the largest gains observed in powrr-peroorming methods such as BSR, AlFeynman, and ITEA. For top-perfie rop-performing SR algorithms, the improvements are more modest due the natural upper limit of R2, but the uplift remains significant. For instance, PAN boosted uDSR from 14th to 5th place in the overall ranking and to 2nd among the standalone SR methods. Fur- thermore, these R² improvements are not accompanied by increased model complexity. some cases, PAN+SR even reduces model complexity, enhancing both parsimony and interpretability.
6

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
TPSR
PySR
MRGP
GP-GOMEA
EPLEX
SBP-GP
Operon
UDSR
AFP FE
ITEA
FEAT
FFX
AFP
yySymNet
EHC
gplearn
DSR
BSR
AlF eynmon
0.00
0.25
0.50
0.75
1.00
20
40
60
500
1000
1500
20000
40000
60000
Test
Solution Rate (%)
Model Size
Training Time (s)
PAN+SR
SR
Figure 2: Results on the ground truth regression problems with n = 1000, SNR oo, and 50. Points indicate the mean test set performance and bars represent the 95% confidence intervals. Training time for PAN+SR includes the runtime of PAN, which averages 325 seconds. AlFeynman fails to complete any run the standalone setting.
high-dimensional SR problems.
be concern for typical variable selection applications, it is crucial to emphasize that PAN's primary objective is to scale up SR methods by reliably identifying superset of the relevant features Sg. this context, minimizing FNs during pre-screening more critical than avoiding FPs.
Similar to our findings in the black-box regression problems, this performance gain is not driven by increased model size, and PAN's average computational overhead of 325 seconds remains insignificant to many SR methods. Remarkably, uDSR benefited from nearly a 6-fold speedup with PAN (from 9573 seconds 1,966 seconds) while almost dou- bling its solution rate (from 36.6% to 71.8%), making the best performer in solution rate. Additionally, PAN elevated several mid-tier performers such as Operon, AFP_FE, AFP, and EHC, enabling them to surpass the 4th place method, GP-GOMEA, in standalone SR solution rate ranking.
Figure 3b illustrates that PAN achieves near 0% FNR across most simulation settings, highlighting its ability identify superset of the true feature set S.. This is crucial to ensure that the pre-screened dataset (y. X used for sub- sequent SR modeling is comprehensive enough to generate the correct expression fu- However, in the most extreme case, where n = 500 and SNR = 0.5, PAN's FNR rises to over 5%, and caution advised when relying on PAN in such cases. On the other hand, the standalone Operon often fails include all relevant features in its final mod- els across all and SNR settings, while PAN consistently lower Operon's 's FNR, enhancing its chance to identify the true function fo- Even with PAN, Operon fails to achieve the best-case FNR set by PAN, particularly under noisy con- ditions. This elevated FNR negatively impacts Operon's solution rate. For example, changing SNR from co 10, PAN+Operon's 's average solution rate drops from 27.4% to 0%, and Operon's solution rate falls from 18.1% to 0%. As La Cava et al. (2021) noted, this limitation persists even when Operon is provided with only the relevant features Xs, and under favorable conditions (n = 100, 000 and SNR = 100), indicating that the issue lies beyond PAN pre- screening. Other performance metrics of this sensitivity analysis are available in Appendix D.5.
Beyond the specific simulation setting of n 1000 and SNR oo, also investigated the sensitivity of PAN+SR across range of sample sizes and SNR. In particular, we evaluated PAN+SR with all combinations of sample size E {500, 1000, 1500, 2000} and SNR E (0.5, 1, 5, 10, 15, 20, oo). Given the extreme computa- tional burden, we select Operon, the best-performing al- gorithm in black-box regression problems, to be the SR module for the sensitivity analysis.
Figure 3a demonstrates that both Operon and PAN+Operon maintain consistently lower FPR across all settings of n and SNR, with negligible differences between them. This low FPR reflects the rare inclusion of irrelevant features in the final symbolic models. noisy settings, we notice significant increase PAN's FPR, from 0% SNR co to 30% at SNR 0.5. While this noise sensitivity could

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
= 500
n 1000
0 =50
n= -0000
30
20
20
10
FPR
1500
2000
n=1500
2000
23
30
20
15
10
10
20
10
5
0.5
20
13
19
potspless
10
SNR
SNR
PAN-5R
5R
PAN
PANHS 58 - PAN
(a) False positive rate (FPR).
(b) False negative rate (FNR).
Figure 3: FPR and FNR of Operon, PAN+Operon, and PAN on the ground truth datasets. PAN refers to the proposed selection method Section 4. Points indicate the mean performance and bars represent the 95% confidence intervals.
DSR, the weakest SR algorithm in Figure 4, increasing its solution rate from 8.2% 14.9% at SNR 10 and from 8.9% to 25.8% at SNR co. These findings highlight the fundamental challenges noise introduces to SR algorithms. To date, SR algorithms have been predominantly developed for noiseless or high-SNR settings, even for "small p" prob- lems. We expect that iterative application of the proposed variable selection method, similar to Ye et al. (2024), along with careful consideration of the challenges in extreme-scale SR, could improve performance in low-SNR settings. This will be explored in future work.
GP-GOMEA
peran
40
4
6
A
4
FEAT
4
AFP
EHC
2
20
DSR
0.5
0.6
0.7
0.8
0.9
1.0
20
40
60
7. Discussion
R²
Test
Solution
PAR noselles
A 5R, noiseless
RAN-5R.SNR-D A SN8-10
In this paper, we introduce PAN+SR, a novel framework designed address the scalability challenges faced by SR methods when applied to high-dimensional datasets. The growing prevalence of big data necessitates tools capa- ble of efficiently handling such complexity, and PAN+SR addresses this need by integrating a nonparametric pre- screening mechanism with SR. This integration enables the framework to focus the model search on relevant subset of features, reducing computational burden and improving accuracy.
Figure 4: Results of selected methods the ground truth problems with 1000, SNR foo, 10), and 50. Points indicate the mean test set performance and bars rep- resent the 95% confidence intervals.
Beyond Operon, we also evaluated several top- pe-perfformng SR methods on the ground-truth problems using = 1000 and SNR E foo, 10). As shown in Figure 4, PAN+SR consistently improves SR methods across all SNR levels, though all SR and their PoN-boosted variants become less accurate at SNR 10, indicating the challenge when noise is present. In particular, GP-GOMEL performs similarly to Operon, with its solution rate dropping to 0% at SNR = 10 for both the standalone and PAN-boosted variants. The best- performing SR algorithm, uDSR, also exhibits vulnerability to noise, with its PAN-boosted solution rate falling from 71.8% to 7.4%. Surprisingly, PAN significantly benefits
The core innovation of PAN+SR lies in its nonparametric variable selection method, which filters the imput dataset to reduce dimensionality before applying SR. A key chal- lenge in this process is minimizing the risk of false nega- tives (FNs), where relevant features are mistakenly excluded. Such omissions can critically impair SR methods, as the success of SR depends having access to the true feature set. To address this issue, we developed variable selec- tion method designed to ensure that the selected features
8

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
form superset of the true feature set, effectively minimiz- ing the FNR. Our approach leverages the characteristies of VIP rankings derived BART, providing a Juning-free, data-driven variable selection criterion capable of retaining relevant features while excluding irrelevant ones. By pre- serving a comprehensive set of candidate features, PAN+SR maximizes the likelihood of identifying the true underlying model.
References
Arnaldo, L., Krawiec, K., and O'Reilly, U.-M. Multiple regression genetic programming. In Proceedings of the 2014 Anmual Conference on Genetic and Evolutionary Computation, GECCO '14, pp. 879-886, New York, NY, USA, 2014. Association for Computing Machinery.
Bleich, J., Kapelner, A., George, E. I., and Jensen, S. T. Variable selection for BART: an application to gene regu- lation. Annals of Applied Statistics, (3):1750-1781, 09 2014.
We evaluated PAN+SR across diverse set of datasets, in- cluding 35 high-dimensional real-world datasets from the PMLB database and 100 modified simulated datasets based on the Feynman Symbolic Regression Database. The results were highly promising: PAN+SR improved the performance of 18 out of 19 SR methods on real datasets and all 19 meth- ods on simulated datasets when noise is absent. These findings underscore the framework's ' potential to enhance the robustness and scalability of SR methods across diverse datasets.
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, GECCO '20, pp. 1562-1570, New York, NY, USA, 2020. Association for Computing Machinery.
Candès, E.,, Fan, Y., Janson, L., and Lv. J. Panning for Gold: 'Model-X' Knockoffs for High Dimensional Controlled Variable Selection. Journal of the Royal Statistical Soci- ety Series B: Statistical Methodology, 80(3):551-577 01 2018.
In addition, we explored the sensitivity of PAN+SR to vary- ing sample sizes and SNR. Our analysis demonstrated that the performance gains achieved by PAN+SR are consis- tent across different sample sizes and remain robust in the presence of noise. Like our extended Feynman database, SDSR (Matsubara et al., 2024) augments the original Feyn- man database with irrelevant features, bringing the synthetic benchmarks closer to real-world scientific process. However, SDSR adds only 1-3 irrelevant variables, while our setup introduces 100-450 irrelevant variables, posing a substan- tially more challenging test for both variable selection and symbolic regression. Nonetheless, SDSR rectifies several physical inconsistencies present in the original Feynman benchmark, such as more realistic treatment of constants and integer-valued variables, and a more careful specifica- tion of sampling ranges. Our investigation extends beyond ground-truth datasets by incorporating black-box datasets, thereby mitigating, to some extent, the limitations inherent in purely simulated data. Still, view SRSD as valuable and complementary benchmark and plan to incorporate its refinements in future evaluations. In summary, PAN+SR provides a significant step forward in enabling SR methods to handle the complexities of modern datasets, offering im- proved performance and scalability across a wide range of applications.
Chipman, H. A., George, E. I., and McCulloch, R. E. BART:
Bayesian additive regression trees. Annals of Applied Statisties, 4(1):266-298 03 2010.
Cranmer, M. Interpretable Machine Learning for Science with PySR and SymbolicRegression.jl. arXiv:2305.011 008522 2023.
de Franca, F. O. and Aldeia, G. S. I. Interac- tion-transformation evolutionary algorithm for symbolic regression. Evolutionary Computation, 29(3):367-390, 09 2021.
Dernrr, E., Kubalík. J., Ancona, N., and Babuška, R. Con- structing parsimonious analytic models for dynamic sys- tems via symbolic regression. Applied Soft Computing, 94:106432 2020.
Dick, G. Genetic programming, standardisation, and stochastic gradient descent revisited: initial findings on srbench. In Proceedings of the Genetic and Evolutionary Computation Conference Companion, GECCO '22, pp. 2265-2273 New York, NY, USA, 2022. Association for Computing Machinery.
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 of which we feel must be specifically highlighted here.
Fan, J., Feng, Y., and and, R. S. Nonparametric indepen- dence screening in sparse ultra-high-dimensionaal additive models. Journal of the American Statistical Association, (06(494):544-551) 25,,4-555,1,111115
Feynman, R. P., Leighton, R. B., and Sands, M. The Feyn- man Lectures on Physics. Basic Books, New York, NY, 2010.
9

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
Friedman, J. H. Multivariate Adaptive Regression Splines. The Annals of Statistics, 19(1):1-67, 1991.
La Cava, W., Orzechowski, P., Burlacu, B., de Franca, F, Virgolin, M., Jin, Y., Kommenda, M., and Moore, Con- temporary symbolic regression methods and their relative performance. In Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks, volume 1, 2021.
Hernandez, A., Balasubramanian, A., Yuan, F., Mason, S. A. M., and Mueller, T. Fast, accurate, and transferable mody-body interatomic potentials by symbolic regression. upj Computational Materials, 5(1): 112, November 2019.
Lafferty, J. and Wasserman, L. Rodeo: Sparse, greedy nonparametric regression. The Annals of Statistics, 36(1): 28-63. 2008.
Jin, Y., Fu, W. Kang, J., Guo, J., and Guo, J. Bayesian Symbolic Regression. arXiv: arXiv:1910.08892 2020.
Kamienny, P.-a., d' Ascoli, S., Lample, G., and Charton, F. End-to-end symbolic regression with transformers. In Koyejo, S., Mohamed, S., Agarwal, A., Belgrave, D., Cho, K., and Oh, A. (eds.), Advances in Neural Informa- fion Processing Systems, volume 35, pp. 10269-10281 10281. Curran Associates, Inc., 2022.
Landajuela, M., Lee, S., Yang, J., Glatt, R., Santiago, C. P., Aravena, I., Mundhenk, T, Mulcahy, G., and Pe- tersen, B. K. A unified framework for deep symbolic regression. In Koyejo, S., Mohamed, S., Agarwal, A., Belgrave, D., Cho, K., and Oh, A. (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 33985-33995 Curran Associates, Inc., 2022.
Kamienny, P.-A., Lample, G., Lamprier, S., and Virgolin, M. Deep generative symbolic regression with Monte-Carlo- tree-search. In Proceedings of the 40th International Conference on Machine Learning, ICML'23, pp. 15655- 15668. JoLR.org, 2023.
Lemos, P., Jeffrey, N., Cranmer, M., Ho, S., and Battaglia, P. Rediscovering orbital mechanies with machine learning. Machine Learning: Science and Technology. 4(4):045002, October 2023.
Kelly, M., Longjohn, R., and Nottingham, K. The UCI Machine Learning Repository, 2013.
Li, W., Li, W., Yu, L., Wu, M., Sun, L., Liu, J., Li, Y., Wei, S. Yusong, D., and Hao, M. A neural-guided dynamic symbolic network for exploring mathematical expressions from data. In Salakhutdinov, R., Kolter, Z., Heller, K., Weller, A., Oliver, N., Scarlett, J.. and Berkenkamp, F. (eds.), Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 28222-28242. PMLR, 21-27 Jul 2024.
Keren, L. S., Liberzon, A., and Lazebnik, T. A computa- tional framework for physics-informed symbolic regres- sion with straightforward integration of domain knowl- edge. Scientific Reports, 13(1):1249 January 2023.
Kronberger, G., Kommenda, M., Promberger, A., and Nickel, F. Predicting friction system performance with symbolic regression and genetic programming with factor variables. In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO '18, pp. 1278-1285, New York, NY, USA, 2018. Association for Computing Machinery.
Liu, C.-Y., Zhang, S., Martinez, D., Li, M., and Senftle, T. P. Using statistical learning to predict interactions between single metal atoms and modified MgO (100) supports. npj Computational Materials, 6 6(2, 2020.
Liu, C.-Y., Ye, S., Li, M., and Senftle, T. P. A rapid feature selection method for catalyst design: Iterative Bayesian additive regression trees (iBART). The Journal of Chemi- cal Physics, 156(16), 2022.
La Cava, W., Spector, L.,, Danai, K., and Lackner, M. Evoly- ing differential equations with developmental linear ge- netic programming and epigenetic hill climbing. In Pro- ceedings of the Companion Publication of the 2014 An- mual Conference on Genetic and Evolutionary Computa- tion, GECCO Comp 14, pp. 141-142, New York, NY, USA, 2014. Association for Computing Machinery.
Liu, Y., Ročková, V., and Wang. Y. Variable selection with ABC Bayesian forests. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 83(3):453- 481, 04 2021.
La Cava, W., Helmuth, T., Spector, L., and Moore, J. H. A Probabilistic and Multi-Objection Analysis of Lexi- case Selection and e-Lexicase Selection. Evolutionaar Computation, 27(3) September 2019a.
Makke, N. and Chawla, S. Interpretable scientific discovery with symbolic regression: a review. Artificial Intelligence Review, 57(1): January 2024.
La Cava, W., Singh, T. R., Taggart, J., Suri, S., and Moore, J. Learning concise representations for regression by evolving networks of trees. In International Conference on Learning Representations, 2019b.
Märtens, M. and Izzo, D. Symbolic regression for space applications: Differentiable cartesian genetic program- ming powered by objti-objective memetic algorithms. arXiv:2206.062/3 2022.
10

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
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.
Udreseu, S.-M., Tan, A., Feng, J., Neto, O., Wu, T., and Tegmark, M. Al Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), Advances in Neural Information Processing Sys- tems, volume 33, pp. 4860-4871. Curran Associates, Inc., 2020.
McConaghy, T. FFX: Fast, Scalable, Deterministic Sym- bolic Regression Technology, pp. 235-260. Springer New York, New York, NY, 2011.
Vanschoren, J., van Rijn, J. N., Bischl, B., and Torgo, L. Openml: networked science in machine learning. SIGKDD Explor. Newsl., 15(2):99-60, June 2014.
Meurer, A., Smith, C. P. Paprocki, M., Čertík, O., Kirpichev,
S. B., Rocklin, M., Kumar, A., Ivanov, S., Moore, K., Singh, S., Rathnayake, T., Vig. S., Granger, B. E.,, Muller, R. P., Bonazzi, F., Gupta, H., Vats, S., Johansson, F,, Pedregosa, F., Curry, M. J., Terrel, A. R., Roučka, v., Saboo, A., Fernando, I., Kulal, S., Cimrman, R., and Scopatz, A. SymPy: symbolic computing in Python. PeerJ Computer Science, 3:3, January 2017.
Verstyuk, S. and Douglas, M. R. Machine learning the gravity equation for international trade. Available at SSRN 4053795, 2022.
Virgolin, M. and Pissis, S. P. Symbolic regression is NP- hard. Transactions Machine Learning Research, 2022.
Petersen, B. K., Larma, M. L.,, Mundhenk, T. N., Santi- ago, C. P., Kim, S. K., and Kim, T. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, 2021.
Virgolin, M., Alderliesten, T., and Bosman, P. A. N. Linear scaling with and within semantic backpropagation-batt genetic programming for symbolic regression. In Pro- ceedings of the Genetic and Evolutionary Computation Conference, GECCO '1,, pp. 1084-1092, New York, NY, USA, 2019. Association for Computing Machinery.
Romano, J. D., Le, T. T., La Cava, W., Gregg, J. T., Gold- berg, D. Chakraborty, P., Ray, N. L., Himmelstein, D., Fu, W., and Moore, H. PMLB v1.0: an open-source dataset collection for benchmarking machine learning methods. Bioinformatics, 38(3):878-880, October 2021.
Virgolin, M., Wang. Z., Alderliesten, T., and Bosman, P. A. N. Machine learning for the prediction of pseudore- alistic pediatric abdominal phantoms for radiation dose reconstruction. Journal of Medical Imaging, 7(4):046501, 2020.
Schmidt, M. D. and Lipson, H. Age-fitness pareto optimiza- tion. In Proceedings of the 12th Annual Conference Genetic and Evolutionary Computation, GECCO 10, pp- 543-544, New York, NY, USA, 2010. Association for Computing Machinery.
Virgolin, M., Alderliesten, T., Witteveen, C., and Bosman, P. A. N. Improving Model- Genetic Programming for Symbolic Regression of Small Expressions. Evolutionary Computation, 29(2):211 1-27,, 06 2021.
Shojaee, P., Meidami, K., Barati Farimani, A., and Reddy, C. Transformer-basee planning for symbolic regression. In Oh, A., Naumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S. (eds.), Advances in Neural Informa- tion Processing Systems, volume 36, pp. 45907-45919. Curran Associates, Inc., 2023.
Wu, T. and Tegmark, M. Toward an artificial intelligence physicist for unsupervised learning. Phys. Rev. E, 100: 033311, September 2019.
Ye, S., Senftle, T. P., and Li, M. Operator-Induced Struc- tural Variable Selection for Identifying Materials Genes. Journal of the American Statistical Association. 119(545): 81-94, 2024.
Stephens, T. gplearn: Genetic Programming in Python. https:: //git wwi.ihhub.cc.m wtwpab/cc//tteeestteeeeeeee gplearn, 2020.
Strogatz, H. Nonlinear Dynamics and Chaos: With Appli- cations to Physics, Biology, Chemistry, and Engineering. CRC Press, 2015.
Tenachi, W., Ibata, R., and Diakogiannis, F. I. Deep sym- bolic regression for physics guided by units constraints: Toward the automated discovery of physical laws. The Astrophysical Journal, 959(2):99, December 2023.
Udreseu, S.-M. and Tegmark, M. AI Feynman: A physics- inspired method for symbolic regression. Science Ad- vances, 6(1 (633) 2020.

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
A. Reproducing the Experiment
The experiment made use of an existing symbolic regression (SR) benchmarking platform, SRBench (La Cava et al., 2021). and changes were made facilitate other functionalities, including signal-to-noies ratio (SNR) tuning, feature pre- sereeennnnng and variable usage accuracy calculation. The README file in our GitHub repository details the complete of commands for reproducing the experiment. Here, provide short summary of the experiment process. Experiments are launched from the experiments folder via the script ana yyze. Py. After installing and configuring the conda environment provided by SRBench, the complete black-box experiment on standalone SR methods can be started via the following command:
cwmmmmttsseng/////e... .RR
python analyze.py
-results /path/to/pmlb/
-n_trials 10
limie_limit. 24:00 -tuned -skip_tuning
To enable PAN pre-screening, the users can either specify the path a pre-run variable selection result or the pre- screening in place. The first option is useful when the users need to compare different SR methods the same dataset:
python analyze. py /path/to/pmlb \ -results //results blackbox/ SR _ARI_VIP VIP
-n_trials 10
-time_limit 24:00
-vs_method BARI_VIP \
-vs_rcsult_pathh creuultssbblcckoo/p//lll. http://www...///o...../ooeetetee wwwpa/AAAAAtttttaieaaaa
feather \
-va_idx_label idx_helst
-tuned -kkip_tuning
If no path is given to -vs resuut.pt.p the PAN pre-screening will be run in place. Similarly, the ground-truth experiment for the standalone SR methods Feynman datasets with sample size of 1000 and an SNR of 10 can be run by the following command:
python analyze.py /path/to/feynna
-results weesulss_feymmann
sornal_to_nciss 10
-n 1000
-sym_data
-r_trials 10
-time_limit 24:00
-tuned -skip tuning
Note that sym_datt enables more performance metric calculations only available for ground-truth problems. To run PAN pre-screening only the Feynman datasets with sample size of 1000 and an SNR of 10, the following command:
python analyze py /ptth/tttttttepeeee fatth/tto/eeyman feynman
-script BART_selection \
-ml BART_VIE
-results .. ../results. /
ssignal_te_nat
10
-n 1000
-sym_data \
-n_trials 10 \
-rep 20
10
-ime 1imit 24:00
The -rep 20 argument instructs the program to run K 20 replications BART for estimating the variable ranking Tju
of the jith feature at the kth run. Users can use other variable selection methods by modifying the BART.selection. script.
12

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
B. Additional Dataset Information
PMLB datasets Black-box datasets and their metadata are available from PMLB under an MIT license and described in detail in Romano et al. (2021). In this experiment, we only focus on disensinnn regression datasets available from PMLB. Specifically, we use PMLB regression datasets satisfying the following criteria:
1. n < 200 and p 2 10, or
2. 2 200 and p 2 20.
Furthermore, datasets that have categorical features (number of unique value I 5) or non-continuoss response variable (proportion of unique value < 0.9) are excluded since they are incorrectly classified as regression task (Dick, 2022). Among the datasets meeting these criteria, we found that two datasets, wa..uttoppiii and 207 Laut ce, are identical, and we only kept 1 a5_o pprcee in our analysis. See Dick (2022) for detailed analysis of the dataset duplication and incorrect problem classification issues of PMLB.
Feynman datasets The original Feynman database described in Udrescu and Tegmark (2020) consists of only the rel- evant features Xs, and a large sample size of n = 105, and is available in Feynman Symbolic Regression Database We extended the Feynman Symbolic Regression Database to include irrelevant features X inr E R.x* for each relevant feature Fj. j E Sp- To take advantage of the SRBench platform, we standardized the Feynman equations PMLB format and included metadata detailing the true model and the units of each variable. The extended Feynman datasets feynman_de t taset_code/ for each sample size considered in the main paper, we set s = 50 and n 100, 000 for all Feynman equations with random state control; we refer to this as the full Feynman datasets. In the experiment, the full Feynman datasets are randomly split into a 75% train/test set. If the train set contains more samples than the desired training sample size n, the train and test sets will be further subsampled so that Ximan has exactly n samples and X has exactly [n/3] samples.
(https://spac..min www.:it...eede htt.ooome//ggaaaaaaaeaamaa.mmm
generated using the Python script provided in .p. To avoid the need to generate different datasets
Users can also generate datasets using other data-generating functions So by supplying CSV file with the ex- pression of fo(-) and an additional CSV file describing the desired uniform distribution (i.e., the lower and upper wtaroed///ymmneuuaaiiioom csv
bounds of the distribution) of each variable in fo(-). See unius.nt............i...... csv for more details.
and feynman.dataset
Sampling Process for Extended Feynman Datasets The sampling process for the extended Feynman datasets is described in the main text and is reproduced here for completeness of the data description in this section.
For each equation fo(-) in the Feymman Lectures on Physics, we generated the relevant features Xs, following Udrescu and Tegmark (2020):
(21.f...., In. iid Unif(ay,b)) b,), for 1 S j S Po:
(6)
where Po = |So| is the number of relevant features, n is the sample size, and aj and bj are the lower and upper bounds for feature ay described generated as follow:
Then, the response variable
y
Ti.po) Ei, for 1 s n,
where e1 sid (O, o²) is an additive Gaussian error, of denotes the sample variance of fo(-), and 0² o3 /2/SNR is the error variance tuned to a prescribed signal-to-noise ratio (SNR). When of = (i.e., SNR = 20), (6) and (7) generate the original Feynman Symbolic Regression Database.
For each relevant feature Zj. j 1,...,,oP Po. we generate 50 copies of irrelevant features following the distribution X Unif(a,, by). Then, the final feature matrix is X [Xs,, X1 ...,, X R=Xp, where
of zy: (x) jimm .... s..m) iid
X = (x) (x) jimm... x jiir)) E R=*x is the irreverent feature matrix induced by the jth relevant feature for j = 1,..., Po.
totaling p = Pu(1 + s) number of features.
In Appendix D.4, we consider sampling processes where features are not iid sampled from a uniform distribution.
13

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
C. Additional Experiment Details
C.1. General Experiment Settings
Experiments were run in heterogeneous cluster composed of nodes with Intel(R) (R) CPU E5-2620 v2 @ 2.60GHz, Intel(R) Xeon(R) CPU E5-2650 v4 @ 2.20GHz, Intel(R) Xeon(R) Gold 6126 CPU 2.60GHz, Intel(R) Xeon(R) Gold 6230 CPU @ 2. 2.1 and AMD EPYC 7642 CPU @ 2.3GHz processors. The training of single method on a single dataset for fixed random seed considered a job. Each job was managed by SLURM Workload Manager to receive one CPU core, 16GB of RAM, and time limit of 24 hours. For the ground-truth problems, each final model was given an additional 5 minutes for each of the following steps: 1) cleaning the model for SymPy parsing, 2) simplifying the cleaned model using SymPy, 3) checking the difference solution criterion of the simplified model, 4) checking the ratio solution criterion of the simplified model, and 5) calculating model size (complexity). When the simplification of the cleaned model exceeded the 5-minute wall clock, steps 3-5 were run on the cleaned model instead.
C.2. Implementation Details of the Proposed Variable Selection Method
The proposed method uses the bartMachine R package for its BART implementation. For each dataset, we fit K = 20 independent BART models and record the ranking rj,k of variable a, 's variable inclusion proportion (VIP) the kth run; the hyperparameters for bartMachine are summarized in Table 2. To cluster the VIP rankings into 2 clusters, we use the hel ust function in R to perform agglomeration clustering (unweighted pair group method with arithmetic mean) on the Euclidean dissimilarity matrix of the VIP rankings. Then, 2; is selected if f.. =k=1 J //K belongs to the low-mean
cluster.
Table 2: Hyperparameters bartMachine.
Parameter
Value
# of trees
20
# of burn-in samples
10,000
# of posterior samples
10,000
D. Additional Results
D.1. Visualization of Average VIP Rankings Ty.
Figure shows the average BART VIP rankings for Feynman equation 1-38-12 with n = 1000. At high SNR, there is a clear separation between the low- and high-mean clusters, and the hypothesized cluster means closely match their actual values. As SNR decreases, irrelevant features tend receive higher rankings, slightly shifting the cluster means incurring more false positives (FPs). Despite this deviation, the cluster means remain far apart, ensuring separation between relevant and irrelevant features.
Figure 6 further demonstrates the clustering accuracy of the proposed method. Regardless of the SNR level, all true features are consistently assigned to the low-mean cluster, which is highly desirable in the PAN+SR framework. While decreasing SNR leads to some misclassification of the irrelevant features, the proposed method ensures that no true features are excluded. This robustness in retaining the true features under varying noise levels makes the proposed method well-suited for PAN+SR framework and nigh-dimensional SR tasks.
14

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
SNR naiseless
SNR =00
20
20
20
20
10
10
100
150
50
100
150
SNR =15
SNR =10
30
30
20
20
10
30
100
150
50
100
150
Frequency
SNR 5
SNR -2
30
30
20
70
10
10
100
150
50
100
SNR =1
SNR 0.5
30
30
20
20
10
10
150
Average VIP Ranking
Cluster
bam-menn
high-mean
Figure 5: Average BART VIP rankings Fy. over K runs on Feynman equation 1-38-12 with n = 1000, Po 4, and p 204. Black vertical dashed lines indicate the cluster means. Red solid vertical lines are the hypothesized cluster means: (1 + po)/2 2.5 and (po + 1 + p)/2 104.5.
15

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
SNR = noiseless
SNR =0
130
130
50
0
50
100
150
50
150
250
SNR 15
SNR 10
150
190
Average VIP Ranking
0
0
50
100
150
50
150
150
SNR 55
SNR -2
150
150
100
0
0
50
100
150
50
100
150
SNR =1
SNR 0.5
150
150
0
150
Average VIP Ranking
True Label
True Feature
Inrelventt eoture
Cluster
lon-trenn
hightrean
Figure 6: Hierarchical clustering accuracy on Feynman equation 1-38-12 with n 1000, Po = 4, and p = 204. Red and teal represent the low- and high-mean clusters, respectively. Circles and triangles represent relevant and irrelevant features, respectively.
16

---

Ab Initio Nomparametric Variable Selection for Scalable Symbolic Regression with Large
D.2. Analysis of Different Nonparametric Variable Selection Methods
= 500
n 1000
100
100
80
80
60
70
VIP Rank
Local
1500
2000
G S.
100
100
G MAX
95
90
90
85
nolseless
20
15
10
2
1
0.5
nolseless
20
15
10
2
0.5
SNR
Figure 7: True positive rate (TPR) on the Feynman datasets for n = 500, 1000, 1500, 2000 and SNR 00, 20, 15, 10, 5, 2, 1, 0.5. Points indicate the mean performance, and bars show the 95% confidence interval. VIP Rank is the proposed method for PAN pre-screening. Local, G.SE, G.MAX, and RF are alternative nonparametric variable selection methods.
PAN pre- screening presents a unique challenge to nonparametric variable selection methods where any missed signals (false negative) will eliminate the correct expression fo(-) from the search space. That is, a true positive rate (TPR) near 100% in the pre-screening phase is necessary to ensure successful SR tasks. Figure compares the average TPR of five nonparametric variable selection methods across various configurations of and SNR the Feynman datasets. VIP Rank, the proposed method, is compared with three BART permutation test-based methods (Local, G.SE, and G.MAX) (Bleich et al., 2014) and the Random Forest (RF) variable selection method in PySR (Cranmer, 2023). Of the three BART permutation best-based methods, BART-Local applies the least stringent selection criteria, while BART-G.MAX is the most stringent, with BART-G.SE offering a balance between the two. The RF implementation requires users to specify the number of selected variables k, which we tuned over 21,2..... 20} using 5-fold cross-validation.
VIP Rank consistently achieves the highest TPR, nearing or reaching 100% across all experimental settings. noiseless conditions (SNR- oo), only VIP Rank attains a perfect TPR of 100%. Although there is slight TPR decline for VIP Rank 500 and SNR < 5, it still outperforms the other methods, particularly results reinforce the need for specialized variable selection method for PAN pre-screening In addition to the four methods considered here, we point readers to Ye et al. (2024). where they analyzed three additional nonparametric variable selection methods and showed that none outperform BART-G.SE in terms TPR.
500 and SNR 0.5. These
Figure illustrates the false positive rate (FPR), a crucial metric for evaluating variable selection accuracy. discussed in the main paper, VIP Rank produces higher FPR under low SNR conditions-a tradeoff made to maintain near near-perfect true positive rate (TPR). While this tradeoff may be undesirable for typical variable selection tasks, it is acceptable for PAN pre-screening, where minimizing false negatives (FNs) is the priority. The three BART permutation-based methods and RF consistently maintain low and robust FPRs across all settings of and SNR. However, as Figure 7 shows, this strict control of FPR comes at the cost of worse TPR performance.
To further evaluate the impact of variable selection methods in the PAN+SR framework, we replaced VIP Rank with BART-G.SS and compared their performance using Operon as the SR method. Operon was chosen for this analysis due to
17

---

Ab Initio Nomparametric Variable Selection for Scalable Symbolic Regression with Large
n 500
n 1000
30
30
20
20
VIP Rank Local
FPR
1
=
G SE
G MAX
30
30
20
0
nolseless
20
15
10
5
2
0.5
noiseless
20
15
10
5
2
SNR
Figure 8: False positive rate (FPR) the Feynman datasets for n 500, 1000, 1500, 2000 and SNR 0o, 20, 15, 10, 5, 2, 1, 0.5. Points indicate the mean performance, and bars show the 95% confidence interval. VIP Rank is the proposed method for PAN pre-screening. Local, G.SE, G.MAX, and RF are alternative nonparametric variable selection methods.
2000
its strong Rº performance in both the black-box and ground-truth experiments. Table 3 summarizes the average test set R² on the Feynman dataset. VIP+SR consistently achieves the highest R² across all experimental settings. For instance, at n 500 and SNR= 20, VIP+SR achieves an average I²2 of 0.892, compared to 0.860 for GSE+SR and 0.846 standalone SR. Under high noise conditions, VIP+SR continues demonstrate better robustness than GSE+SR. At n 500 and SNR= 0.5, VIP+SR scores 0.145, slightly outperforming GSE+SR (0.142) and standalone (0.142). This trend is consistent across different different sample sizes
D.3. Effect of Different Clustering Algorithms
The proposed VIP Rank variable selection method be implemented using various off-the-shelf clustering algorithms. However, due to the class imbalance nature of the variable selection problem, not all clustering algorithms are suitable. In this ablation study, we examine the effect of clustering algorithms on TPR and FPR performances of VIP Rank. We elected 10 clustering algorithms available in srikit-leare v1. v1.5.7: 7: agglomerative hierarchical clustering (AHC), k-means++, Gaussian mixture model (GMM), Birch, Mean Shift, Affinity Propagation, Spectral, OPTICS, HDBSCAN, and DBSCAN.
As illustrated Figure 9, the first 5 clustering algorithms (AHC, k-mean++, GMM, Birch, Mean Shift) achieve the highest TPR across all simulation settings with indistinguishable differences. Affinity Propagation also has similar TPR compared with the top 5 algorithms but lacks behind in noisy (e.g., SNR 0.5) and small-n (e.g., 500) settings. The rest of the pack has significantly worse TPR and are thus not suitable in VIP Rank.
Since the top 5 algorithms have indistinguishable TPR, we elect one with the least FPR. As shown Figure 10, AHC has significantly lower FPR than the rest of the top algorithms across most simulation settings. Combine with its near 100% TPR, AHC is capable of identifying more compact feature set that has high probability containing all relevant features.
18

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
n 500
1000
100
100
75
75
50
50
AHC
25
k-means++
GMM
TPR (%)
Birch
1500
n =000
MeanShift
100
100
AffinilyPropggnnn
Spectral
OPTICS
75
75
HDBSCAN
DBSCAN
50
25
25
noiseless
20
15
10
0.5
noiseless
20
15
10
0.5
SNR
Figure 9: True positive rate of various ablations of clustering algorithm.
10000
20
(%)
k-means++
1500
2000
GMM
40
Birch
MeanShift
20
20
nolseless
20
15
10
0.5
noiseless
20
15
10
05
SNR
Figure 10: False positive rate of various ablations clustering algorithm.
19

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
Table 3: Average test set R².
The highest value in each experimental setting is in bold.
noiseless
20
15
10
5
2
0.5
500
VIP+SR
0.974
0.892
0.870
0.837
0.730
0.525
0.335
0.145
GSE+SR
0.948
0.860
0.859
0.818
0.710
0.510
0.327
0.142
SR
0.915
0.846
0.840
0.792
0.702
0.506
0.322
0.142
" = 1000
VIP+SR
0.984
0.919
0.901
0.867
0.774
0.586
0.406
0.229
GSE+SR
0.971
0.914
0.897
0.851
0.774
0.574
0.405
0.229
SR
0.942
0.883
0.867
0.825
0.747
0.580
0.393
0.227
n 1500
VIP+SR
0.990
0.928
0.909
0.874
0.792
0.612
0.433
0.260
GSE+SR
0.961
0.910
0.899
0.866
0.781
0.600
0,428
0.257
SR
0.956
0.895
0.878
0.856
0.761
0.592
0.426
0.255
2000
VIP+SR
0.990
0.935
0.914
0.887
0.805
0.619
0.448
0.277
GSE+SR
0.963
0.918
0.905
0.872
0.787
0.617
0.445
0.272
SR
0.960
0.907
0.892
0.855
0.781
0.611
0.437
0.272
D.4. Effect of Noisy, Duplicated, and Correlated Predictors
In addition to the extensive simulation settings deseribed in Section 5.2, we further evaluate VIP Rank under alternative predictor structures that challenge common modeling assumptions:
Baseline:
21....,2p Zp
Unif 1)
Noisy X: Independent Gaussian noise is added to each predictor with variance equal to 1/5 of the signal variance
Duplicated X: A redundant feature is added:
where and 22 are relevant predictors
Correlated X: ........ Unif(O, 1) with an autocorrelation structure: Psy 0.9
The response variable y is generated according to the Friedman equation (1991):
y 10sin(TZ112)) (00x=-05))² + 10xx + 5xs + E,
E E~ NO0
We fix 1000, p 100, SNR 10, and repeat each scenarios for 100 trials. Table 4 reports the average TPR and FPR. VIP Rank consistently identifies all relevant features across all scenarios, demonstrating strong robustness to noise, redundancy, and correlation among predictors.
Table 4: Average performance in each scenario across 100 trials.
Scenario
TPR
FPR
Baseline
100%
10.58%
Noisy X
100%
26.42%
Duplicated X
100%
11.11%
Correlated X
100%
15.98%
20

---

Ab Initio Nonparametric Variable Selection for Scalable Symbolic Regression with Large
D.5. Additional Performance Metries for Operon PAN+Operon
Figures 11, 12, 13 show additional metries not discussed in the main paper. Although PAN+Operon's is solution rate plummeted from ~27% at SNR 0% at SNR 20 across all n, Figure 11 shows there is still improvement in R on test set across all and SNR, while improving model interpretability evidenced by the uniformly lower model size in Figure 12.
= 500
1000
1.00
0.75
0.50
0.25
1500
n =
1.00
0.75
0.50
0.25
noiseless
20
15
10
5
2
0.5
noiseless
20
15
10
5
2
0.5
SNR
PAN+SR +
SR
Figure 11: R on test set with Operon as the SR module. Points indicate the average R² on test set and bars represent the 95% confidence intervals.
21

---

Ab Initio Nomparametric Variable Selection for Scalable Symbolic Regression with Large
500
1000
80
Model Size
40
n 15000
2000
60
40
noiseless
20
15
10
2
0.5
noiseless
20
15
10
2
1
0.5
SNR
PAN+SR
SR
Figure 12: Model size with Operon the SR module. Points indicate the average model size and bars represent the 95% confidence intervals.
500
= 1000
30
20
10
Solution Rate (%))
n 1500
n =000
20
noiseless
20
15
10
0.5
noiseless
20
15
10
0.5
SNR
PAN+5R
SR
Figure 13: Solution rate with Operon the SR module. Points indicate the average solution rate and bars represent the 95% confidence intervals.
22


Reviews:
[{"review_id": "WIQp2YKToI", "reviewer": "Reviewer_pZx6", "summary": "The authors are interested in the problem of discovering mathematical equations from raw data. One of the largest bottleneck for SR methods is that it's extremely hard to scale the equation search to > 10 variables. This is because each additional variable considered combinatorially increases the search space which makes the search less efficient. The authors hypothesize that, given a large number of variables, we can analyze certain statistical properties to cluster the variables based on their relevance and prune the irrelevant variables, which should considerably increase the search efficiency. Towards this extend, the authors propose PAN-SR, a framework which modifies the BART variable selection method to select a set of relevant variables to run an off-the-shelf SR algorithm on.\n\nThe authors demonstrate that, for almost all models on SRBench, preprocessing the input with PAN-SR improves performance (although there is a higher delta improvement for lower performing models compared to higher performing models).", "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": "We appreciate the reviewers' thoughtful, constructive, and positive feedback on our work. We are glad that the motivation behind our approach was found to be clearly presented (pCrp, gjaN, pZx6), and that our focus on minimizing false negatives (FNs) in variable selection (VS) for symbolic regression (SR) was considered both important and well justified (pCrp, gjaN). We're pleased that reviewers found our method to be novel and effective (gjaN, GHAd, pZx6), and that the empirical evaluation was described as solid, convincing, and comprehensive (gjaN, GHAd, pZx6). We also appreciate the recognition of PAN+SR's novel direction to improve scalability in SR (pZx6), as well as our contribution to SR benchmarking and software development (gjaN).\n\nWe'd like to take this opportunity to address your specific concerns and suggestions. We will include necessary revisions to clarify our work and strengthen the manuscript.\n\n### Concern 1: Scalability challenge being solved by PAN+SR\nWe agree that PAN+SR does not completely resolve the scalability challenge, but we believe it takes a significant step forward. In addition to the performance improvements presented, we'd like to share further concrete evidence, which will be included in the final paper for clarity.  \n\nFor context, the average values of $p$ and $p_0$ in the Feynman datasets are 186.15 and 3.65, respectively. On average, PAN reduces $p$ to 5.23 (97% reduction) in the best-case scenario (no noise), and to 63.1 (66% reduction) in the worst-case scenario (SNR=0.5). Furthermore, scalable SR methods such as TPSR [2] and DySymNet also benefit from PAN+SR, as shown in Figures 1 & 2. We chose not to include NeSymReS [1] as a baseline because TPSR [2] already employs a pre-trained NeSymReS backbone in our implementation.\n\n### Concern 2: Existence of a clear separation under noisy settings\nAs you rightly pointed out, separating signal from noise becomes harder in extremely noisy settings (Figures 5 & 6 in Appendix D.1). Nonetheless, relevant predictors still tend to cluster around the low-mean cluster, which explains the consistently high true positive rate (TPR) despite heavy noise (Figure 7 in Appendix D.2).  \n\n### Concern 3: PAN's role in improving performance\nIn Appendix D.2, we compared PAN against 4 other nonparametric VS methods on the high-dimensional Feynman database. As shown in Figures 7 and 8 of Appendix D.2, PAN consistently achieves the highest TPR across all settings, albeit at the cost of frequent false positives. Since identifying TPs is more critical in the SR pre-screening context, we believe PAN offers the most practical, effective, and safe solution.  \n\n### Concern 4: PAN's robustness to different sampling distributions for the predictors\nYour observation about the distributional assumptions in SRBench and the suggestion to consider [6] are also much appreciated. Indeed, each variable in the Feynman dataset is sampled from a uniform distribution, a common design choice in empirical studies aimed at ensuring even coverage of the input space to support the generalizability and robustness of the study. As [6] highlights, however, the sampling range of variables can influence SR performance. We agree with many points raised in [6] and find it to be a valuable reference for the related work section.\n\nHowever, we do not expect sampling range to impact the pre-screening performance of PAN+SR. This is because **tree-based methods like BART are invariant to monotonic transformations of the input features**. Thus, the pre-screening result should remain robust to variation in the sampling range.  \n\nAdditionally, while [6] enhances the Feynman dataset by adding 1-3 irrelevant variables, our experimental setup considers a **substantially more challenging high-dimensional regime**, adding $50p_0$ irrelevant variables (ranging **from 100 to 450**). We also note that while the original treatment of constants and integer-valued variables as real-valued variables does violate their physical meanings—as [6] rightly points out—this choice inadvertently increases the difficulty of the task, thereby providing a more stringent test for both pre-screening and SR modeling. For these reasons, we believe additional experiments based on [6] would not significantly change the conclusion of our study, though we acknowledge its importance and will include it in the discussion of future directions. \n\nHowever, motivated by other reviewers' comments, we will \n\n1. add an **ablation study to compare the effect of different clustering algorithms**, and\n2. include a **high correlation** experiment using the Friedman equation. See response to gjaN for details.\n\nWe believe these additions will strengthen the motivation behind PAN and broaden our experimental evaluation.\n \nThank you again for your time, constructive feedback, and thoughtful suggestions. We appreciate the opportunity to clarify our work and believe the revisions and additional discussions will help strengthen the paper.", "claims_and_evidence": "Claim: PAN+SR is a general purpose algorithm to increase SR scalability.\n\nComments: I'm not 100% convinced that PAN+SR \"solves\" the scalability problem but the authors present a novel direction of improvement for SR scalability which is pretty exciting. Previous work has generally handled this scalability challenge by (1) inducing programs with neural networks (which exposes an out of distribution problem) `[1, 2]` and (2) by using LLMs to induce programs  `[3, 4, 5]`  which are expensive to run. PAN+SR circumvents this problem by pre-selecting the variables but I'm not sure if a clean variable clustering exists in extremely noisy settings. Regardless, PAN+SR provides empirical justification that the model scales better than baseline models.\n\n\nClaim: PAN is the most performant pre-screening strategy.\n(Not really an explicit claim but an implicit one)\n\nComments: I generally agree that the algorithm is well motivated but it's not completely certain whether the increased performance is a result of pre-screening in general or pre-screening specifically with PAN+SR. Some additional experiments using other variable selection methods (e.g: BART) would be extremely insightful here.\n\n\n`[1]`: https://proceedings.mlr.press/v139/biggio21a/biggio21a.pdf\n\n`[2]`: https://github.com/deep-symbolic-mathematics/TPSR\n\n`[3]`: https://arxiv.org/abs/2404.18400\n\n`[4]`: https://arxiv.org/abs/2409.09359\n\n`[5]`: https://ai-2-ase.github.io/papers/52_InceptionSR_AAAI_25.pdf", "methods_and_evaluation_criteria": "The authors use a modified version of SRBench to simulate sampling equations with noisy irrelevant variables. These variables are drawn from the same data distribution to increase the task hardness. Overall, I found the evaluation criteria to be well motivated. \n\nOne small nitpick: The data for each variable in SRBench (ground-truth) is sampled from a normally distributed random variable. Would PAN+SR's performance suffer if the empirical data is not normally distributed?  Specifically, we know that physical laws tend to have very diverse data distributions `[6]`. e.g.: Newton's law and Coulomb's law have very similar equation sketches but the scale at which they operate is extremely different. Specifically, an additional experiment of PAN+SR's performance on `[6]` would be extremely helpful!\n\nOverall, I believe this paper will lead to great discussions at the conference and am in favor of **accepting** this paper. PAN+SR presents a new and refreshing direction for scaling SR methods by pre-screening the variable set using statistical heuristics that is empirically validated.\n\n\n`[6]`: https://arxiv.org/abs/2206.10540", "theoretical_claims": ".", "experimental_designs_or_analyses": "The methodology is sufficiently sound. The authors utilize and describe the evaluation setup proposed in SRBench.", "relation_to_broader_scientific_literature": ".", "essential_references_not_discussed": ".", "other_strengths_and_weaknesses": "I generally found the paper to be polished and easy to read. I think a bit more time can be devoted (in the appendix maybe) to go through how the pre-screening strategy would work on a small example, but otherwise this was a pretty interesting read.", "other_comments_or_suggestions": ".", "questions_for_authors": ".", "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": "1oy1DRuz9j", "reviewer": "Reviewer_GHAd", "summary": "The authors propose a feature selection preprocessing step to enhance the performance of symbolic regression algorithms. They introduce the method and evaluate its usage on SRBench across a number of algorithms and datasets.", "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": "We appreciate the reviewers' thoughtful, constructive, and positive feedback on our work. We are glad that the motivation behind our approach was found to be clearly presented (pCrp, gjaN, pZx6), and that our focus on minimizing false negatives (FNs) in variable selection (VS) for symbolic regression (SR) was considered both important and well justified (pCrp, gjaN). We're pleased that reviewers found our method to be novel and effective (gjaN, GHAd, pZx6), and that the empirical evaluation was described as solid, convincing, and comprehensive (gjaN, GHAd, pZx6). We also appreciate the recognition of PAN+SR's novel direction to improve scalability in SR (pZx6), as well as our contribution to SR benchmarking and software development (gjaN). \n\nWe'd like to take this opportunity to address your specific concerns and suggestions. We will include necessary revisions to clarify our work and strengthen the manuscript.\n\n### Concern 1: Motivation behind HAC\nWe fully agree with you that the usage of HAC is not fully motivated, and further justification is warranted in the main paper. We chose to skip this motivation primarily due to the page limit. Although not included in the paper, we tested several clustering algorithms, including HAC, $k$-means, affinity propagation, Gaussian mixture model (GMM), spectral clustering, mean shift, DBSCAN, and BIRCH. We found that **HAC consistently outperformed others**. A key reason for choosing HAC is its **robustness to class imbalance**, which is intrinsic to sparse regression problems where the number of relevant predictors ($p_0$) is much smaller than the number of irrelevant ones ($p-p_0$). In imbalanced data, large clusters can dominate centroid positions and overshadow the density signals of smaller clusters, making centroid-based and density-based methods less suitable for this task. In contrast, HAC starts with each point as its own cluster and merges based purely on pairwise distances. This allows small clusters (i.e., relevant predictors) to remain distinguishable. \n\nWe appreciate you highlighting this gap and will revise the manuscript to include a **dedicated explanation of our motivation and intuition in the method section**. We also plan to **add an ablation study** to compare the effect of different clustering algorithms on selection accuracy to justify our design choice further. We believe these additions will strengthen the motivation behind PAN and provide useful insights into the challenges posed by class imbalance in VS.\n\n### Concern 2: Performance under correlated structure (by pCrp and gjaN)\nWe appreciate your insightful comment on the role of correlated, confounding, and redundant variables in real-world settings. While our strategy of favoring false positives (FPs) over false negatives (FNs) may lead to more FPs in the presence of correlated predictors, we view this as a desired property. When uncertainty is high, it is preferable to retain a broader pool of potentially relevant features rather than risk excluding true signals. Furthermore, SR algorithms typically possess implicit regularization, which helps to further filter out irrelevant variables during model construction. \n\nAs you rightly noted, PAN+SR already **performs well on the real-world datasets**, where **27/35 (77%)** of datasets have at least 1 pair of predictors with **correlation >= 0.85**. \n\nIn addition, **motivated by your comment**, we conducted a pilot experiment using the Friedman equation under a **high correlation structure**. Predictors $x_1,\\ldots,x_p$ are drawn from a multivariate Uniform(0,1) distribution with an autocorrelation structure: $\\rho_{ij} = 0.9^{|i-j|}$, and the responses were generated as \n\n$$y = 10\\sin(\\pi x_1x_2) + 20(x_3 – 0.5)^2 + 10x_4 + 5x_5 + \\varepsilon, \\qquad\\varepsilon \\sim N(0,\\sigma^2).$$ \n\nWe fixed $n=1000$, $p=100$, and ran 100 trials per SNR level:  \n\n| Metric | SNR=10 | SNR=5 | SNR=1 |\n|-|-|-|-|\n| TPR    | 100%   | 100%   | 99.8%  |\n| FPR    | 15.98% | 23.26% | 33.23% |\n\nEven so, **PAN consistently identified all 5 relevant predictors**, except for 1 run under SNR=1. We plan to include the discussion of correlated structure in the paper to help motivate future research focused on more realistic simulation settings. \n\nWe are also more than happy to incorporate your suggestions on redefining PAN outside of the abstract and to motivate the difficulty of variable selection using the provided reference. Thank you again for your time, constructive feedback, and thoughtful suggestions. We appreciate the opportunity to clarify our work and believe the revisions and additional discussions will help strengthen the paper.", "claims_and_evidence": "The authors do a good job overall of making evidence-based claims. Evidence is convincing as it is over multiple algorithms, datasets, levels of noise, across metrics, and in comparison to other feature selection strategies (appendix).", "methods_and_evaluation_criteria": "Overall the methods make sense and the evaluation is commendable - the authors test their algorithm in combination with a number of SR algorithms and across real-world and synthetic datasets in a robust way. They also compare their feature selection strategy with others in terms of false positive/negative rates (of feature selection) with different levels of noise. \n\nThe only part of the methodology that didn't make total sense to me was the use of hierarchical agglomerative clustering to partition the feature space into two groups. It seems like overkill, and it isn't fully motivated. if you're just clustering the feature ranks over many trials, can't you use a much simpler threshold finder than HAC?", "theoretical_claims": "There are no proofs", "experimental_designs_or_analyses": "Authors appear to follow recommended benchmarking practices. There aren't statistical tests for differences but the effect sizes are reported and laid out visually.", "relation_to_broader_scientific_literature": "Most work in SR develops a new method and then reports its benchmark results on SRBench in comparison to previous results. The authors instead evaluate a broadly applicable preprocessing strategy that synergizes with many methods. \n\nThe relevant literature is appropriately cited in symbolic regression, and reviewed in terms of feature selection as well.", "essential_references_not_discussed": "There are some SR-based feature selection research papers that might be worth including but it is ancillary.", "other_strengths_and_weaknesses": "Overall I thought the paper was well done. \n\nOne weakness, in general, is that the strategy of preferring False Positives to false negatives may be less optimal on real-world data where there is a lot of cross-correlation, confounding, and redundancy in the dataset. It's a bit easier to optimize for that criteria on synthetic benchmarks of known physical systems than it might be for unknown systems. The comparison on real-world datasets mitigates some of this worry but it is worth mentioning.", "other_comments_or_suggestions": "- redefine PAN outside of the abstract\n- the authors motivate the problem by mentioning SR is NP-hard, but they might as well also note that feature selection is hard as well (http://www.jmlr.org/proceedings/papers/v40/Foster15.pdf)", "questions_for_authors": "None apart from those mentioned thus far", "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": "hx5zaltZCn", "reviewer": "Reviewer_gjaN", "summary": "This paper proposes a rank-clustering PAN strategy for screening for relevant features before running symbolic regression (SR) methods on very high dimensional data, where the goal is to minimize the false negative rate (avoiding missing important variables). The idea is to repeatedly run BART and use the ranking of feature importance measures to cluster relevant and irrelevant features. Then, only the features identified as relevant are fed into a symbolic regression algorithm. Through extensive simulation studies, this strategy is shown to be effective in improving the performance and reduce the running time of all popular downstream SR algorithms.", "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": "We appreciate the reviewers' thoughtful, constructive, and positive feedback on our work. We are glad that the motivation behind our approach was found to be clearly presented (pCrp, gjaN, pZx6), and that our focus on minimizing false negatives (FNs) in variable selection (VS) for symbolic regression (SR) was considered both important and well justified (pCrp, gjaN). We're pleased that reviewers found our method to be novel and effective (gjaN, GHAd, pZx6), and that the empirical evaluation was described as solid, convincing, and comprehensive (gjaN, GHAd, pZx6). We also appreciate the recognition of PAN+SR's novel direction to improve scalability in SR (pZx6), as well as our contribution to SR benchmarking and software development (gjaN). \n\nWe'd like to take this opportunity to address your specific concerns and suggestions. We will include necessary revisions to clarify our work and strengthen the manuscript.\n\n### Concern 1: Motivation and theoretical justification for the proposed method (PAN)\nWe agree that the motivation for PAN can be more clearly articulated. PAN is motivated by the observation that VIP rankings, unlike raw VIP values, tend to exhibit a clear **2-group structure** separating relevant from irrelevant predictors. This **2-group structure** is obscured in raw VIPs due to their bounded nature (summing to 1), whereas the ranking scale is wider and more interpretable. In Section 4 (paragraphs 5-8) and Appendix D.1, we show that under the uniform assumption, the expected mean rank $\\bar{r}_{j\\cdot}$ equals $(1+p_0)/2$ for relevant predictors and $(p_0+1+p)/2$ for irrelevant ones, naturally forming 2 clusters. \n\nThis led us to **frame VS as a clustering problem** over mean VIP rankings $\\bar{r}_{j\\cdot}$. Although not included in the paper, we tested several clustering algorithms, including HAC, $k$-means, affinity propagation, Gaussian mixture model (GMM), spectral clustering, mean shift, DBSCAN, and BIRCH. We found that **HAC consistently outperformed others**. A key reason is its **robustness to class imbalance**, which is intrinsic to sparse regression problems where $p_0 \\ll p$. In imbalanced data, large clusters can dominate centroid positions and overshadow the density signals of smaller clusters, making centroid-based and density-based methods less suitable for this task. In contrast, HAC starts with each point as its own cluster and merges based purely on pairwise distances. This allows small clusters (i.e., relevant predictors) to remain distinguishable. \n\nWe plan to add an ablation study to compare the effect of different clustering algorithms on selection accuracy to justify our design choice further. We believe these additions will strengthen the motivation behind PAN and provide useful insights into the challenges posed by class imbalance in VS.\n\n### Concern 2: Performance under correlated predictors (also raised by pCrp and GHAd)\nWe appreciate your suggestion regarding correlation in the predictors. While the standard practice in SR benchmarks focuses on independence among the predictors, we agree that correlated inputs represent a realistic and important challenge. That said, our current experimental design already introduces several challenges **beyond existing SR benchmarks**:\n\n1. high dimensionality\n2. irrelevant features\n3. high output noises (8 levels)\n4. low sample sizes (4 levels)\n\nleading to $8 \\cdot 4 \\cdot 100 \\cdot 10=\\textbf{32000}$ simulation settings. Although we're not able to test all suggested settings, they are not uncovered by our work. First, it's worth noting that PAN+SR already **performs well on the real-world datasets** (see Figure 1), where **27/35 (77%)** of datasets have at least 1 pair of predictors with **correlation >= 0.85**. \n\nIn addition, **motivated by your comment**, we conducted a pilot experiment using the Friedman equation under a **high correlation structure**. Predictors $x_1,\\ldots,x_p$ are drawn from a multivariate Uniform(0,1) with an autocorrelation structure: $\\rho_{ij} = 0.9^{|i-j|}$, and the responses were generated as \n\n$$y = 10\\sin(\\pi x_1x_2) + 20(x_3 – 0.5)^2 + 10x_4 + 5x_5 + \\varepsilon, \\qquad\\varepsilon \\sim N(0,\\sigma^2).$$ \n\nWe fixed $n=1000$, $p=100$, and ran 100 trials per SNR level:\n\n| Metric | SNR=10 | SNR=5 | SNR=1 |\n|-|-|-|-|\n| TPR    | 100%   | 100%   | 99.8%  |\n| FPR    | 15.98% | 23.26% | 33.23% |\n\nEven so, **PAN consistently identified all 5 relevant predictors**, except for 1 run under SNR=1. We plan to include the discussion of correlated structure in the paper to help motivate future research focused on more realistic simulation settings.\n\nThank you again for your time, constructive feedback, and thoughtful suggestions. We appreciate the opportunity to clarify our work and believe the revisions and additional discussions will help strengthen the paper.", "claims_and_evidence": "The claims for the benefits of the proposed method are mostly supported by empirical evidence. It would be nice to have some theoretical insights.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "This submission does not contain much theoretical aspects.", "experimental_designs_or_analyses": "I have checked the experimental setups, including data generating process, algorithms used, and their implementations. Existing experiments are solid, but it would further benefit from a more comprehensive design, which I elaborate in the \"Questions\" section.", "relation_to_broader_scientific_literature": "The key contributions include:\n- Formalize the goal of variable selection in SR (minimizing FNR) and its distinction from standard variable selection.\n- Propose a strategy to conduct variable selection, which boosts existing methods, and may also inspire future development in this direction.\n- Add to the software and benchmark development of the SR field.", "essential_references_not_discussed": "I haven't found such a case.", "other_strengths_and_weaknesses": "Please see my comments and questions.", "other_comments_or_suggestions": "- While the problem considered is well motivated and the argument on limiting FNR is solid, the proposed solution seems rather heuristic. It is not clear why the authors settle down to this specific solution. Some theoretical justifications will be helpful in establishing this method. \n- The experiments are limited to i.i.d. covariates where the importance of features can be readily assessed by ranking and clustering. However, this kind of ranking strategy might be problematic when features are correlated. It will be helpful to extend the experiments.", "questions_for_authors": "Similar to my comments, \n- Is there theoretical justification for the method? Why do you pick this specific strategy?\n- How do the methods work for correlated features?", "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": "A1yiujL2O0", "reviewer": "Reviewer_pCrp", "summary": "This paper proposes a variable selection method for input variables related to the output. The proposed method is used for data preprocessing in symbolic regression and can improve the accuracy and speed of symbolic regression. Specifically, the authors propose a method called PAN+SR, which combines a key idea of ab initio nonparametric variable selection with SR to efficiently pre-screen large input spaces and reduce search complexity while maintaining accuracy. In the paper, the authors emphasize the importance of FNR, as any erroneous exclusion of related variables will lead to the failure of subsequent symbolic regression. The authors conducted numerous experiments in the paper to demonstrate the effectiveness of the proposed method, including the improvement of regression accuracy and the robustness of the method to noise.", "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": "We appreciate the reviewers' thoughtful, constructive, and positive feedback on our work. We are glad that the motivation behind our approach was found to be clearly presented (pCrp, gjaN, pZx6), and that our focus on minimizing false negatives (FNs) in variable selection (VS) for symbolic regression (SR) was considered both important and well justified (pCrp, gjaN). We're pleased that reviewers found our method to be novel and effective (gjaN, GHAd, pZx6), and that the empirical evaluation was described as solid, convincing, and comprehensive (gjaN, GHAd, pZx6). We also appreciate the recognition of PAN+SR's novel direction to improve scalability in SR (pZx6), as well as our contribution to SR benchmarking and software development (gjaN). \n\nWe'd like to take this opportunity to address your specific concerns and suggestions. We will include necessary revisions to clarify our work and strengthen the manuscript.\n\n### Concern 1: Lack of comparison with other VS methods\nContrary to the impression that no comparisons were made, our paper includes a **comprehensive evaluation of the proposed VS method (PAN) against 4 strong nonparametric VS methods**. These results are presented in **Appendix D.2** and referenced in **Section 4 paragraph 3**, where we note PAN consistently achieves the highest true positive rate (TPR) across all settings. We'll make this pointer more prominent in the revision to avoid confusion. \n\nSecond, while the suggested papers share some surface-level similarities with our work and are worth citing, their focus is meaningfully different. Specifically, they use genetic programming (GP)—an SR method—as a tool for VS. In contrast, our work develops a **model-free VS** method that serves as a pre-screening step to improve the scalability and performance of SR methods, including but not limited to GP-based SR. \n\nMoreover, GP-based VS methods have several scalability and expressive power limitations. In addition to the number of input variables $p$, their effectiveness and tractability also depend heavily on: \n\n1. The complexity and size of the operator set used for expression constructions\n2. The maximum tree depth allowed\n\nAs these grow, GP's search space **expands combinatorially**, making it less scalable and less suitable for high-dimensional problems. Furthermore, the reliance on a **pre-defined operator set** makes them inherently not model-free, while our proposed approach is designed to be **modular, scalable, and model-free**. \n\nFinally, we respectfully disagree with the characterization that \"the contribution of the paper lies in proposing a variable selection algorithm.\" Our work goes beyond this—we address scaling SR to extreme dimensions, where the synergy between VS and SR is central. Our key message, supported by comprehensive experiments, is that a carefully designed model-free VS step consistently improves downstream SR across diverse methods. \n\n### Concern 2: Performance under correlated and/or noisy predictors (also raised by gjaN and GHAd).  \nWe also appreciate your thoughtful suggestion regarding noise and correlation in predictors. While standard practice in SR benchmarks focuses on no/low noise in the output and independence in the predictors, we agree that noisy and/or correlated inputs represent a realistic and important challenge. That said, our current experimental design already introduces several challenges **beyond existing SR benchmarks**: \n\n1. high dimensionality\n2. irrelevant features\n3. high output noises (8 levels)\n4. low sample sizes (4 levels)\n\nleading to $8 \\cdot 4 \\cdot 100 \\cdot 10=\\textbf{32000}$ simulation settings. Although we're not able to test all suggested settings, they're not uncovered by our work. First, it's worth noting that PAN+SR already **performs well on the real-world datasets** (see Figure 1), where **27/35 (77%)** of datasets have at least 1 pair of predictors with **correlation >= 0.85**. \n\nIn addition, **motivated by your comment**, we conducted a pilot experiment using the Friedman equation (also used in [3]) under a **high correlation structure**. Predictors $x_1,\\ldots,x_p$ are drawn from a multivariate Uniform(0,1) distribution with an autocorrelation structure: $\\rho_{ij} = 0.9^{|i-j|}$, and the responses were generated as \n\n$$y = 10\\sin(\\pi x_1x_2) + 20(x_3 – 0.5)^2 + 10x_4 + 5x_5 + \\varepsilon, \\qquad\\varepsilon \\sim N(0,\\sigma^2).$$ \n\nWe fixed $n=1000$, $p=100$, and ran 100 trials per SNR level:\n\n| Metric | SNR=10 | SNR=5 | SNR=1 |\n|-|-|-|-|\n| TPR    | 100%   | 100%   | 99.8%  |\n| FPR    | 15.98% | 23.26% | 33.23% |\n\nEven so, **PAN consistently identified all 5 relevant predictors**, except for 1 run under SNR=1. We plan to include the discussion of correlated structure in the paper to help motivate future research focused on more realistic simulation settings. \n\nThank you again for your time and constructive feedback. We appreciate the opportunity to clarify our work and believe these revisions will help strengthen the paper.", "claims_and_evidence": "The supporting evidence is insufficient, as the data used in the paper is self constructed and lacks comparison with similar methods.", "methods_and_evaluation_criteria": "Yes，The proposed methods and/or evaluation criteria (e.g., benchmark datasets) make sense for the problem or application at hand.", "theoretical_claims": "The paper has no theoretical proof, and its contribution lies in the design of the algorithm. However, I think the author's emphasis on the importance of low FNR is correct.", "experimental_designs_or_analyses": "There are some issues with the experimental design of the paper. Firstly, the contribution of the paper lies in proposing a variable selection method, rather than a symbolic regression algorithm. The paper lacks comparison with other variable selection algorithms. In addition, in the noise robustness experiment, the paper only considered the presence of noise in the output variable y, which is inconsistent with reality. In practice, both input and output variables may have noise.", "relation_to_broader_scientific_literature": "The key contributions of the paper are related to the correlation analysis and information theory in other literature. In previous studies, correlation coefficients and mutual information between variables were often used to determine whether variables were correlated.", "essential_references_not_discussed": "The contribution of the paper lies in proposing a variable selection algorithm rather than a symbolic regression algorithm, but the author only discussed the symbolic regression algorithm in Related Works and did not discuss other variable selection algorithms. Moreover, there was no comparison with other variable selection algorithms in the experiment, which is puzzling. Here are three relevant literature for the author's reference.\n\n1.Q. Chen, M. Zhang and B. Xue, \"Feature Selection to Improve Generalization of Genetic Programming for High-Dimensional Symbolic Regression,\" in IEEE Transactions on Evolutionary Computation, vol. 21, no. 5, pp. 792-806, Oct. 2017, doi: 10.1109/TEVC.2017.2683489.\n\n2.Q. Chen, B. Xue, B. Niu and M. Zhang, \"Improving generalisation of genetic programming for high-dimensional symbolic regression with feature selection,\" 2016 IEEE Congress on Evolutionary Computation (CEC), Vancouver, BC, Canada, 2016, pp. 3793-3800, doi: 10.1109/CEC.2016.7744270.\n\n3.B. Al-Helali, Q. Chen, B. Xue and M. Zhang, \"Genetic Programming for Feature Selection Based on Feature Removal Impact in High-Dimensional Symbolic Regression,\" in IEEE Transactions on Emerging Topics in Computational Intelligence, vol. 8, no. 3, pp. 2269-2282, June 2024, doi: 10.1109/TETCI.2024.3369407.", "other_strengths_and_weaknesses": "The paper has done a great job, with standardized writing and clear discourse. However, there are the following suggestions for the author to further improve.\n1.The paper constructed a new high-dimensional SR dataset for testing variable selection, which is a good work. However, the current constructed dataset is still too simple. It is possible to consider constructing some more difficult redundant variables, such as constructing new redundant variables like x3=x1+x2.\n2. The paper only considered adding noise to the output variable y in the experiment, and it is recommended that the author also add noise to the input variable in the future.", "other_comments_or_suggestions": "All comments and suggestions are provided above.", "questions_for_authors": "All comments and suggestions are provided above. Here, I would like to emphasize the two issues that concern me the most: first, the comparison with other variable selection methods, and second, adding noise to the input variables as well.", "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: Overall, the reviewers showed positive impressions for this work and offered 1 accept and 3 weak accepts.

Many of the reviewers valued the proposed method and evaluation results that support the benefits of using the method.
To the best of my knowledge, it is rare to see discussions focused on feature selection in symbolic regression tasks as it's usually part of proposed SR methods. The proposed method is method agnostic as it's a preprocessing approach, and I believe such discussions should be very important for the community.

While I support accepting this work for ICML 2025, I want to make the following suggestions (2 and 3 are based on our internal discussion):
1. Update references section. I found many of papers cited as arXiv preprints published at some venues 
2. Add an explicit section in the appendix that explains the dataset construction and sampling process, especially for irrelevant data e.g., how exactly the irrelevant features were generated such as range of values, distributions, whether or not the range is overlapped with those of the original features
3. In the camera-ready or future studies, consider using SRSD datasets (dummy variables ver.) for evaluations or at least follow their procedure to generate dummy (irrelevant) variables as it's more challenging (much larger range and independent from a given equation)
Decision: Accept (poster)