系数校准作为符号回归中的选择压力

# Coefficient Calibration as Selection Pressure in Symbolic Regression
Source: [https://arxiv.org/html/2610.10931](https://arxiv.org/html/2610.10931)
Veronica Guidettiand Federica Mandreoli††thanks:We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support\.††thanks:M\. Billa, V\. Guidetti, and F\. Mandreoli are with the Department of Physics, Informatics and Mathematics, University of Modena and Reggio Emilia, 41125 Modena, Italy \(e\-mail: mattia\.billa@unimore\.it; veronica\.guidetti@unimore\.it; federica\.mandreoli@unimore\.it\)\.

###### Abstract

In memetic symbolic regression, candidate structures are compared after coefficient calibration, so the calibration protocol itself contributes to evolutionary selection\. Standard centralized calibration evaluates each structure at its pooled\-sample optimum, ignoring how stable this calibration is under covariate shifts, and can thus favor structures whose fit relies on sample\-specific coefficients\.

We propose Dirichlet\-Sinkhorn Constant Averaging \(DSCA\), a calibration strategy that partitions the optimization data into equally sized subsets with different covariate distributions, calibrates each candidate independently on every partition, and evaluates it at the mean of the resulting parameters\. We show that the excess loss of DSCA relative to centralized calibration vanishes at the population level for correctly specified, identifiable expressions, whereas under misspecification it persists when partition\-specific calibrations do not aggregate to the pooled optimum\.

On synthetic benchmarks and ten real\-world datasets, DSCA improves functional recovery and the accuracy\-complexity trade\-off over centralized Broyden–Fletcher–Goldfarb–Shanno and Levenberg\-Marquardt calibration, under selection by negative log\-likelihood and by the Akaike and Bayesian information criteria\. Mechanism analyses associate the DSCA excess loss with the generalization gap and show that the effect is not reproduced by repeated centralized fitting\. These results indicate that controlled heterogeneous calibration provides a complementary source of selection pressure in symbolic\-regression search\.

###### Index Terms:

Genetic Programming, Symbolic Regression, Optimization

## IIntroduction

Symbolic regression \(SR\) aims to discover mathematical models that capture the underlying relationships within data, prioritizing both predictive accuracy and structural interpretability\[[1](https://arxiv.org/html/2610.10931#bib.bib1),[2](https://arxiv.org/html/2610.10931#bib.bib2)\]\. Unlike traditional regression, which optimizes parameters for a predefined architecture, SR involves a bi\-level optimization task\. The algorithm must simultaneously search the discrete space of algebraic structures and the continuous space of numerical coefficients\. Genetic programming represents the most prominent framework for this task, utilizing evolutionary mechanisms to navigate the combinatorial space of mathematical expressions\[[3](https://arxiv.org/html/2610.10931#bib.bib3)\]\.

Within this bi\-level framework, the identification of optimal numerical constants acts as a critical bottleneck\. Historically, genetic programming relied on evolutionary variation to adapt constants, a process that is inefficient and often results in premature convergence or structural bloat\[[4](https://arxiv.org/html/2610.10931#bib.bib4),[5](https://arxiv.org/html/2610.10931#bib.bib5)\]\. To alleviate this, modern approaches frequently hybridize the evolutionary search with continuous local optimizers, such as the Levenberg\-Marquardt or Broyden\-Fletcher\-Goldfarb\-Shanno algorithms\[[6](https://arxiv.org/html/2610.10931#bib.bib6),[7](https://arxiv.org/html/2610.10931#bib.bib7),[8](https://arxiv.org/html/2610.10931#bib.bib8)\]\.

Accurate coefficient calibration substantially improves the evaluation of a fixed candidate structure, but it also affects the outer structural search\. In a memetic SR system, candidate expressions are compared only after their numerical parameters have been adapted to the observed training sample\. Consequently, flexible or misspecified structures may become competitive by exploiting sample\-specific variation through finely tuned coefficients\. This effect is particularly relevant in SR, where most structures visited during evolution are necessarily imperfect approximations of the underlying relationship and coefficient estimation is frequently affected by ill\-conditioning\[[9](https://arxiv.org/html/2610.10931#bib.bib9),[10](https://arxiv.org/html/2610.10931#bib.bib10)\]\.

The resulting generalization\-complexity trade\-off is typically addressed at the selection level, through explicit parsimony pressure or information\- and description\-based criteria\[[11](https://arxiv.org/html/2610.10931#bib.bib11)\]\. Such mechanisms determine how candidate expressions are compared after their numerical coefficients have already been calibrated on the observed sample\. Their penalized likelihood structure usually depends linearly on the number of parameters used by the model, a structure which is not suitable for misspecified models typical of SR and may not be sufficient to guide evolution toward the right accuracy\-complexity trade\-off\[[12](https://arxiv.org/html/2610.10931#bib.bib12),[13](https://arxiv.org/html/2610.10931#bib.bib13)\]\.

Here we investigate a complementary intervention: we modify the coefficient\-calibration stage through which an individual receives the parameter values used for its fitness evaluation\. The objective is therefore not to improve the training optimum itself, but to make the evaluation of candidate structures less favorable to solutions whose empirical performance relies on highly sample\-specific coefficient calibration\.

To achieve this, we draw on ideas from distributed and federated parameter estimation, where local models are fitted on distinct data sources and subsequently aggregated into a common parameterization\[[14](https://arxiv.org/html/2610.10931#bib.bib14)\]\. In those settings, statistical heterogeneity is generally a constraint that the learning procedure must accommodate\. In our centralized setting, we instead deliberately induce controlled heterogeneity and exploit it during candidate evaluation\. Intuitively, an expression whose good empirical fit depends on finely tuned, sample\-specific coefficients should be less compatible with parameter estimates obtained across different views of the covariate space\. Aggregating such local estimates can therefore reduce the competitiveness of these candidates during evolutionary selection, providing an additional source of regularization against overfitting that is complementary to the outer model\-selection criterion\.

### I\-AOur Contributions

Our contributions are threefold:

- •A new calibration strategy for SR\.We introduce Dirichlet\-Sinkhorn Constant Averaging \(DSCA\), which replaces centralized coefficient calibration with calibration over controlled covariate\-heterogeneous partitions followed by parameter averaging\. Because DSCA modifies only the individual\-evaluation stage, it can be readily integrated into standard SR pipelines without changing the underlying search procedure or outer model\-selection criterion\.
- •A theoretical characterization of its selection effect\.We characterize the excess loss induced by DSCA relative to centralized calibration and show how this additional calibration\-stage penalty can alter evolutionary selection, providing an implicit regularization mechanism\. At population level, the excess vanishes for correctly specified and identifiable expressions, whereas under misspecification it can persist when environment\-specific calibrations do not aggregate to the pooled optimum\. We further characterize its dependence on covariate heterogeneity and its finite\-sample behavior\.
- •An empirical evaluation of effectiveness and mechanism\.Across controlled synthetic benchmarks and ten real\-world datasets, we compare DSCA with standard centralized calibration methods under multiple model\-selection criteria\. DSCA improves functional recovery and the accuracy–complexity trade\-off, while mechanism analyses show that its excess loss is associated with the generalization gap and that the observed effect is not reproduced by repeated centralized fitting\.

The remainder of this article is organized as follows\. Section[II](https://arxiv.org/html/2610.10931#S2)reviews memetic approaches, model selection in SR, and distributed inference\. Section[III](https://arxiv.org/html/2610.10931#S3)presents DSCA, its partitioning and averaging steps, and an analysis of the excess loss it induces\. Section[IV](https://arxiv.org/html/2610.10931#S4)describes the experimental setup, and Section[V](https://arxiv.org/html/2610.10931#S5)reports results on the synthetic and real\-world tracks\. Finally, Section[VI](https://arxiv.org/html/2610.10931#S6)concludes and outlines directions for future work\.

## IIRelated Work

### II\-AMemetic Approaches and Parameter Optimization in Symbolic Regression

Discovering SR models requires navigating two interrelated search spaces\. The algorithm must identify an appropriate algebraic structure and simultaneously determine the optimal numerical coefficients for that structure\. In traditional genetic programming, parameters were often initialized as ephemeral random constants that remained fixed throughout the evolutionary process\[[4](https://arxiv.org/html/2610.10931#bib.bib4),[15](https://arxiv.org/html/2610.10931#bib.bib15),[16](https://arxiv.org/html/2610.10931#bib.bib16),[5](https://arxiv.org/html/2610.10931#bib.bib5)\]\. This paradigm limited the precision of the resulting models, as the algorithm was forced to rely on the stochastic combination of existing constants via crossover to approximate optimal values\.

To overcome these structural limitations, researchers introduced memetic algorithms, which hybridize population\-based evolutionary search with localized continuous optimization strategies\[[17](https://arxiv.org/html/2610.10931#bib.bib17),[18](https://arxiv.org/html/2610.10931#bib.bib18)\]\. Topchy and Punch provided foundational contributions to this area by demonstrating that embedding a local gradient descent search for numeric leaf values significantly accelerated the convergence of genetic programming\[[6](https://arxiv.org/html/2610.10931#bib.bib6)\]\. Their work also revealed a critical secondary effect: the integration of local search inherently steered the global evolutionary process toward algebraic structures that were mathematically more amenable to parametric tuning\.

Subsequent advancements have formalized the application of continuous optimization within discrete expression trees\. Kommenda and colleagues conducted extensive analyses on the use of nonlinear least squares techniques, specifically employing the Levenberg\-Marquardt algorithm\[[19](https://arxiv.org/html/2610.10931#bib.bib19),[15](https://arxiv.org/html/2610.10931#bib.bib15),[16](https://arxiv.org/html/2610.10931#bib.bib16),[7](https://arxiv.org/html/2610.10931#bib.bib7)\]\. To supply the optimizer with exact gradient information, they incorporated automatic differentiation, a methodology that proved faster and more accurate than standard numerical or symbolic differentiation\. This approach has become a core component in high\-performance genetic programming frameworks such as Operon, which relies on exact function derivatives for local parameter identification\[[20](https://arxiv.org/html/2610.10931#bib.bib20),[9](https://arxiv.org/html/2610.10931#bib.bib9)\]\. Other recent approaches, such as PySR\[[8](https://arxiv.org/html/2610.10931#bib.bib8)\], formalize this process into a dedicated evolve, simplify, and optimize cycle that periodically refines constants via algorithms like BFGS\[[21](https://arxiv.org/html/2610.10931#bib.bib21)\]\. Concurrently, researchers have tackled this problem using Bayesian frameworks to handle both model selection and uncertainty estimation\[[22](https://arxiv.org/html/2610.10931#bib.bib22),[23](https://arxiv.org/html/2610.10931#bib.bib23),[24](https://arxiv.org/html/2610.10931#bib.bib24)\]\.

Recent comparative studies have evaluated a wider array of continuous optimizers in the context of SR\. For example, dos Reis et al\. evaluated multiple constant optimization methods, noting that no single technique is universally optimal across all problem types\[[25](https://arxiv.org/html/2610.10931#bib.bib25)\]\. Furthermore, Harrison et al\. systematically compared first\- and second\-order gradient\-based coefficient optimization in SR and found second\-order approaches to be consistently more effective, despite their higher per\-step computational cost\[[10](https://arxiv.org/html/2610.10931#bib.bib10)\]\. This distinction is particularly relevant in SR, where individual expressions generally contain relatively few numerical parameters and nonlinear least\-squares or quasi\-Newton methods are therefore computationally viable\. Recent benchmarking further shows that the quality of coefficient optimization can itself affect the structures favored by evolution: poorly calibrated candidates may remain artificially short because introducing additional constants provides little fitness advantage when those constants cannot be optimized effectively\[[25](https://arxiv.org/html/2610.10931#bib.bib25)\]\.

However, integrating continuous optimization into genetic programming presents algorithmic difficulties\. Kronberger highlighted that the evolutionary process frequently generates overparameterized expressions containing mathematically redundant constants\[[9](https://arxiv.org/html/2610.10931#bib.bib9)\]\. This overparameterization leads to rank\-deficient and ill\-conditioned Jacobian matrices during the nonlinear least squares optimization process, which can degrade optimizer efficiency and reliability\. These observations emphasize that coefficient calibration and structural search cannot be regarded as independent components: the quality and stability of the inner optimization directly influence which expressions remain competitive during evolution\.

### II\-BModel Selection and Complexity Control in Symbolic Regression

Effective coefficient calibration does not by itself resolve the model\-selection problem\. Improving the fit of the numerical parameters makes a given expression more competitive, but does not determine whether its structure provides a useful approximation of the underlying relationship\. Once their coefficients have been calibrated, symbolic expressions must therefore still be compared according to both empirical fit and structural complexity\.

Genetic programming has traditionally addressed this problem through explicit parsimony pressure and multi\-objective accuracy\-complexity optimization, while more recent work has investigated likelihood\- and description\-based selection criteria such as AIC, AICc, BIC, and MDL\[[26](https://arxiv.org/html/2610.10931#bib.bib26),[27](https://arxiv.org/html/2610.10931#bib.bib27),[11](https://arxiv.org/html/2610.10931#bib.bib11),[28](https://arxiv.org/html/2610.10931#bib.bib28)\]\. These criteria act at the outer selection stage by correcting empirical fit according to some notion of model complexity\. AIC and BIC use the number of fitted parameters through their familiar2​p2pandp​log⁡Np\\log Ncorrections, whereas recent MDL formulations for SR additionally account for structural complexity and coefficient precision through an evidence\-related description\-length objective\.

The standard forms of these corrections nevertheless rely on assumptions that are not automatically satisfied by symbolic models\. In particular, parameter\-count penalties such as those underlying AIC and BIC are derived under regular parametric asymptotics, including local identifiability and nonsingular Fisher information\. Symbolic expressions can violate these conditions through redundant coefficients, equivalent parameterizations, inactive components, and highly correlated or weakly identified parameter directions\[[9](https://arxiv.org/html/2610.10931#bib.bib9)\]\. In singular models, nominal parameter count no longer provides an adequate measure of effective statistical complexity; for evidence\-based quantities, the standard Laplace expansion is instead replaced by geometry\-dependent corrections such as those involving the real log canonical threshold\[[29](https://arxiv.org/html/2610.10931#bib.bib29),[30](https://arxiv.org/html/2610.10931#bib.bib30)\]\.

Misspecification raises a distinct problem\. The model\-identification and evidence\-based interpretations underlying criteria such as BIC and Bayesian or description\-length model comparison are most naturally formulated in an*M\-closed*setting, in which the data\-generating distribution is represented within the candidate model class\[[31](https://arxiv.org/html/2610.10931#bib.bib31)\]\. Symbolic\-regression search is instead naturally*M\-open*: most individuals encountered during evolution are misspecified approximations, and the true relationship need not belong to the symbolic hypothesis space, let alone to the finite population available at a particular generation\. AIC is less directly tied to this assumption because its target is predictive risk, but its standard parameter\-count correction still relies on regular asymptotics\.

Despite these limitations, selection criteria have proved useful in SR, particularly for controlling formula bloat and favoring parsimonious expressions\. AIC and BIC have been successfully used to guide evolutionary selection, while recent experiments with richer MDL\-based penalties show that their stronger complexity pressure can also lead to premature convergence or excessively simple expressions when used directly as evolutionary fitness\[[11](https://arxiv.org/html/2610.10931#bib.bib11),[28](https://arxiv.org/html/2610.10931#bib.bib28)\]\. More generally, these criteria have been more directly effective at controlling the complexity of the formulas retained by the search than at preventing misspecified, sample\-sensitive structures from becoming competitive after coefficient calibration\. DSCA addresses this complementary mechanism: rather than adding another outer complexity penalty, it modifies the coefficient calibration through which an individual is evaluated and can therefore be combined with AIC, BIC, MDL, or other selection criteria\.

### II\-CDistributed Inference, Federated Learning, and Subset\-Based Model Aggregation

Partitioning large datasets to compute local estimators independently and subsequently aggregating the results has emerged as a highly scalable paradigm for statistical learning\. Theoretical analyses demonstrate that averaging local solutions, such as empirical risk minimizers, can achieve the minimax optimal statistical convergence rates of a centralized model, provided the number of data partitions is bounded appropriately\[[32](https://arxiv.org/html/2610.10931#bib.bib32),[33](https://arxiv.org/html/2610.10931#bib.bib33),[34](https://arxiv.org/html/2610.10931#bib.bib34)\]\. The central idea is that parameter information extracted independently from different subsets can be combined into a single global estimate, provided that the local estimators are sufficiently compatible\.

This mathematical architecture forms the foundation of federated learning, where the objective is to collaboratively train a global model across decentralized clients while strictly preserving data privacy\[[14](https://arxiv.org/html/2610.10931#bib.bib14)\]\. A critical challenge in federated networks is the presence of non\-identically distributed data across participating clients\[[35](https://arxiv.org/html/2610.10931#bib.bib35)\]\. To benchmark algorithmic resilience to statistical heterogeneity, researchers often synthesize non\-ideal partitions using Dirichlet distributions\. The concentration parameter of the Dirichlet distribution provides a convenient way to control the degree of statistical heterogeneity of feature and label among simulated clients\.

Translating these principles to SR, early federated frameworks aggregated fitness across clients but lacked global parameter estimation, restricting models to evaluating randomly sampled formulas with fixed ephemeral random constants\[[36](https://arxiv.org/html/2610.10931#bib.bib36)\]\. To resolve this fundamental limitation, Billa et al\. developed Bayesian Federated Symbolic Regression, which fuses local Gaussian posteriors for marginal likelihood\-based model selection without sharing raw data\[[24](https://arxiv.org/html/2610.10931#bib.bib24)\]\. DSCA starts from the same general idea of estimating numerical parameters locally and combining information across heterogeneous data sources, but uses it for a different purpose\. In a genuine federated problem, heterogeneity is imposed by the participating clients and aggregation seeks to recover an effective common model despite that heterogeneity\. In our centralized setting, the partitions can instead be designed deliberately\. We exploit this freedom to construct subsets with different covariate marginals while preserving the conditional distribution of the response, and use arithmetic aggregation not as a communication\-efficient approximation to centralized optimization, but as a controlled perturbation of candidate evaluation during evolutionary selection\.

## IIIMethodology

The proposed approach, termed Dirichlet\-Sinkhorn Constant Averaging \(DSCA\), modifies the coefficient\-calibration stage through which candidate expressions are evaluated during the evolutionary search\.

Let𝒟=\{\(xi,yi\)\}i=1N\\mathcal\{D\}=\\\{\(x\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{N\}denote the input dataset\. Standard memetic calibration fits the coefficients of each candidate once on the available training data\. DSCA replaces this single calibration step with two stages:

- •Once, before the search, the dataset𝒟\\mathcal\{D\}is used to createKKequally sized subsets with controlled covariate heterogeneity\.
- •During the search, every candidate is calibrated independently on each subset and the resulting parameter vectors are averaged to obtain the parameterization used for fitness evaluation\.

The following subsections describe these two contributions in detail\.

### III\-ASurrogate\-Based Dirichlet\-Sinkhorn Partitioning

Starting from𝒟\\mathcal\{D\}, the objective is to generateKKmutually exclusive, equally sized subsets𝒟1,…,𝒟K\\mathcal\{D\}\_\{1\},\\dots,\\mathcal\{D\}\_\{K\}\.

Our mechanism is inspired by a recent work where Dirichlet allocation is used to simulate client heterogeneity in federated SR\[[24](https://arxiv.org/html/2610.10931#bib.bib24)\]\. Here, we use a similar approach in a centralized setting to leverage controlled data heterogeneity as a probe of model calibration\. Notably, unlike\[[24](https://arxiv.org/html/2610.10931#bib.bib24)\], sample allocation is not driven by target observation quantiles, but is a function of covariates only\. As we will show in Section[III\-C](https://arxiv.org/html/2610.10931#S3.SS3), this ensures thatp⁡\(y∣x\)p\(y\\mid x\)stays identical across subsets\. This property guarantees that the population penalty vanishes for correctly specified and identifiable models, while misspecified models can incur a persistent penalty that, in the near\-homogeneous perturbative regime, scales with partition heterogeneity\.

Algorithm 1Surrogate\-Based Dirichlet\-Sinkhorn Partitioning0:Training set

𝒟\\mathcal\{D\}of size

NN, surrogate fraction

ρ\\rho, partitions

KK, strata

SS, heterogeneity

α\\alpha
1:Draw

𝒟sur⊂𝒟\\mathcal\{D\}\_\{\\text\{sur\}\}\\subset\\mathcal\{D\}with

\|𝒟sur\|=⌊ρ​N⌋\|\\mathcal\{D\}\_\{\\text\{sur\}\}\|=\\lfloor\\rho N\\rfloor; set

𝒟opt=𝒟∖𝒟sur\\mathcal\{D\}\_\{\\text\{opt\}\}=\\mathcal\{D\}\\setminus\\mathcal\{D\}\_\{\\text\{sur\}\},

N′=\|𝒟opt\|N^\{\\prime\}=\|\\mathcal\{D\}\_\{\\text\{opt\}\}\|\.

2:Fit ridge surrogate

ggon

𝒟sur\\mathcal\{D\}\_\{\\text\{sur\}\}\(penalty by LOO CV\)\.

3:Score

ui←g⁡\(xi\)u\_\{i\}\\leftarrow g\(x\_\{i\}\)for all

i∈𝒟opti\\in\\mathcal\{D\}\_\{\\text\{opt\}\}\.

4:Sort

𝒟opt\\mathcal\{D\}\_\{\\text\{opt\}\}by

uuand bin into

SSequal strata\.

5:Initialize

M∈ℝS×KM\\in\\mathbb\{R\}^\{S\\times K\}where each row

Ms,∗∼Dir​\(α​𝟏K\)M\_\{s,\*\}\\sim\\text\{Dir\}\(\\alpha\\mathbf\{1\}\_\{K\}\)\.

6:Define row targets

𝐫=\(N′/S\)​𝟏S\\mathbf\{r\}=\(N^\{\\prime\}/S\)\\mathbf\{1\}\_\{S\}and column targets

𝐜=\(N′/K\)​𝟏K\\mathbf\{c\}=\(N^\{\\prime\}/K\)\\mathbf\{1\}\_\{K\}\.

7:foriteration

1,2,…,501,2,\\dots,50do

8:Scale rows:

Ms,k←Ms,k×\(rs/∑j=1KMs,j\)M\_\{s,k\}\\leftarrow M\_\{s,k\}\\times\(r\_\{s\}/\\sum\_\{j=1\}^\{K\}M\_\{s,j\}\)
9:Scale columns:

Ms,k←Ms,k×\(ck/∑i=1SMi,k\)M\_\{s,k\}\\leftarrow M\_\{s,k\}\\times\(c\_\{k\}/\\sum\_\{i=1\}^\{S\}M\_\{i,k\}\)
10:endfor

11:Initialize

Mint=⌊M⌋M\_\{\\text\{int\}\}=\\lfloor M\\rfloor\.

12:Row deficits

dsr=rs−∑j=1KMints,jd^\{r\}\_\{s\}=r\_\{s\}\-\\sum\_\{j=1\}^\{K\}M\_\{\\text\{int\}\_\{s,j\}\}; column deficits

dkc=ck−∑i=1SMinti,kd^\{c\}\_\{k\}=c\_\{k\}\-\\sum\_\{i=1\}^\{S\}M\_\{\\text\{int\}\_\{i,k\}\}\.

13:Fractional remainders

F=M−MintF=M\-M\_\{\\text\{int\}\}\.

14:for

s=1s=1to

SSdo

15:while

dsr\>0d^\{r\}\_\{s\}\>0do

16:Select column

k∗k^\{\*\}with

dk∗c\>0d^\{c\}\_\{k^\{\*\}\}\>0maximizing

Fs,k∗F\_\{s,k^\{\*\}\}\.

17:

Mints,k∗←Mints,k∗\+1M\_\{\\text\{int\}\_\{s,k^\{\*\}\}\}\\leftarrow M\_\{\\text\{int\}\_\{s,k^\{\*\}\}\}\+1
18:

dsr←dsr−1d^\{r\}\_\{s\}\\leftarrow d^\{r\}\_\{s\}\-1;

dk∗c←dk∗c−1d^\{c\}\_\{k^\{\*\}\}\\leftarrow d^\{c\}\_\{k^\{\*\}\}\-1
19:

Fs,k∗←−1F\_\{s,k^\{\*\}\}\\leftarrow\-1\{Prevent immediate re\-selection\}

20:endwhile

21:endfor

22:Group

𝒟opt\\mathcal\{D\}\_\{\\text\{opt\}\}by strata and shuffle indices\.

23:Slice and deal to

KKpartitions according to

MintM\_\{\\text\{int\}\}\.

24:return

KKpartitions and the evaluation set

𝒟opt\\mathcal\{D\}\_\{\\text\{opt\}\}\.

Given𝒟\\mathcal\{D\}, the partitioning procedure is therefore guided by four main parameters: i\) the fractionρ\\rhoof𝒟\\mathcal\{D\}which is used to define the covariate function driving sample stratification for allocation, ii\) the number of subsetsKK, iii\) the number of strataSS, and iv\) the concentration parameter,α\\alpha, which controls data heterogeneity across partitions\. Algorithm[1](https://arxiv.org/html/2610.10931#alg1)details the procedure, which is organized into six main stages:

1. 1\.Surrogate Scoring \(steps 1\-3\):A fractionρ\\rhoof𝒟\\mathcal\{D\}is set aside as𝒟sur\\mathcal\{D\}\_\{\\text\{sur\}\}and used to fit a ridge regressiongg, predictingyyfromxx, with the penalty selected by leave\-one\-out cross\-validation\. The subset𝒟sur\\mathcal\{D\}\_\{\\text\{sur\}\}is then discarded\. The remainingN′=⌈\(1−ρ\)​N⌉N^\{\\prime\}=\\lceil\(1\-\\rho\)N\\rceilsamples form𝒟opt\\mathcal\{D\}\_\{\\text\{opt\}\}, which will be partitioned into𝒟1,…,𝒟K\\mathcal\{D\}\_\{1\},\\ldots,\\mathcal\{D\}\_\{K\}\. To this aim,xi∈𝒟optx\_\{i\}\\in\\mathcal\{D\}\_\{\\text\{opt\}\}receives the scoreui=g⁡\(xi\)u\_\{i\}=g\(x\_\{i\}\)\. The surrogate scoring functionggis a deterministic function of the covariates alone, and fittingggon held\-out data keeps it independent of the targets later used for calibration\. Surrogate quality affects how informative the induced heterogeneity is, not the validity of the argument \(see Section[III\-C](https://arxiv.org/html/2610.10931#S3.SS3)\): a poor surrogate simply yields nearly random, weakly heterogeneous partitions\. Any covariate\-measurable score can be used, including a single input variable correlated with the target or the output of a shallow tree ensemble\.
2. 2\.Stratification \(step 4\):𝒟opt\\mathcal\{D\}\_\{\\text\{opt\}\}is sorted byuuand divided intoSSequal\-sized strata ofN′/SN^\{\\prime\}/Ssamples each\.
3. 3\.Dirichlet Sampling \(step 5\):AnS×KS\\times Kallocation matrixMMis initialized to define the allocation of each stratum across partitions\. Each row is drawn from a symmetric Dirichlet distribution,Ms,∗∼Dirichlet​\(α​𝟏K\)M\_\{s,\*\}\\sim\\text\{Dirichlet\}\(\\alpha\\mathbf\{1\}\_\{K\}\), so thatMs,kM\_\{s,k\}is the fraction of stratumsstentatively assigned to partitionkk, and each row describes how stratumssis spread over theKKpartitions\. Largeα\\alphaspreads each stratum more evenly across partitions, so that the𝒟k\\mathcal\{D\}\_\{k\}become increasingly similar in covariate composition to𝒟opt\\mathcal\{D\}\_\{\\mathrm\{opt\}\}; smallα\\alphaconcentrates each stratum in fewer partitions and produces stronger covariate heterogeneity\.
4. 4\.Matrix Balancing \(steps 6\-10\):MMis subjected to the Sinkhorn\-Knopp algorithm\[[37](https://arxiv.org/html/2610.10931#bib.bib37)\]to enforce the marginal constraints: RowsMs,∗M\_\{s,\*\}must sum toN′/SN^\{\\prime\}/Sso that every sample is used, and columnsM∗,kM\_\{\*,k\}toN′/KN^\{\\prime\}/Kso that partitions are equally sized\. The matrix is updated by alternating row and column scaling: Ms,k←Ms,k​N′/S∑j=1KMs,jMs,k←Ms,k​N′/K∑i=1SMi,kM\_\{s,k\}\\leftarrow M\_\{s,k\}\\frac\{N^\{\\prime\}/S\}\{\\sum\_\{j=1\}^\{K\}M\_\{s,j\}\}\\quad M\_\{s,k\}\\leftarrow M\_\{s,k\}\\frac\{N^\{\\prime\}/K\}\{\\sum\_\{i=1\}^\{S\}M\_\{i,k\}\}\(1\)This iterative proportional fitting converges to a balanced matrix of fractional allocations, typically in fewer than ten iterations\.
5. 5\.Integer Rounding \(steps 11\-21\):The balanced matrix is converted to integers byMint=⌊M⌋M\_\{\\text\{int\}\}=\\lfloor M\\rfloor, and the residual samples are greedily assigned by incrementing the cells with the largest fractional remainders until the integer margins match the targets exactly\.
6. 6\.Slicing \(steps 22\-23\):Sample indices are grouped by stratum, shuffled, and sliced mutually exclusively according toMintM\_\{\\text\{int\}\}\.

Since partitioning is performed once, before the evolutionary search begins, and the surrogate fit is a single ridge regression on at mostρ​N\\rho Npoints, the overhead is negligible relative to the search itself\.

### III\-BLocal Optimization and Averaging

Before numerical calibration, each candidate expressionf⁡\(⋅,θ\)f\(\\cdot,\\theta\)is simplified symbolically\. This preprocessing removes exact algebraic redundancies and structurally non\-identifiable parameterizations whenever they can be detected, thereby reducing trivial non\-identifiabilities before parameter vectors are optimized and averaged\. For example, a parameterization such asθ1​θ2​eθ3​h​\(x\)\\theta\_\{1\}\\theta\_\{2\}e^\{\\theta\_\{3\}\}h\(x\)depends only on a single effective multiplicative coefficient and should therefore be collapsed before numerical optimization to, e\.g\.,θ1​h​\(x\)\\theta\_\{1\}h\(x\)\.

During the evolutionary search,f⁡\(⋅,θ\)f\(\\cdot,\\theta\)is calibrated independently on each of theKKgenerated partitions\. To this end, we apply a gradient\-based optimizer, here BFGS\[[21](https://arxiv.org/html/2610.10931#bib.bib21)\], to each individual partition, yieldingKKlocally calibrated parameter vectorsθ^1,…,θ^K\\hat\{\\theta\}\_\{1\},\\ldots,\\hat\{\\theta\}\_\{K\}, whereθ^k\\hat\{\\theta\}\_\{k\}denotes the numerical solution returned on partitionkk\.

TheKKsets of local constants are aggregated by computing their arithmetic mean\. This produces a single, averaged set of constants for the candidate expression:

θ¯=1K​∑k=1Kθ^k\.\\bar\{\\theta\}=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\hat\{\\theta\}\_\{k\}\\,\.\(2\)Finally, the fitness of the candidate is computed on𝒟opt\\mathcal\{D\}\_\{\\text\{opt\}\}, the union of theKKpartitions, fixing parameter calibration toθ¯\\bar\{\\theta\}\.

### III\-CTheoretical analysis of DSCA implications

Replacing centralized calibration with the average ofKKlocally calibrated parameter vectors changes the loss assigned to a candidate expression\. We analyze this effect in five steps\. First, we isolate the excess loss that is directly relevant to evolutionary selection\. We then formalize the partitioned data as covariate\-reweighted population environments and characterize this excess under correct specification and misspecification\. Next, we connect the population result with the empirical finite\-sample behavior\. We then study how the surrogate scoring functionggand the heterogeneity parameterα\\alphadetermine whether, and how strongly, the mechanism acts\. Finally, we discuss how the assumptions of the theoretical analysis relate to the numerical implementation\. Proofs and full derivations are given in the supplementary material\.

Throughout this section, we use the following notation\. Let\(X,Y\)\(X,Y\)denote an input\-target variables andf⁡\(X,θ\)f\(X,\\theta\)a candidate symbolic expression with parameter vectorθ\\theta\. We interpretf⁡\(X,θ\)f\(X,\\theta\)as a model for the conditional mean ofYYgivenXX\. Equivalently, under additive Gaussian noise,Y=f⁡\(X,θ\)\+εY=f\(X,\\theta\)\+\\varepsilon, whereε∼𝒩⁡\(0,σ2\)\\varepsilon\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}\), thus the candidate induces the conditional modelpθ​\(y∣x\)=𝒩⁡\(f⁡\(x,θ\),σ2\)p\_\{\\theta\}\(y\\mid x\)=\\mathcal\{N\}\(f\(x,\\theta\),\\sigma^\{2\}\)\. For fixed noise variance, maximizing this likelihood with respect toθ\\thetais equivalent to minimizing mean squared error\. We therefore describe coefficient calibration through its empirical loss:Lk​\(θ\)L\_\{k\}\(\\theta\)denotes the mean loss on partition𝒟k\\mathcal\{D\}\_\{k\}, whileLN′​\(θ\)L\_\{N^\{\\prime\}\}\(\\theta\)denotes the corresponding loss on the complete optimization set𝒟opt\\mathcal\{D\}\_\{\\mathrm\{opt\}\}\. Their population counterparts, introduced below, are denoted byΛk​\(θ\)\\Lambda\_\{k\}\(\\theta\)andΛ⁡\(θ\)\\Lambda\(\\theta\)\.

#### III\-C1DSCA and selection\-relevant excess loss

Because the partitions are disjoint and equally sized,

LN′​\(θ\)=1K​∑k=1KLk​\(θ\)L\_\{N^\{\\prime\}\}\(\\theta\)=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}L\_\{k\}\(\\theta\)\(3\)holds exactly for everyθ\\theta\. Letθ^k=arg⁡minθ​Lk​\(θ\)\\hat\{\\theta\}\_\{k\}=\\arg\\min\_\{\\theta\}L\_\{k\}\(\\theta\)be the local empirical optima,θ¯=K−1​∑kθ^k\\bar\{\\theta\}=K^\{\-1\}\\sum\_\{k\}\\hat\{\\theta\}\_\{k\}their average, andθ^N′=arg⁡minθ​LN′​\(θ\)\\hat\{\\theta\}\_\{N^\{\\prime\}\}=\\arg\\min\_\{\\theta\}L\_\{N^\{\\prime\}\}\(\\theta\)the centralized optimum\.

DSCA evaluates the candidate atθ¯\\bar\{\\theta\}, whereas centralized calibration would evaluate the same structure atθ^N′\\hat\{\\theta\}\_\{N^\{\\prime\}\}\. The selection\-relevant excess loss introduced by DSCA is therefore

ΔN′:=LN′​\(θ¯\)−LN′​\(θ^N′\)\.\\Delta\_\{N^\{\\prime\}\}:=L\_\{N^\{\\prime\}\}\(\\bar\{\\theta\}\)\-L\_\{N^\{\\prime\}\}\(\\hat\{\\theta\}\_\{N^\{\\prime\}\}\)\.\(4\)By construction,ΔN′≥0\\Delta\_\{N^\{\\prime\}\}\\geq 0under exact minimization\. It measures how much additional empirical loss is assigned to a fixed candidate structure because DSCA uses the average of the locally calibrated parameters rather than the centralized optimum\. This is the quantity through which DSCA can modify evolutionary selection\.

#### III\-C2Covariate\-induced environments and population behavior

DSCA operates empirically on theKKfinite partitions𝒟1,…,𝒟K\\mathcal\{D\}\_\{1\},\\ldots,\\mathcal\{D\}\_\{K\}\. To determine whether the resulting excessΔN′\\Delta\_\{N^\{\\prime\}\}is purely a finite\-sample effect or persists at population level, we associate each finite partition𝒟k\\mathcal\{D\}\_\{k\}with a population environmentPkP\_\{k\}induced by the same covariate\-reweighting mechanism\.

LetPPdenote the original population law of\(X,Y\)\(X,Y\), with densityp⁡\(x,y\)=p⁡\(y∣x\)​p​\(x\)p\(x,y\)=p\(y\\mid x\)p\(x\)\. Given the fitted surrogategg, letA∈\{1,…,K\}A\\in\\\{1,\\ldots,K\\\}denote the population environment assigned to an observation and defineπk​\(u\):=Pr⁡\(A=k∣g⁡\(x\)=u\)\\pi\_\{k\}\(u\):=\\Pr\(A=k\\mid g\(x\)=u\)\. Since assignment depends on the observation only throughg⁡\(X\)g\(X\),Pr⁡\(A=k∣X=x,Y=y\)=πk​\(g⁡\(x\)\)\\Pr\(A=k\\mid X=x,Y=y\)=\\pi\_\{k\}\(g\(x\)\)\. Equal\-size balancing givesPr⁡\(A=k\)=1/K\\Pr\(A=k\)=1/K, and hence environmentkkhas densitypk​\(x,y\)=p⁡\(x,y∣A=k\)=K​πk​\(g⁡\(x\)\)​p​\(x,y\)p\_\{k\}\(x,y\)=p\(x,y\\mid A=k\)=K\\,\\pi\_\{k\}\(g\(x\)\)\\,p\(x,y\)\. Consequently,K−1​∑kpk​\(x,y\)=p⁡\(x,y\)K^\{\-1\}\\sum\_\{k\}p\_\{k\}\(x,y\)=p\(x,y\)and, more importantly,

pk​\(y∣x\)=p⁡\(y∣x\)for every​k\.p\_\{k\}\(y\\mid x\)=p\(y\\mid x\)\\qquad\\text\{for every \}k\\,\.\(5\)Thus the environments may differ in their covariate marginals while preserving the same conditional lawp⁡\(y∣x\)p\(y\\mid x\)\. This requires only that, once fitted, the allocation score be a function ofxx; the quality of the surrogate allocation functionggaffects the informativeness of the induced heterogeneity, not the validity of \([5](https://arxiv.org/html/2610.10931#S3.E5)\)\.

The population counterparts ofLkL\_\{k\}andLN′L\_\{N^\{\\prime\}\}areΛk​\(θ\):=𝔼Pk​\[\(Y−f⁡\(X,θ\)\)2\],\\Lambda\_\{k\}\(\\theta\):=\\mathbb\{E\}\_\{P\_\{k\}\}\\\!\\left\[\(Y\-f\(X,\\theta\)\)^\{2\}\\right\],andΛ⁡\(θ\):=𝔼P​\[\(Y−f⁡\(X,θ\)\)2\]\\Lambda\(\\theta\):=\\mathbb\{E\}\_\{P\}\\\!\\left\[\(Y\-f\(X,\\theta\)\)^\{2\}\\right\], with minimizersθk∗\\theta\_\{k\}^\{\*\}andθ∗\\theta^\{\*\}, respectively\. Under the environment construction conditions, the population counterpart ofΔN′\\Delta\_\{N^\{\\prime\}\}is thenΔ:=Λ⁡\(θ¯∗\)−Λ⁡\(θ∗\)\\Delta:=\\Lambda\(\\bar\{\\theta\}^\{\*\}\)\-\\Lambda\(\\theta^\{\*\}\), whereθ¯∗=K−1​∑kθk∗\\bar\{\\theta\}^\{\*\}=K^\{\-1\}\\sum\_\{k\}\\theta\_\{k\}^\{\*\}\. We now inspect, under identifiability assumptions, howΔ\\Deltachanges between well specified and misspecified models\.

Correctly specified models\.A candidate is correctly specified for squared\-error calibration if there exists a parameter vectorθ∗\\theta^\{\*\}such thatf⁡\(X,θ∗\)=𝔼⁡\[Y∣X\]f\(X,\\theta^\{\*\}\)=\\mathbb\{E\}\[Y\\mid X\]almost surely\. Because the induced environments modify only the covariate marginal while preservingp⁡\(y\|x\)p\(y\|x\), the same optimum is retained in every environment\. Under regularity and support conditions \(details can be found in the supplementary material\),θk∗=θ∗\\theta\_\{k\}^\{\*\}=\\theta^\{\*\}for allkkandθ¯∗=θ∗\\bar\{\\theta\}^\{\*\}=\\theta^\{\*\}, leading toΔ=0\\Delta=0\. This is the main reason why we define allocation based onxxonly: target\-dependent partitioning can alterp⁡\(y∣x\)p\(y\\mid x\)and therefore move the optimum even for a correctly specified candidate\.

Misspecified models\.A candidate is misspecified when no parameter vectorθ\\thetareproduces the conditional mean𝔼⁡\[Y∣X\]\\mathbb\{E\}\[Y\\mid X\]\. Its environment\-specific optimumθk∗\\theta\_\{k\}^\{\*\}is then a covariate\-weighted pseudo\-true parameter\. Reweighting different regions of covariate space can change the approximation trade\-off made by the candidate and therefore produce different preferred parameter values across environments\. Under a local quadratic approximation, the supplementary material shows that

Δ=12​‖1K​∑k=1K\(Hk−H¯\)​\(θk∗−θ¯∗\)‖H¯−12,\\Delta=\\frac\{1\}\{2\}\\left\\\|\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\(H\_\{k\}\-\\bar\{H\}\)\(\\theta\_\{k\}^\{\*\}\-\\bar\{\\theta\}^\{\*\}\)\\right\\\|\_\{\\bar\{H\}^\{\-1\}\}^\{2\},\(6\)whereHk=∇2Λk​\(θk∗\)H\_\{k\}=\\nabla^\{2\}\\Lambda\_\{k\}\(\\theta\_\{k\}^\{\*\}\)is the local curvature in environmentkk,H¯=K−1​∑kHk\\bar\{H\}=K^\{\-1\}\\sum\_\{k\}H\_\{k\}is the average curvature, and‖v‖H¯−12=v⊤​H¯−1​v\\\|v\\\|\_\{\\bar\{H\}^\{\-1\}\}^\{2\}=v^\{\\top\}\\bar\{H\}^\{\-1\}v\. ThusΔ\\Deltadepends jointly on how the preferred parameters shift across environments and on how their local curvatures differ\. This shows that misspecification alone is not sufficient to always ensure a population\-level DSCA penalty\. If allHkH\_\{k\}coincide or if by chance the exact cancellation∑k=1K\(Hk−H¯\)​dk=0\\sum\_\{k=1\}^\{K\}\(H\_\{k\}\-\\bar\{H\}\)d\_\{k\}=0holds, thenΔ=0\\Delta=0even when theθk∗\\theta\_\{k\}^\{\*\}differ substantially\. Absent such a symmetry or compensating structure, heterogeneous curvatures and displacements generically produce a nonzero imbalance, so thatθ¯∗≠θ∗\\bar\{\\theta\}^\{\*\}\\neq\\theta^\{\*\}and a persistent population excess arises\.

Overall, the population effect characterized above is governed by specification, identifiability, and sensitivity of the preferred calibration to the induced environments, rather than by parameter count itself\. DSCA is therefore not an explicit complexity penalty: candidates with many coefficients need not incur a population penalty, while misspecified candidates can incur one when heterogeneous calibration moves the averaged parameters away from the pooled optimum\. Structural complexity must consequently remain controlled by the outer selection mechanism\.

#### III\-C3Finite\-sample behavior

The preceding results describe the population effect, whereas DSCA operates on finite empirical partitions\. We therefore ask whether the observed excessΔN′\\Delta\_\{N^\{\\prime\}\}is only a sampling discrepancy or instead converges to the population quantityΔ\\Delta\. Under standard regularity conditions and assuming identifiable models, fixedKK, and largeN′N^\{\\prime\},

ΔN′=\{Op​\(1/N′\),correctly specified,Δ\+Op​\(1/N′\),misspecified;\\Delta\_\{N^\{\\prime\}\}=\\begin\{cases\}O\_\{p\}\(1/N^\{\\prime\}\),&\\text\{correctly specified\},\\\\\[2\.0pt\] \\Delta\+O\_\{p\}\(1/\\sqrt\{N^\{\\prime\}\}\),&\\text\{misspecified \};\\end\{cases\}\(7\)whereOp​\(⋅\)O\_\{p\}\(\\cdot\)denotes stochastic order, i\.e\.,XN′=Op​\(aN′\)X\_\{N^\{\\prime\}\}=O\_\{p\}\(a\_\{N^\{\\prime\}\}\)means thatXN′/aN′X\_\{N^\{\\prime\}\}/a\_\{N^\{\\prime\}\}remains bounded in probability\. Thus, under correct specification the entire DSCA excess is a second\-order finite\-sample effect, whereas whenΔ\>0\\Delta\>0a nonzero population excess persists and only its sampling fluctuations vanish, at a slower rate\.

#### III\-C4The role of environment definition and heterogeneity

The DSCA effect depends on two complementary aspects of environment construction: the surrogate allocation functionggdetermines which regions of covariate space are contrasted, whereasα\\alphacontrols the strength of the resulting heterogeneity\.

Surrogate informativeness\.For a misspecified candidate, an informative surrogate should induce environments that favor different calibrations\. Under squared\-error loss, the relevant interaction for environmentkkis

rk:=𝔼P​\[δk​\(X\)​∇θ𝔼​\[\(Y−f⁡\(X,θ\)\)2∣X\]θ=θ∗\],r\_\{k\}:=\\mathbb\{E\}\_\{P\}\\\!\\left\[\\delta\_\{k\}\(X\)\\nabla\_\{\\theta\}\\mathbb\{E\}\\\!\\left\[\(Y\-f\(X,\\theta\)\)^\{2\}\\mid X\\right\]\_\{\\theta=\\theta^\{\*\}\}\\right\],\(8\)whereδk​\(X\)=πk​\(g⁡\(X\)\)−1/K\\delta\_\{k\}\(X\)=\\pi\_\{k\}\(g\(X\)\)\-1/Kis the covariate reweighting induced by the surrogate\. Thusrkr\_\{k\}measures whether that reweighting is aligned with regions in which the candidate prefers a different calibration: smallrkr\_\{k\}implies little first\-order displacement of the local optimum, whereas nonzerorkr\_\{k\}induces a calibration shift\. The supplementary material formalizes this relation perturbatively\.

Choosingggas a predictor ofYYprovides a candidate\-independent, response\-informed ordering of covariate space, making it more likely to contrast regions with different approximation trade\-offs/estimated noise levels\. It is not optimal for every candidate, since the relevant response also depends on candidate\-specific parameter sensitivity, but it provides a common response\-informed partitioning rule that can be constructed once before the search\.

Heterogeneity level\.Within DSCA,α\\alphais the concentration parameter of the initial Dirichlet allocation and therefore controls the strength of the induced covariate heterogeneity\. Asα→∞\\alpha\\to\\infty, the allocation probabilities approach1/K1/K, the environments become homogeneous,Pk→PP\_\{k\}\\to P, andΔ\\Deltavanishes\. More precisely, the perturbation analysis in the supplementary material gives

Δ=O⁡\(\(K​α\+1\)−2\)as​α→∞\.\\Delta=O\\\!\\left\(\(K\\alpha\+1\)^\{\-2\}\\right\)\\qquad\\text\{as \}\\alpha\\to\\infty\.Hence nonzero heterogeneity is necessary for DSCA to induce a population\-level selection effect, although the result does not imply thatΔ\\Deltais monotone inα\\alphaor identify an optimal heterogeneity level\.

Heterogeneity also cannot be increased indefinitely, since excessively concentrated environments may provide insufficient covariate support for local identification\. Parameter variation can then reflect weak identification or induced near\-singularity rather than genuine sensitivity to covariate reweighting\.

From theoretical assumptions to practical implementation\.The preceding analysis is stated in terms of exact empirical minimizers and identifiable population optima\. In the implementation, parameter estimates are obtained numerically with BFGS and therefore need not attain the corresponding global empirical minima\. As a diagnostic of obvious optimization failure, we checked the ordering implied by local minimization,

Lk​\(θ¯\)≥Lk​\(θ^k\),L\_\{k\}\(\\bar\{\\theta\}\)\\geq L\_\{k\}\(\\hat\{\\theta\}\_\{k\}\),\(9\)which holds for 98\.9% of the evaluated candidate\-partition pairs\. Violations are therefore rare, indicating that the observed DSCA averaging effect is not driven by obvious failures of the local optimizer\.

Candidate expressions are also symbolically simplified before calibration\. This preprocessing removes exact algebraic redundancies whenever they can be identified symbolically, thereby strengthening the correspondence with the identifiability assumption used above and improving the conditioning of the subsequent calibration step\. Residual near\-singular or weakly identified parameter directions may still arise from the local statistical geometry of the candidate expression\.

## IVExperimental Setup

### IV\-AGenetic Programming Framework

For our experimental evaluations, we employ a standard tree\-based Genetic Programming algorithm as the foundational framework\. Because the contribution of this work concerns the coefficient\-calibration mechanism used during individual evaluation rather than the evolutionary operators themselves, we keep the GP configuration fixed across all compared methods\. The complete set of evolutionary hyperparameters is detailed in Table[I](https://arxiv.org/html/2610.10931#S4.T1)and explained in the following\.

The initial population is generated using a random initialization method with decay\. In this approach, tree depth is regulated by two hyperparameters: parsimony \(representing the initial operator probability at the root node\) and a parsimony decay \(which reduces this probability as tree depth increases\)\. During the evolutionary cycle, parent selection is performed via tournament selection\. The population is subsequently evolved using a uniformly distributed set of genetic operators to maintain structural diversity\.

The critical distinguishing feature of our implementation is the integration of the proposed coefficient calibration, which transforms the standard GP into a memetic algorithm\. After variation generates a new individual, its expression is symbolically simplified and its free numerical parameters are calibrated before fitness evaluation\. The calibration rule is the only component changed across the compared methods\. The memetic phase proceeds as follows for each newly generated individual:

1. 1\.Parameter Extraction:After symbolic simplification, the tree is traversed to identify and extract all free constants and parameterizable operator coefficients\.
2. 2\.Local Search:The selected calibration procedure is applied to the extracted parameters\. Centralized baselines optimize them on the complete optimization set, whereas DSCA performs the partition\-wise optimization and averaging procedure described in Section[III](https://arxiv.org/html/2610.10931#S3)\.
3. 3\.Parameter Injection:Following a Lamarckian evolutionary scheme, the resulting calibrated parameters are directly injected back into the individual’s symbolic tree, permanently updating its genetic material\.

Only after this inner calibration step terminates, the next generation of individuals is selected\. Specifically, each calibrated individual is scored according to one of the three selection criteria and the survival strategy is applied to select the best\-scoring individuals from the combined parent\-offspring pool\.

Coefficient calibration therefore determines the parameter values at which an individual is evaluated, whereas the outer criterion determines how calibrated individuals compete during evolution\. To test whether the effect of DSCA depends on a particular outer selection rule, we evaluate each calibration strategy under three criteria commonly used in SR: the negative log\-likelihood \(NLL\), AIC, and BIC\[[11](https://arxiv.org/html/2610.10931#bib.bib11),[38](https://arxiv.org/html/2610.10931#bib.bib38),[28](https://arxiv.org/html/2610.10931#bib.bib28)\]\. These criteria span unpenalized likelihood fitting and two standard levels of parameter\-count penalization,2​k2kandk​ln⁡Nk\\ln N, wherekkis the number of free parameters andNNthe number of training samples\. Recent comparative studies indicate that they form a compact and competitive set, since more elaborate alternatives are not uniformly superior when used directly to guide evolutionary search\[[11](https://arxiv.org/html/2610.10931#bib.bib11),[38](https://arxiv.org/html/2610.10931#bib.bib38),[28](https://arxiv.org/html/2610.10931#bib.bib28)\]\. The purpose of this comparison is thus not to benchmark model\-selection criteria, but to assess whether the effect of DSCA is robust across them\.

TABLE I:Evolutionary and numerical optimization hyperparameters\.ParameterValue/SettingPopulation size \(μ\\mu\)100100Offspring size \(λ\\lambda\)100100Number of generations100100Initialization methodRandom with decayParsimony0\.80\.8Parsimony decay0\.850\.85Maximum tree size7575nodesFunction set\{\+,−,×,÷,\(⋅\)2,\(⋅\)3,\\\{\+,\-,\\times,\\div,\(\\cdot\)^\{2\},\(\\cdot\)^\{3\},sin,exp,ln\}\\sin,\\exp,\\ln\\\}Terminal setIndependent variablesX1,…,XnX\_\{1\},\\dots,X\_\{n\}, ERCs∈\[−10,10\]\\in\[\-10,10\]Selection methodTournament selection \(sizek=3k=3\)Genetic OperatorsCrossover, Subtree mutation, Insert node, Delete node, Mutate leaf, Mutate operatorSelection criteria\[NLL, AIC, BIC\]Survival strategy\(μ\+λ\)\(\\mu\+\\lambda\)selection \(best100100survive\)Numerical Optimizer SettingsOptimization frequencyApplied to all newly generated individualsProposed Method: DSCABase optimizerBFGS \(settings as in Baseline 1\)Number of partitions \(KK\)55Number of strata \(SS\)1010Heterogeneity \(α\\alpha\)11Stratification scoreRidge predictiong⁡\(x\)g\(x\)Surrogate fraction \(ρ\\rho\)15%15\\%of the training setRidge penaltyLeave\-one\-out CV over\{10−3,…,103\}\\\{10^\{\-3\},\\dots,10^\{3\}\\\}Baseline 1: BFGS \(via SciPy ‘minimize’\)Jacobian estimation2\-point finite differenceTolerance \(tol\)None \(relies on defaults\)Baseline 2: Levenberg\-Marquardt \(via SciPy ‘least\_squares’\)Jacobian estimation2\-point finite differenceCost change \(ftol\)1×10−81\\times 10^\{\-8\}Step change \(xtol\)1×10−81\\times 10^\{\-8\}Gradient norm \(gtol\)1×10−81\\times 10^\{\-8\}
### IV\-BBaselines

We compare against two established centralized calibration methods: Levenberg\-Marquardt \(LM\), the standard local search in SR frameworks such as Operon\[[20](https://arxiv.org/html/2610.10931#bib.bib20)\], and BFGS, the default optimizer in PySR\[[8](https://arxiv.org/html/2610.10931#bib.bib8)\]\. Both optimize the same within\-individual objective as the local fits of DSCA, so differences in evolutionary behavior can be attributed to the calibration strategy rather than to a change in loss function\. Throughout, prediction errors are assumed independent and normally distributed, under which minimizing the NLL with respect toθ\\thetais equivalent to least\-squares fitting\.

To ensure a fair comparison, all baseline optimizers are integrated into the memetic framework in exactly the same manner as our proposed approach\. Table[I](https://arxiv.org/html/2610.10931#S4.T1)details the hyperparameters used for these algorithms, which follow standard default configurations\. In addition, our ablation experiments include a multi\-start centralized BFGS control, in which the same structure is optimized repeatedly on the pooled data\. This will help distinguish the effect of heterogeneous averaging from the simpler effect of allocating multiple optimization attempts to each candidate\.

One asymmetry deserves to be stated explicitly\. DSCA reserves15%15\\%of the training samples to fit the stratification surrogate and never uses them again, so it optimizes constants and evaluates fitness onN′=0\.85​NN^\{\\prime\}=0\.85Nsamples, whereas the two baselines use allNN\. We kept the nominal training size fixed across methods rather than equalizingN′N^\{\\prime\}, so the comparison is conservative with respect to our own method\.

### IV\-CSynthetic Track

The first track evaluates model recall on synthetic datasets where the true generating equations are known, allowing us to explicitly control the injected noise\. Table[II](https://arxiv.org/html/2610.10931#S4.T2)shows the generating functions tested\. We selected these from widely used benchmarks to ensure a diverse representation of functional forms and operators, including polynomials, rational expressions, exponentials, and trigonometric functions\.

We simulate two noise regimes by adopting the methodology from\[[3](https://arxiv.org/html/2610.10931#bib.bib3)\]\. Specifically, we add Gaussian noise to the target values:

yn​o​i​s​e∼𝒩⁡\(y,σy​r​a​t​i​o1−r​a​t​i​o\)y\_\{noise\}\\sim\\mathcal\{N\}\\left\(y,\\sigma\_\{y\}\\sqrt\{\\frac\{ratio\}\{1\-ratio\}\}\\right\)\(10\)whereyyis the true target value,σy\\sigma\_\{y\}is the standard deviation of the target, andr​a​t​i​oratiodetermines the noise level\. For our experiments, we use noise ratios of 0\.025 and 0\.1, corresponding to the*easy*and*hard*levels in\[[3](https://arxiv.org/html/2610.10931#bib.bib3)\], and both 100 and 200 samples for the training set\.

We consider an equation successfully recovered if it achieves anR2\>0\.99R^\{2\}\>0\.99on the noiseless test set\. Some studies verify exact symbolic equivalence by checking if the ratio or difference relative to the generating function is constant \(typically using libraries like SymPy\)\. We initially tested this methodology but observed too many false negatives, as computer algebra systems frequently fail to simplify structurally equivalent expressions to a common syntactic form\.

TABLE II:Generating functions and variable distributions for the synthetic benchmark track\.
### IV\-DReal\-World Track

For the second track, we selected 10 real\-world datasets from the PMLB collection\[[42](https://arxiv.org/html/2610.10931#bib.bib42)\], which serves as a standard benchmark within the SR community\[[43](https://arxiv.org/html/2610.10931#bib.bib43),[3](https://arxiv.org/html/2610.10931#bib.bib3)\]\. We chose these specific datasets to ensure a diverse representation of sample sizes, dimensionalities, ratios of binary to continuous variables, and overall difficulty levels for SR algorithms\. Table[III](https://arxiv.org/html/2610.10931#S4.T3)details the characteristics of each dataset, including the total number of instances and the test set sizes\.

To ensure statistical robustness, we conduct 10 independent runs for each dataset and sample size regime\. To prevent data leakage and simulate these different regimes, we employ a splitting strategy for each run\. First, we reserve 25% of the data, sampled uniformly at random, to serve as the test set\. From the remaining 75%, we again sample uniformly at random to extract the necessary number of instances to form the training set\. Specifically, we evaluate three training set size regimes: 100, 200, and 500 samples\.

TABLE III:Characteristics of the selected PMLB datasets\.

## VResults

### V\-ASynthetic Track

Table[IV](https://arxiv.org/html/2610.10931#S5.T4)summarizes the recovery rates on the synthetic benchmarks\. For every combination of sample size, noise level, and selection criterion, it reports on how many of the five generating functions DSCA attains a higher \(win\), equal \(tie\), or lower \(loss\) success rate than each baseline over 10 seeds\. DSCA wins in 64 of the 120 function\-level comparisons, ties in 48, and loses in 8\. It never loses on more than one function in any configuration, and it never loses to BFGS in the 200\-sample regime\.

Per\-function success rates are reported in the supplementary material\. Ties are largely driven byf3f\_\{3\}, on which all methods fail almost systematically\. All losses occur on the rational functionsf4f\_\{4\}andf5f\_\{5\}, where the three optimizers perform comparably and differences are typically within one or two seeds\. The wins are concentrated on the polynomial functionsf1f\_\{1\}andf2f\_\{2\}, and this holds in both noise regimes\.

It is worth noting where the advantage is largest\. Onf1f\_\{1\}andf2f\_\{2\}, which have well\-conditioned parameterizations, DSCA gains most under the unregularized NLL criterion, where the baselines have no defense against overparameterized alternatives\. Under BIC the gap narrows, since the two mechanisms then act on partly overlapping failure modes\. This pattern recurs, more clearly, on the real\-world track\.

Table[V](https://arxiv.org/html/2610.10931#S5.T5)presents the aggregated success counts across the different synthetic datasets\. These results demonstrate the superior performance of our proposed approach compared to the baselines, maintaining this advantage both overall and within each individual model selection criterion\.

TABLE IV:Win/tie/loss counts of DSCA against each baseline over the synthetic datasets, by sample size, noise level, and selection criterion\.TABLE V:Aggregate success counts on synthetic datasets by sample size and selection criterion\.
### V\-BReal\-World Track

#### V\-B1Bayesian Statistical Analysis

Before comparing the algorithms, we aggregate each configuration by taking the median over the 10 independent seeds of every dataset, and we clipR2R^\{2\}values below−1\-1so that degenerate fits do not dominate the comparison\. The same aggregation is used for the analysis reported in the remainder of this section\.

We compare DSCA against BFGS and LM using the Bayesian Wilcoxon signed\-rank test\[[44](https://arxiv.org/html/2610.10931#bib.bib44),[45](https://arxiv.org/html/2610.10931#bib.bib45)\]\. The test provides posterior probabilities that DSCA is better, practically equivalent, or worse than each baseline\. We perform the comparison separately for predictive accuracy, measured byR2R^\{2\}, and model complexity, measured by the number of expression\-tree nodes\. Practical equivalence is defined through Regions of Practical Equivalence \(ROPEs\):\[−0\.01,0\.01\]\[\-0\.01,0\.01\]for differences inR2R^\{2\}and\[−2,2\]\[\-2,2\]nodes for differences in complexity\. Differences within these intervals are treated as practically negligible\.

Figure[1](https://arxiv.org/html/2610.10931#S5.F1)reports the posterior comparison for predictive accuracy and model complexity\. In the 100\- and 200\-sample regimes, DSCA generally has a high posterior probability of outperforming both centralized baselines\. As the sample size increases, this advantage narrows\. This convergence in predictive accuracy should not be interpreted as a consequence of a vanishing population DSCA effect: Eq\. \([7](https://arxiv.org/html/2610.10931#S3.E7)\) shows that, for misspecified candidates withΔ\>0\\Delta\>0, increasingN′N^\{\\prime\}only reduces the sampling fluctuations around the persistent excessΔ\\Delta\. Rather, with more data the search has less room to retain structures whose apparent advantage depends mainly on finite\-sample calibration effects\. The convergence is particularly visible under BIC, whose strong complexity pressure and M\-closed interpretation can substantially restrict the set of competitive structures\. AtN=500N=500under BIC, DSCA also closely tracks its centralized BFGS counterpart in predictive accuracy, suggesting that residual differences may partly reflect the shared base optimizer rather than the averaging mechanism\. Importantly, in Figure[1](https://arxiv.org/html/2610.10931#S5.F1)the same convergence is not observed for complexity: DSCA retains a high posterior probability of producing smaller expressions across configurations\. This divergence suggests a non\-trivial interplay between heterogeneous calibration and explicit complexity control: stronger outer penalization reduces the predictive differences between methods without eliminating the structural effect induced by DSCA\.

Fig\. 1:Bayesian Wilcoxon signed\-rank test results evaluating predictive accuracy \(top panels\) and model complexity \(bottom panels\) against the baselines\.
#### V\-B2Absolute Performance Differences

While rank\-based evaluations effectively illustrate relative positioning, they inherently hide the absolute magnitudes of performance gaps between algorithms\. To provide a concrete measure of these performance variations, Tables[VI](https://arxiv.org/html/2610.10931#S5.T6)and[VII](https://arxiv.org/html/2610.10931#S5.T7)detail the absolute differences between the evaluated algorithms and a standard baseline\. We selected LM as the baseline, given its widespread adoption within the SR community\.

To ensure results are comparable across datasets that yield different baseline accuracies and equation sizes, we aggregate the scores by first calculating the median performance across seeds for both the evaluated algorithm and the baseline, and then computing the difference\. For a given algorithmaa, evaluated on datasetdd, we define the aggregated differenceΔa\\Delta\_\{a\}as:

Δa=1ND∑d=1ND\\displaystyle\\Delta\_\{a\}=\\frac\{1\}\{N\_\{D\}\}\\sum\_\{d=1\}^\{N\_\{D\}\}\(medi=110​\(Scorea,d,i\)−medi=110​\(ScoreLM,d,i\)\)\\displaystyle\\Big\(\\text\{med\}\_\{i=1\}^\{10\}\(\\text\{Score\}\_\{a,d,i\}\)\-\\text\{med\}\_\{i=1\}^\{10\}\(\\text\{Score\}\_\{\\text\{LM\},d,i\}\)\\Big\)\(11\)whereNDN\_\{D\}is the total number of datasets\.

TABLE VI:Predictive accuracy \(R2R^\{2\}\) differences \(Δa\\Delta\_\{a\}\) relative to the Levenberg\-Marquardt baseline \(Mean±\\pmStd\)\.TABLE VII:Model complexity differences \(Δa\\Delta\_\{a\}\) relative to the Levenberg\-Marquardt baseline \(Mean±\\pmStd\)\.BFGS returns solutions that are remarkably similar to LM: its accuracy gains are limited to the unregularized and weakly penalized criteria in the smaller samples, and the expressions it retains are essentially of the same size as the baseline throughout\. DSCA instead improves on LM along both axes at once\. It is the more accurate of the two methods in every configuration in which a difference is present, with gaps that are practically relevant under NLL and AIC at 100 and 200 samples, and it does so while removing a substantial number of nodes from the final equations, up to roughly 24 in the lowest sample regime\. The size reduction is the more systematic of the two effects, as it persists across all three criteria and at every sample size, whereas the accuracy advantage concentrates where the parameter estimates are least constrained by the data\. As in the previous analyses, the differences narrow at 500 samples and under BIC, where all methods essentially agree in accuracy and the remaining advantage of DSCA is one of parsimony\.

#### V\-B3Computational Efficiency

Training time is an important secondary metric for SR\. Since DSCA performsKKindependent inner optimizations per candidate, it is expected to be more computationally intensive\. Figure[2](https://arxiv.org/html/2610.10931#S5.F2)reports the training time of all evaluated algorithms111All experiments were run on a SLURM\-managed cluster whose nodes have two Intel Xeon Gold 5320 processors \(2\.2 GHz, 52 cores in total\) and 512 GB of shared memory\. To ensure a fair comparison, every job was allocated four CPU cores and 16 GiB of RAM\.\.

Despite the additional local optimizations, DSCA remains competitive with, or faster than, the centralized baselines in the lower\-sample regimes\. With 100 samples, it runs up to 20% faster than LM and BFGS under NLL and AIC\. With 200 samples, the differences are mixed, ranging from about 9% faster than BFGS under AIC to 12–20% slower depending on the criterion\. Only with 500 samples does the overhead of the inner loops fully materialize, with DSCA running 15–26% slower than the baselines across all criteria\.

This non\-linear scaling behavior occurs because numerical optimization is only one component of the SR pipeline, whose runtime is dominated by other operations such as symbolic simplification and offspring generation\. As these operations are faster on smaller expressions, the simpler populations maintained by DSCA accelerate the overall evolutionary process and, in low\-sample regimes, fully offset the cost of the inner optimization loops\.

Fig\. 2:Algorithmic training time \(in seconds\) across different selection criteria and sample size regimes\.

### V\-CMechanism Analysis: Selection\-Relevant Excess and Generalization

The preceding experiments show that DSCA tends to favor smaller expressions while retaining competitive predictive accuracy, particularly in low\-sample regimes\. We now examine whether this behavior is consistent with the mechanism characterized in Section[III\-C](https://arxiv.org/html/2610.10931#S3.SS3), namely the selection\-relevant excess lossΔN′\\Delta\_\{N^\{\\prime\}\}defined in \([4](https://arxiv.org/html/2610.10931#S3.E4)\)\. The objective is not to establish a causal identification of the mechanism, but to test two empirical implications of the theory: whether the structures favored by DSCA exhibit smaller out\-of\-sample deterioration under a common centralized calibration, and whetherΔN′\\Delta\_\{N^\{\\prime\}\}is associated with that deterioration across structures\.

We pool the final populations produced by DSCA, LM, and BFGS across all real\-world datasets, all three selection criteria, and all ten seeds in theN=200N=200regime, and remove duplicate structures\. Each surviving structure is then calibrated twice on the same optimization data: once centrally with BFGS, yieldingθ^N′\\hat\{\\theta\}\_\{N^\{\\prime\}\}and the centralized training lossLN′​\(θ^N′\)L\_\{N^\{\\prime\}\}\(\\hat\{\\theta\}\_\{N^\{\\prime\}\}\), and once with DSCA, yieldingθ¯\\bar\{\\theta\}andLN′​\(θ¯\)L\_\{N^\{\\prime\}\}\(\\bar\{\\theta\}\)\. Their difference isΔN′\\Delta\_\{N^\{\\prime\}\}\. Since each structure is scored by the same two procedures regardless of its origin, this design separates the properties of the structure from the properties of the optimizer that produced it, and lets us ask whetherΔN′\\Delta\_\{N^\{\\prime\}\}is informative about structures that DSCA itself never selected\.

Because the pooled datasets differ in target scale, all MSE quantities are divided by the median training MSE of all models on the same dataset, so thatΔN′\\Delta\_\{N^\{\\prime\}\}and the generalization gap are expressed relative to the typical error level of each dataset\.

Fig\. 3:Generalization gap of the pooled final populations, by the optimizer that produced each structure \(200\-sample regime, centralized re\-fit\)\.#### V\-C1Structures found under centralized optimization generalize worse

Figure[3](https://arxiv.org/html/2610.10931#S5.F3)reports the generalization gap of the pooled structures after all of them have been recalibrated centrally\. This common refit is important: differences in the gap then reflect the structures retained by the different evolutionary procedures rather than the particular coefficient values produced during their original runs\. Across the three outer selection criteria, structures originating from centralized LM and BFGS runs exhibit substantially larger median generalization gaps and heavier upper tails than structures originating from DSCA\. Combined with the predictive\-accuracy results of Section[V](https://arxiv.org/html/2610.10931#S5)\-B, this suggests that DSCA tends to retain structures whose fitted performance is less specific to the training sample, rather than merely producing structures with uniformly poorer in\-sample fit\.

#### V\-C2The excess loss is a signal of overfitting

Figure[4](https://arxiv.org/html/2610.10931#S5.F4)relatesΔN′\\Delta\_\{N^\{\\prime\}\}to the generalization gap over the same pooled set\. The association is strong and stable across criteria, with Pearson correlations between0\.720\.72and0\.830\.83, and it holds within each generating algorithm, not only across them: structures produced by LM with a largeΔN′\\Delta\_\{N^\{\\prime\}\}generalize as poorly as structures produced by any other optimizer with a comparableΔN′\\Delta\_\{N^\{\\prime\}\}\. This is the property the mechanism requires\.ΔN′\\Delta\_\{N^\{\\prime\}\}is computed without access to test data, from the training sample alone, and is available at no additional cost during the search, since DSCA usesLN′​\(θ¯\)L\_\{N^\{\\prime\}\}\(\\bar\{\\theta\}\)as its fitness in any case\.

![Refer to caption](https://arxiv.org/html/2610.10931v1/gen_gap_corr_norm.png)Fig\. 4:Excess lossΔN′\\Delta\_\{N^\{\\prime\}\}against the generalization gap for the pooled final populations, colored by the optimizer that produced the structure\.

### V\-DMechanism Controls: Optimization Effort and Partition Heterogeneity

The mechanism analyzed in Section[III\-C](https://arxiv.org/html/2610.10931#S3.SS3)attributes the behavior of DSCA to the averaging of local optima obtained on partitions drawn from environments that differ in their covariate distribution\. Two competing explanations must be excluded\. First, the advantage might arise simply because DSCA performs multiple optimization attempts per candidate rather than a single centralized fit\. Second, partitioning itself might be sufficient, irrespective of whether the subsets exhibit the controlled covariate heterogeneity targeted by DSCA\.

We therefore compare four variants under the 200\-sample regime\.*Multistart*performsKKindependent runs of the base optimizer on the full training set from random starting points and retains the best solution, so it matches the number of optimization attempts performed by DSCA without partitioning the data\.RbatchR\_\{\\text\{batch\}\}partitions the data uniformly at random\. The remaining two variants use the proposed Dirichlet\-Sinkhorn scheme withα=100\\alpha=100, which yields nearly homogeneous partitions, and withα=1\\alpha=1, which is the configuration used throughout the paper\. Each partition\-based configuration is evaluated withK=3K=3andK=5K=5\.

To prevent datasets of differing intrinsic difficulty from dominating the averages, we report a regret gap rather than raw scores\. For a given dataset and selection criterion, the gap measures how far a configuration’sR2R^\{2\}falls below the bestR2R^\{2\}attained by any configuration, and how far its complexity rises above the smallest complexity attained by any configuration\. A gap of zero denotes the best observed behavior and lower is better throughout\.

Figures[5](https://arxiv.org/html/2610.10931#S5.F5)and[6](https://arxiv.org/html/2610.10931#S5.F6)report the two gaps averaged over the real\-world datasets \(200 samples\)\. On predictive accuracy the multi\-start variant is the strongest or tied\-strongest configuration under every criterion, consistent with the advantage of selecting the best solution among several independent optimizations of the same pooled objective; unlike DSCA, Multistart does not intentionally modify the parameter vector at which the candidate is evaluated\. The partition\-based variants give up a small amount ofR2R^\{2\}relative to it, most visibly atK=3K=3under the unregularized criteria, and raisingKKto55removes most of this deficit\.

The ordering reverses on model complexity, and this is where the specific contribution of the proposed scheme becomes visible\. Multi\-start generally returns by far the largest expressions, and increasing the number of restarts does not reduce their size\. This is the expected consequence of a procedure that improves the fit of every candidate uniformly, including the overparameterized ones, and it confirms that the parsimony observed in Section[V\-B](https://arxiv.org/html/2610.10931#S5.SS2)cannot be attributed to the additional optimization effort alone\. Both partitioning strategies produce markedly smaller models, and among them the heterogeneous stratification is consistently the most parsimonious atK=5K=5, ahead of both random partitioning and the homogeneous split\. The three strategies are much closer atK=3K=3under NLL, which suggests that few partitions provide too little contrast between local optima for the averaging penalty to discriminate\.

The effect ofKKis asymmetric across the two axes\. Moving from three to five partitions leaves the accuracy gap essentially unchanged, or slightly improves it, while substantially reducing the complexity gap for the partition\-based variants; the same pattern is not observed for Multistart\. Taken together, the two figures support the reading advanced in Section[III\-C](https://arxiv.org/html/2610.10931#S3.SS3): a larger number of narrower environments increases the dispersion of the local optima and hence the loss incurred by averaging, and this penalty falls selectively on candidates whose parameters are unstable across environments\. The configuration used in the main experiments,α=1\\alpha=1withK=5K=5, gives up a small amount ofR2R^\{2\}relative to a pure restart strategy in exchange for a reduction of roughly 2 to 15 nodes\.

![Refer to caption](https://arxiv.org/html/2610.10931v1/ablation_regret_gap_multi.png)Fig\. 5:Ablation study of the predictive accuracy \(R2R^\{2\}\) regret gap under different partitioning strategies\.![Refer to caption](https://arxiv.org/html/2610.10931v1/ablation_complexity_gap_multi.png)Fig\. 6:Ablation study of the model complexity regret gap under different partitioning strategies\.

## VIConclusions

This paper introduced Dirichlet\-Sinkhorn Constant Averaging, a coefficient\-calibration strategy for genetic programming\-based SR that repurposes distributed parameter estimation as a source of selection pressure during evolutionary search\. DSCA calibrates each candidate independently on subsets whose covariate distributions differ by a controlled amount, and evaluates it using the average of the resulting parameter vectors\. Rather than improving the standard local search step, DSCA deliberately changes the parameterization at which candidate structures are compared during evolution\. The analysis identifies the excess lossΔN′\\Delta\_\{N^\{\\prime\}\}relative to centralized calibration as the quantity through which DSCA changes selection\. Its population counterpart vanishes for correctly specified and identifiable expressions, whereas misspecified candidates can incur a persistent penalty when their environment\-specific optima do not aggregate to the pooled optimum\.

The experiments are consistent with this interpretation\. DSCA achieves higher functional\-recovery rates than centralized BFGS and Levenberg\-Marquardt on the synthetic benchmarks and a more favorable accuracy\-complexity trade\-off on the real\-world datasets\. Its advantage is clearest in low\-sample regimes, where sample\-specific calibration and overfitting are most likely, and narrows as the sample size grows\. After a common centralized re\-fit, structures retained by DSCA exhibit smaller generalization gaps, andΔN′\\Delta\_\{N^\{\\prime\}\}is positively associated with this gap across structures and algorithms of origin\. Finally, multiple centralized restarts do not reproduce DSCA’s reduction in expression size, while score\-based heterogeneous partitioning yields the most parsimonious models atK=5K=5\. Together, these results support the role of controlled covariate heterogeneity rather than additional optimization effort alone\.

Two directions appear particularly important for future work\. First, symbolic simplification removes exact algebraic redundancies but does not prevent weakly identified or near\-singular parameterizations\. Structural repair operators that eliminate such degeneracies before calibration could make parameter averaging more reliable\. Second, the interaction between heterogeneous calibration and outer selection criteria such as AIC and BIC is explored only empirically here\. Characterizing how they jointly shape evolutionary pressure remains an important direction for further study\.

## References

- \[1\]G\. Kronberger, B\. Burlacu, M\. Kommenda, S\. M\. Winkler, and M\. Affenzeller,*Symbolic regression*\. Chapman and Hall/CRC, 2024\.
- \[2\]N\. Makke and S\. Chawla, “Interpretable scientific discovery with symbolic regression: a review,”*Artificial Intelligence Review*, vol\. 57, no\. 1, p\. 2, 2024\.
- \[3\]F\. O\. de Franca*et al\.*, “SRBench\+\+: Principled benchmarking of symbolic regression with domain\-expert interpretation,”*IEEE Transactions on Evolutionary Computation*, vol\. 29, no\. 4, pp\. 1127–1137, 2025\.
- \[4\]J\. R\. Koza, “Genetic programming as a means for programming computers by natural selection,”*Statistics and Computing*, vol\. 4, no\. 2, pp\. 87–112, Jun 1994\. \[Online\]\. Available:[https://doi\.org/10\.1007/BF00175355](https://doi.org/10.1007/BF00175355)
- \[5\]D\. Sobania, M\. Briesch, D\. Wittenberg, and F\. Rothlauf, “Analyzing optimized constants in genetic programming on a real\-world regression problem,” in*Proceedings of the Genetic and Evolutionary Computation Conference Companion*, ser\. GECCO ’22\. New York, NY, USA: Association for Computing Machinery, 2022, p\. 606–607\. \[Online\]\. Available:[https://doi\.org/10\.1145/3520304\.3528896](https://doi.org/10.1145/3520304.3528896)
- \[6\]A\. Topchy and W\. F\. Punch, “Faster genetic programming based on local gradient search of numeric leaf values,” in*Proceedings of the 3rd Annual Conference on Genetic and Evolutionary Computation*, ser\. GECCO’01\. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc\., 2001, p\. 155–162\.
- \[7\]M\. Kommenda, B\. Burlacu, G\. Kronberger, and M\. Affenzeller, “Parameter identification for symbolic regression using nonlinear least squares,”*Genetic Programming and Evolvable Machines*, vol\. 21, no\. 3, pp\. 471–501, Sep 2020\. \[Online\]\. Available:[https://doi\.org/10\.1007/s10710\-019\-09371\-3](https://doi.org/10.1007/s10710-019-09371-3)
- \[8\]M\. Cranmer, “Interpretable machine learning for science with PySR and SymbolicRegression\.jl,” arXiv:2305\.01582 \[astro\-ph\.IM\], 2023\. \[Online\]\. Available:[https://arxiv\.org/abs/2305\.01582](https://arxiv.org/abs/2305.01582)
- \[9\]G\. Kronberger, “Local optimization often is ill\-conditioned in genetic programming for symbolic regression,” in*2022 24th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing \(SYNASC\)*, 2022, pp\. 304–310\.
- \[10\]J\. Harrison, M\. Virgolin, T\. Alderliesten, and P\. Bosman, “Mini\-batching, gradient\-clipping, first\- versus second\-order: What works in gradient\-based coefficient optimisation for symbolic regression?” in*Proceedings of the Genetic and Evolutionary Computation Conference*, ser\. GECCO ’23\. New York, NY, USA: Association for Computing Machinery, 2023, p\. 1127–1136\. \[Online\]\. Available:[https://doi\.org/10\.1145/3583131\.3590368](https://doi.org/10.1145/3583131.3590368)
- \[11\]A\. Soltani, G\. Kronberger, F\. Olivetti de Franca, M\. Billa, and A\. Lucantonio, “A comparative study of model selection criteria for symbolic regression,” in*Proceedings of the Genetic and Evolutionary Computation Conference*, 2026, pp\. 826–833\.
- \[12\]S\. Watanabe, “Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory,”*Journal of Machine Learning Research*, vol\. 11, pp\. 3571–3594, 2010\.
- \[13\]A\. Vehtari, A\. Gelman, and J\. Gabry, “Practical Bayesian model evaluation using leave\-one\-out cross\-validation and WAIC,”*Statistics and Computing*, vol\. 27, no\. 5, pp\. 1413–1432, 2017\.
- \[14\]B\. McMahan, E\. Moore, D\. Ramage, S\. Hampson, and B\. A\. y\. Arcas, “Communication\-Efficient Learning of Deep Networks from Decentralized Data,” in*Proceedings of the 20th International Conference on Artificial Intelligence and Statistics*, ser\. Proceedings of Machine Learning Research, A\. Singh and J\. Zhu, Eds\., vol\. 54\. PMLR, 20–22 Apr 2017, pp\. 1273–1282\. \[Online\]\. Available:[https://proceedings\.mlr\.press/v54/mcmahan17a\.html](https://proceedings.mlr.press/v54/mcmahan17a.html)
- \[15\]M\. Kommenda, M\. Affenzeller, G\. Kronberger, and S\. M\. Winkler, “Nonlinear least squares optimization of constants in symbolic regression,” in*Computer Aided Systems Theory \- EUROCAST 2013*, R\. Moreno\-Díaz, F\. Pichler, and A\. Quesada\-Arencibia, Eds\. Berlin, Heidelberg: Springer Berlin Heidelberg, 2013, pp\. 420–427\.
- \[16\]M\. Kommenda, G\. Kronberger, S\. Winkler, M\. Affenzeller, and S\. Wagner, “Effects of constant optimization by nonlinear least squares minimization in symbolic regression,” in*Proceedings of the 15th Annual Conference Companion on Genetic and Evolutionary Computation*, ser\. GECCO ’13 Companion\. New York, NY, USA: Association for Computing Machinery, 2013, p\. 1121–1128\. \[Online\]\. Available:[https://doi\.org/10\.1145/2464576\.2482691](https://doi.org/10.1145/2464576.2482691)
- \[17\]P\. Moscato*et al\.*, “On evolution, search, optimization, genetic algorithms and martial arts: Towards memetic algorithms,”*Caltech concurrent computation program, C3P Report*, vol\. 826, no\. 1989, p\. 37, 1989\.
- \[18\]L\. Trujillo*et al\.*,*Local Search is Underused in Genetic Programming*\. Cham: Springer International Publishing, 2018, pp\. 119–137\. \[Online\]\. Available:[https://doi\.org/10\.1007/978\-3\-319\-97088\-2\_8](https://doi.org/10.1007/978-3-319-97088-2_8)
- \[19\]D\. W\. Marquardt, “An algorithm for least\-squares estimation of nonlinear parameters,”*Journal of the Society for Industrial and Applied Mathematics*, vol\. 11, no\. 2, pp\. 431–441, 1963\. \[Online\]\. Available:[http://www\.jstor\.org/stable/2098941](http://www.jstor.org/stable/2098941)
- \[20\]B\. Burlacu, G\. Kronberger, and M\. Kommenda, “Operon c\+\+: an efficient genetic programming framework for symbolic regression,” in*Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion*, ser\. GECCO ’20\. New York, NY, USA: Association for Computing Machinery, 2020, p\. 1562–1570\. \[Online\]\. Available:[https://doi\.org/10\.1145/3377929\.3398099](https://doi.org/10.1145/3377929.3398099)
- \[21\]C\. G\. Broyden, “The convergence of a class of double\-rank minimization algorithms 1\. general considerations,”*IMA Journal of Applied Mathematics*, vol\. 6, no\. 1, pp\. 76–90, 1970\.
- \[22\]G\. F\. Bomarito, P\. E\. Leser, N\. C\. M\. Strauss, K\. M\. Garbrecht, and J\. D\. Hochhalter, “Bayesian model selection for reducing bloat and overfitting in genetic programming for symbolic regression,” in*Proceedings of the Genetic and Evolutionary Computation Conference Companion*, ser\. GECCO ’22\. New York, NY, USA: Association for Computing Machinery, 2022, p\. 526–529\. \[Online\]\. Available:[https://doi\.org/10\.1145/3520304\.3528899](https://doi.org/10.1145/3520304.3528899)
- \[23\]G\. Bomarito, P\. Leser, N\. Strauss, K\. Garbrecht, and J\. Hochhalter, “Automated learning of interpretable models with quantified uncertainty,”*Computer Methods in Applied Mechanics and Engineering*, vol\. 403, p\. 115732, 2023\. \[Online\]\. Available:[https://www\.sciencedirect\.com/science/article/pii/S0045782522006879](https://www.sciencedirect.com/science/article/pii/S0045782522006879)
- \[24\]M\. Billa, V\. Guidetti, L\. La Rocca, and F\. Mandreoli, “Mining trustworthy symbolic regression models in federated settings,” in*2025 IEEE International Conference on Data Mining \(ICDM\)*, 2025, pp\. 1055–1064\.
- \[25\]L\. G\. dos Reis, V\. L\. P\. S\. Caminha, and T\. J\. P\. Penna, “Benchmarking symbolic regression constant optimization schemes,” arXiv:2412\.02126, 2024\. \[Online\]\. Available:[https://arxiv\.org/abs/2412\.02126](https://arxiv.org/abs/2412.02126)
- \[26\]G\. Schwarz, “Estimating the dimension of a model,”*The Annals of Statistics*, vol\. 6, no\. 2, pp\. 461–464, 1978\.
- \[27\]R\. E\. Kass and A\. E\. Raftery, “Bayes factors,”*Journal of the American Statistical Association*, vol\. 90, no\. 430, pp\. 773–795, 1995\.
- \[28\]G\. Kronberger, F\. O\. de França, D\. J\. Bartlett, H\. Desmond, and P\. G\. Ferreira, “Guiding multi\-objective genetic programming with description length improves symbolic regression solutions,” arXiv:2605\.22374, 2026\. \[Online\]\. Available:[https://arxiv\.org/abs/2605\.22374](https://arxiv.org/abs/2605.22374)
- \[29\]S\. Watanabe,*Algebraic Geometry and Statistical Learning Theory*\. Cambridge University Press, 2009\.
- \[30\]——, “A widely applicable Bayesian information criterion,”*Journal of Machine Learning Research*, vol\. 14, pp\. 867–897, 2013\.
- \[31\]J\. M\. Bernardo and A\. F\. M\. Smith,*Bayesian Theory*\. Chichester: Wiley, 1994\.
- \[32\]Y\. Zhang, J\. Duchi, and M\. Wainwright, “Divide and conquer kernel ridge regression,” in*Proceedings of the 26th Annual Conference on Learning Theory*, ser\. Proceedings of Machine Learning Research, S\. Shalev\-Shwartz and I\. Steinwart, Eds\., vol\. 30\. Princeton, NJ, USA: PMLR, 12–14 Jun 2013, pp\. 592–617\. \[Online\]\. Available:[https://proceedings\.mlr\.press/v30/Zhang13\.html](https://proceedings.mlr.press/v30/Zhang13.html)
- \[33\]——, “Communication\-efficient algorithms for statistical optimization,”*J\. Mach\. Learn\. Res\.*, vol\. 14, no\. 1, p\. 3321–3363, Jan\. 2013\.
- \[34\]J\. D\. Rosenblatt and B\. Nadler, “On the optimality of averaging in distributed statistical learning,”*Information and Inference: A Journal of the IMA*, vol\. 5, no\. 4, pp\. 379–404, 12 2016\. \[Online\]\. Available:[https://doi\.org/10\.1093/imaiai/iaw013](https://doi.org/10.1093/imaiai/iaw013)
- \[35\]T\.\-M\. H\. Hsu, H\. Qi, and M\. Brown, “Measuring the effects of non\-identical data distribution for federated visual classification,” arXiv:1909\.06335, 2019\. \[Online\]\. Available:[https://arxiv\.org/abs/1909\.06335](https://arxiv.org/abs/1909.06335)
- \[36\]J\. Dong, J\. Zhong, W\.\-N\. Chen, and J\. Zhang, “An efficient federated genetic programming framework for symbolic regression,”*IEEE Transactions on Emerging Topics in Computational Intelligence*, vol\. 7, no\. 3, pp\. 858–871, 2023\.
- \[37\]R\. Sinkhorn and P\. Knopp, “Concerning nonnegative matrices and doubly stochastic matrices,”*Pacific Journal of Mathematics*, vol\. 21, no\. 2, pp\. 343–348, 1967\.
- \[38\]F\. W\. Ramlan, G\. Kronberger, C\. O’Riordan, and J\. McDermott, “Comparative analysis of model selection criteria for symbolic regression using genetic programming,” in*Computational Intelligence*, F\. Marcelloni, K\. Madani, N\. van Stein, and J\. Filipe, Eds\. Cham: Springer Nature Switzerland, 2026, pp\. 91–108\.
- \[39\]Y\. Jin, W\. Fu, J\. Kang, J\. Guo, and J\. Guo, “Bayesian symbolic regression,” arXiv:1910\.08892, 2019\. \[Online\]\. Available:[https://arxiv\.org/abs/1910\.08892](https://arxiv.org/abs/1910.08892)
- \[40\]Q\. Chen, B\. Xue, and M\. Zhang, “Generalisation and domain adaptation in gp with gradient descent for symbolic regression,” in*2015 IEEE congress on evolutionary computation \(CEC\)*\. IEEE, 2015, pp\. 1137–1144\.
- \[41\]E\. J\. Vladislavleva, G\. F\. Smits, and D\. Den Hertog, “Order of nonlinearity as a complexity measure for models generated by symbolic regression via pareto genetic programming,”*IEEE Transactions on Evolutionary Computation*, vol\. 13, no\. 2, pp\. 333–349, 2008\.
- \[42\]R\. S\. Olson, W\. La Cava, P\. Orzechowski, R\. J\. Urbanowicz, and J\. H\. Moore, “Pmlb: a large benchmark suite for machine learning evaluation and comparison,”*BioData mining*, vol\. 10, no\. 1, p\. 36, 2017\.
- \[43\]W\. La Cava*et al\.*, “Contemporary symbolic regression methods and their relative performance,” in*Proc\. NeurIPS Track on Datasets and Benchmarks*, 2021\.
- \[44\]A\. Benavoli, G\. Corani, F\. Mangili, M\. Zaffalon, and F\. Ruggeri, “A bayesian wilcoxon signed\-rank test based on the dirichlet process,” in*International conference on machine learning*\. PMLR, 2014, pp\. 1026–1034\.
- \[45\]A\. Benavoli, G\. Corani, J\. Demšar, and M\. Zaffalon, “Time for a change: a tutorial for comparing multiple classifiers through bayesian analysis,”*Journal of Machine Learning Research*, vol\. 18, no\. 77, pp\. 1–36, 2017\.

## Appendix ATheoretical Analysis of Averaging Under Heterogeneity

Throughout this appendix,NNdenotes the size of the optimization set𝒟opt\\mathcal\{D\}\_\{\\mathrm\{opt\}\}, denotedN′N^\{\\prime\}in the main text, andn=N/Kn=N/K\. The distinctive effect of DSCA arises from the interaction between heterogeneous partitions and parameter averaging\. For a fixed candidate expression, DSCA optimizes its constants independently onKKpartitions, averages the resulting parameter estimates, and evaluates the averaged parameter vector on the complete training set\. This procedure can increase the loss assigned to candidates whose preferred parameters vary across partitions\. Two quantities are useful for understanding this effect\. The first measures disagreement among locally optimal solutions and explains the mechanism induced by averaging\. The second measures the actual increase in loss relative to centralized constant optimization and is therefore the quantity directly relevant to evolutionary selection\.

*Empirical decomposition\.*For a candidate expressionf⁡\(⋅,θ\)f\(\\cdot,\\theta\), let

Lk​\(θ\):=1n​∑i∈𝒟kℓ⁡\(zi,θ\),LN​\(θ\):=1N​∑i=1Nℓ⁡\(zi,θ\),L\_\{k\}\(\\theta\):=\\frac\{1\}\{n\}\\sum\_\{i\\in\\mathcal\{D\}\_\{k\}\}\\ell\(z\_\{i\},\\theta\),\\qquad L\_\{N\}\(\\theta\):=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\ell\(z\_\{i\},\\theta\),\(12\)wherezi=\(xi,yi\)z\_\{i\}=\(x\_\{i\},y\_\{i\}\)andn=N/Kn=N/K\. Because theKKpartitions are disjoint and equally sized,

LN​\(θ\)=1K​∑k=1KLk​\(θ\)L\_\{N\}\(\\theta\)=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}L\_\{k\}\(\\theta\)\(13\)holds exactly for everyθ\\theta\. The equal\-size constraint is enforced by the Sinkhorn scaling; with unequal partitions the same identity holds with weightsnk/Nn\_\{k\}/N\. Let

θ^k:=argminθLk\(θ\),θ¯:=1K∑k=1Kθ^k,θ^N:=arg⁡minθ​LN​\(θ\)\.\\begin\{array\}\[\]\{cc\}\\hat\{\\theta\}\_\{k\}:=\\arg\\min\_\{\\theta\}L\_\{k\}\(\\theta\),\\qquad\\bar\{\\theta\}:=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\hat\{\\theta\}\_\{k\},\\\\\[15\.0pt\] \\hat\{\\theta\}\_\{N\}:=\\arg\\min\_\{\\theta\}L\_\{N\}\(\\theta\)\.\\end\{array\}\(14\)Evaluating \([13](https://arxiv.org/html/2610.10931#A1.E13)\) atθ¯\\bar\{\\theta\}gives

LN​\(θ¯\)=1K​∑k=1KLk​\(θ^k\)\+DN,DN:=1K​∑k=1K\[Lk​\(θ¯\)−Lk​\(θ^k\)\]\.\\begin\{array\}\[\]\{cc\}L\_\{N\}\(\\bar\{\\theta\}\)=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}L\_\{k\}\(\\hat\{\\theta\}\_\{k\}\)\+D\_\{N\},\\\\\[15\.0pt\] D\_\{N\}:=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\bigl\[L\_\{k\}\(\\bar\{\\theta\}\)\-L\_\{k\}\(\\hat\{\\theta\}\_\{k\}\)\\bigr\]\.\\end\{array\}\(15\)If eachθ^k\\hat\{\\theta\}\_\{k\}is a global minimizer of its local objective, thenDN≥0D\_\{N\}\\geq 0\. The quantity that measures the loss added by DSCA relative to centralized parameter optimization is

ΔN:=LN​\(θ¯\)−LN​\(θ^N\)≥0\.\\Delta\_\{N\}:=L\_\{N\}\(\\bar\{\\theta\}\)\-L\_\{N\}\(\\hat\{\\theta\}\_\{N\}\)\\geq 0\.\(16\)
ΔN\\Delta\_\{N\}is an analytical quantity rather than an additional component of the algorithm: DSCA assigns the fitnessLN​\(θ¯\)L\_\{N\}\(\\bar\{\\theta\}\), whereasΔN\\Delta\_\{N\}measures its excess over the loss obtained by centrally calibrating the same structure\. Define also

FN:=LN​\(θ^N\)−1K​∑k=1KLk​\(θ^k\)\.F\_\{N\}:=L\_\{N\}\(\\hat\{\\theta\}\_\{N\}\)\-\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}L\_\{k\}\(\\hat\{\\theta\}\_\{k\}\)\.\(17\)Then

DN=ΔN\+FN\\boxed\{\\;D\_\{N\}=\\Delta\_\{N\}\+F\_\{N\}\\;\}\(18\)Under exact minimization,DND\_\{N\},ΔN\\Delta\_\{N\}, andFNF\_\{N\}are non\-negative\. The decomposition distinguishes disagreement among local empirical optima from the component that actually changes the loss relative to centralized calibration\. In particular,DN\>0D\_\{N\}\>0does not implyΔN\>0\\Delta\_\{N\}\>0\.

*Covariate\-induced environments\.*LetPPdenote the data\-generating law and letu=g⁡\(x\)u=g\(x\)be the stratification score of Section[III](https://arxiv.org/html/2610.10931#S3), withgga function of the covariates alone, so thatuuis measurable with respect toσ⁡\(x\)\\sigma\(x\)\. At the population level, letπk​\(u\)\\pi\_\{k\}\(u\)denote the probability of assigning an observation with scoreuuto partitionkk\. The balancing constraints imply

∑k=1Kπk​\(u\)=1,𝔼P​\[πk​\(g⁡\(X\)\)\]=1K\.\\sum\_\{k=1\}^\{K\}\\pi\_\{k\}\(u\)=1,\\qquad\\mathbb\{E\}\_\{P\}\[\\pi\_\{k\}\(g\(X\)\)\]=\\frac\{1\}\{K\}\.\(19\)Environmentkkis the lawPkP\_\{k\}from which𝒟k\\mathcal\{D\}\_\{k\}is drawn, with density

pk​\(x,y\)=K​πk​\(g⁡\(x\)\)​p​\(x,y\)\.p\_\{k\}\(x,y\)=K\\,\\pi\_\{k\}\(g\(x\)\)\\,p\(x,y\)\.\(20\)Consequently,

1K​∑k=1Kpk​\(x,y\)=p⁡\(x,y\),pk​\(y∣x\)=p⁡\(y∣x\)\.\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}p\_\{k\}\(x,y\)=p\(x,y\),\\qquad p\_\{k\}\(y\\mid x\)=p\(y\\mid x\)\.\(21\)Thus the induced environments differ in their covariate marginals while preserving the conditional data\-generating mechanism\.

*Population quantities\.*Define

Λk​\(θ\):=𝔼Pk​\[ℓ⁡\(Z,θ\)\],Λ⁡\(θ\):=𝔼P​\[ℓ⁡\(Z,θ\)\]=1K​∑k=1KΛk​\(θ\),\\begin\{array\}\[\]\{cc\}\\Lambda\_\{k\}\(\\theta\):=\\mathbb\{E\}\_\{P\_\{k\}\}\[\\ell\(Z,\\theta\)\],\\\\\[15\.0pt\] \\Lambda\(\\theta\):=\\mathbb\{E\}\_\{P\}\[\\ell\(Z,\\theta\)\]=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\Lambda\_\{k\}\(\\theta\),\\end\{array\}\(22\)the last equality following from \([21](https://arxiv.org/html/2610.10931#A1.E21)\)\. Let

θk∗:=arg⁡minθ​Λk​\(θ\),θ¯∗:=1K​∑k=1Kθk∗,θ∗:=arg⁡minθ⁡Λ⁡\(θ\)\.\\begin\{array\}\[\]\{ccc\}\\theta\_\{k\}^\{\*\}:=\\arg\\min\_\{\\theta\}\\Lambda\_\{k\}\(\\theta\),\\\\\[10\.0pt\] \\bar\{\\theta\}^\{\*\}:=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\theta\_\{k\}^\{\*\},\\\\\[10\.0pt\] \\theta^\{\*\}:=\\arg\\min\_\{\\theta\}\\Lambda\(\\theta\)\.\\end\{array\}\(23\)Define the selection\-relevant population excess

Δ:=Λ⁡\(θ¯∗\)−Λ⁡\(θ∗\)≥0\.\\Delta:=\\Lambda\(\\bar\{\\theta\}^\{\*\}\)\-\\Lambda\(\\theta^\{\*\}\)\\geq 0\.\(24\)This is the population counterpart ofΔN\\Delta\_\{N\}and is the quantity whose behavior is characterized below\.

###### Proposition 1\(Correct specification under covariate reweighting\)\.

Suppose that the partitioning score is a function of the covariates only,u=g⁡\(x\)u=g\(x\), and that there exists an identifiableθ∗\\theta^\{\*\}satisfying

f⁡\(X,θ∗\)=𝔼⁡\[Y∣X\]f\(X,\\theta^\{\*\}\)=\\mathbb\{E\}\[Y\\mid X\]PP\-almost surely\. Ifπk​\(g​\(x\)\)\>0\\pi\_\{k\}\(g\(x\)\)\>0forPP\-almost everyxxand everykk, thenθk∗=θ∗\\theta\_\{k\}^\{\*\}=\\theta^\{\*\}for everykk, and consequentlyΔ=0\\Delta=0\.

Proof\.Because the partitioning score depends only onxx, \([20](https://arxiv.org/html/2610.10931#A1.E20)\) factorizes as

pk​\(x,y\)=K​πk​\(g⁡\(x\)\)​p​\(x\)​p​\(y∣x\)=pk​\(x\)​p​\(y∣x\),pk​\(x\):=K​πk​\(g⁡\(x\)\)​p​\(x\),\\begin\{array\}\[\]\{ll\}p\_\{k\}\(x,y\)=K\\,\\pi\_\{k\}\(g\(x\)\)\\,p\(x\)\\,p\(y\\mid x\)=p\_\{k\}\(x\)\\,p\(y\\mid x\),\\\\\[10\.0pt\] p\_\{k\}\(x\):=K\\,\\pi\_\{k\}\(g\(x\)\)\\,p\(x\),\\end\{array\}\(25\)sopk​\(y∣x\)=p⁡\(y∣x\)p\_\{k\}\(y\\mid x\)=p\(y\\mid x\)for every environment\. The local population risk can therefore be written as

Λk\(θ\)=𝔼Pk\[H\(p\(⋅∣X\)\)\]\+𝔼Pk\[KL\(p\(⋅∣X\)∥pθ\(⋅∣X\)\)\]\.\\Lambda\_\{k\}\(\\theta\)=\\mathbb\{E\}\_\{P\_\{k\}\}\\\!\\left\[\\mathrm\{H\}\\bigl\(p\(\\cdot\\mid X\)\\bigr\)\\right\]\+\\mathbb\{E\}\_\{P\_\{k\}\}\\\!\\left\[\\mathrm\{KL\}\\bigl\(p\(\\cdot\\mid X\)\\,\\\|\\,p\_\{\\theta\}\(\\cdot\\mid X\)\\bigr\)\\right\]\.\(26\)The first term is independent ofθ\\theta, so allθ\\theta\-dependence of every local objective is carried by

gθ\(x\):=KL\(p\(⋅∣x\)∥pθ\(⋅∣x\)\)≥0\.g\_\{\\theta\}\(x\):=\\mathrm\{KL\}\\bigl\(p\(\\cdot\\mid x\)\\,\\\|\\,p\_\{\\theta\}\(\\cdot\\mid x\)\\bigr\)\\geq 0\.Under correct specification,gθ∗​\(x\)=0g\_\{\\theta^\{\*\}\}\(x\)=0forPP\-almost everyxx\. Sincepkp\_\{k\}is obtained fromppthrough the non\-negative bounded weightK​πk​\(g⁡\(x\)\)∈\[0,K\]K\\pi\_\{k\}\(g\(x\)\)\\in\[0,K\], everyPP\-null set is alsoPkP\_\{k\}\-null and

𝔼Pk​\[gθ∗​\(X\)\]=0≤𝔼Pk​\[gθ​\(X\)\]∀θ,∀k\.\\mathbb\{E\}\_\{P\_\{k\}\}\[g\_\{\\theta^\{\*\}\}\(X\)\]=0\\leq\\mathbb\{E\}\_\{P\_\{k\}\}\[g\_\{\\theta\}\(X\)\]\\qquad\\forall\\theta,\\ \\forall k\.\(27\)Henceθ∗\\theta^\{\*\}minimizes every local population risk\. Identifiability givesθk∗=θ∗\\theta\_\{k\}^\{\*\}=\\theta^\{\*\}for everykk, and thereforeθ¯∗=θ∗\\bar\{\\theta\}^\{\*\}=\\theta^\{\*\}andΔ=0\\Delta=0\. □\\square

*Relation with squared\-error fitting\.*Under squared\-error loss,

Λk​\(θ\)=𝔼Pk​\[Var​\(Y∣X\)\]⏟independent of​θ\+𝔼Pk​\[\(𝔼⁡\[Y∣X\]−f⁡\(X,θ\)\)2\]\.\\begin\{array\}\[\]\{lll\}\\Lambda\_\{k\}\(\\theta\)=&\\underbrace\{\\mathbb\{E\}\_\{P\_\{k\}\}\\bigl\[\\mathrm\{Var\}\(Y\\mid X\)\\bigr\]\}\_\{\\text\{independent of \}\\theta\}\\\\\[20\.0pt\] &\+\\mathbb\{E\}\_\{P\_\{k\}\}\\\!\\left\[\\bigl\(\\mathbb\{E\}\[Y\\mid X\]\-f\(X,\\theta\)\\bigr\)^\{2\}\\right\]\.\\end\{array\}\(28\)The first term is independent ofθ\\theta\. Proposition[1](https://arxiv.org/html/2610.10931#Thmproposition1)therefore extends directly to squared\-error calibration whenever the conditional mean is correctly specified; correct specification of the complete conditional noise distribution is not required\.

*Misspecification and selection\-relevant instability\.*For a misspecified candidate, the local pseudo\-true parameter is

θk∗=argminθ𝔼Pk\[KL\(p\(⋅∣X\)∥pθ\(⋅∣X\)\)\],\\theta\_\{k\}^\{\*\}=\\arg\\min\_\{\\theta\}\\mathbb\{E\}\_\{P\_\{k\}\}\\\!\\left\[\\mathrm\{KL\}\\bigl\(p\(\\cdot\\mid X\)\\,\\\|\\,p\_\{\\theta\}\(\\cdot\\mid X\)\\bigr\)\\right\],\(29\)a weighted KL projection\. If improving the approximation in one region of covariate space necessarily worsens it in another, the optimal resolution of that trade\-off depends on the weights assigned to those regions, and the local pseudo\-true parameters may therefore differ\. Misspecification alone is insufficient for a population\-level DSCA effect\. A misspecified structure may share the same pseudo\-true optimum across the induced environments, and even when theθk∗\\theta\_\{k\}^\{\*\}differ, their arithmetic mean may still coincide with the pooled pseudo\-true parameter\. The selection\-relevant condition is therefore

θ¯∗≠θ∗\.\\bar\{\\theta\}^\{\*\}\\neq\\theta^\{\*\}\.\(30\)
*Local quadratic interpretation\.*Suppose that, near the relevant optima,

Λk​\(θ\)=ck\+12​\(θ−θk∗\)⊤​Hk​\(θ−θk∗\),Hk≻0\.\\Lambda\_\{k\}\(\\theta\)=c\_\{k\}\+\\frac\{1\}\{2\}\(\\theta\-\\theta\_\{k\}^\{\*\}\)^\{\\top\}H\_\{k\}\(\\theta\-\\theta\_\{k\}^\{\*\}\),\\qquad H\_\{k\}\\succ 0\.\(31\)Let

H¯:=1K​∑k=1KHk,dk:=θk∗−θ¯∗\.\\bar\{H\}:=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}H\_\{k\},\\qquad d\_\{k\}:=\\theta\_\{k\}^\{\*\}\-\\bar\{\\theta\}^\{\*\}\.Since∑kdk=0\\sum\_\{k\}d\_\{k\}=0, the pooled first\-order condition gives

θ¯∗−θ∗=−H¯−11K∑k=1K\(Hk−H¯\)dk\.\\bar\{\\theta\}^\{\*\}\-\\theta^\{\*\}=\-\\bar\{H\}^\{\-1\}\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\(H\_\{k\}\-\\bar\{H\}\)d\_\{k\}\.\(32\)Consequently,

Δ=12​‖1K​∑k=1K\(Hk−H¯\)​dk‖H¯−12\.\\Delta=\\frac\{1\}\{2\}\\left\\\|\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\(H\_\{k\}\-\\bar\{H\}\)d\_\{k\}\\right\\\|\_\{\\bar\{H\}^\{\-1\}\}^\{2\}\.\(33\)Thus disagreement among the preferred local optima is insufficient for a selection\-relevant population penalty\. In the quadratic approximation,Δ\\Deltaarises from the interaction between their displacement and differences in local curvature\. If allHkH\_\{k\}coincide, thenθ¯∗=θ∗\\bar\{\\theta\}^\{\*\}=\\theta^\{\*\}andΔ=0\\Delta=0even when the local optima differ substantially\. For non\-quadratic risks, higher\-order differences in local objective shape can play an analogous role\.

Homogeneous\-limit behavior\.Asα→∞\\alpha\\to\\infty, the allocation probabilities converge toπk​\(u\)=1/K\\pi\_\{k\}\(u\)=1/K, soPk→PP\_\{k\}\\to Pfor everykkand consequentlyθk∗→θ∗\\theta\_\{k\}^\{\*\}\\to\\theta^\{\*\}andΔ→0\\Delta\\to 0\. To characterize the rate, write

πk​\(u\)=K−1\+εk​\(u\),∑kεk​\(u\)≡0,\\pi\_\{k\}\(u\)=K^\{\-1\}\+\\varepsilon\_\{k\}\(u\),\\qquad\\sum\_\{k\}\\varepsilon\_\{k\}\(u\)\\equiv 0,and letrαr\_\{\\alpha\}denote the typical size of the perturbations\. Then

Λk=Λ\+K​Gk,Gk​\(θ\):=𝔼P​\[εk​\(g⁡\(X\)\)​ℓ​\(Z,θ\)\],∑kGk≡0\.\\begin\{array\}\[\]\{lll\}\\Lambda\_\{k\}=\\Lambda\+K\\,G\_\{k\},\\\\\[10\.0pt\] G\_\{k\}\(\\theta\):=\\mathbb\{E\}\_\{P\}\[\\varepsilon\_\{k\}\(g\(X\)\)\\ell\(Z,\\theta\)\],\\\\\[10\.0pt\] \\sum\_\{k\}G\_\{k\}\\equiv 0\.\\end\{array\}Let

H:=∇2Λ​\(θ∗\)H:=\\nabla^\{2\}\\Lambda\(\\theta^\{\*\}\)and assumeH≻0H\\succ 0\. Under smoothness, perturbing the stationarity condition∇Λk​\(θk∗\)=0\\nabla\\Lambda\_\{k\}\(\\theta\_\{k\}^\{\*\}\)=0gives

θk∗−θ∗=−KH−1∇Gk\(θ∗\)\+O\(rα2\)=O\(rα\)\.\\theta\_\{k\}^\{\*\}\-\\theta^\{\*\}=\-KH^\{\-1\}\\nabla G\_\{k\}\(\\theta^\{\*\}\)\+O\(r\_\{\\alpha\}^\{2\}\)=O\(r\_\{\\alpha\}\)\.\(34\)Because

∑k=1K∇Gk​\(θ∗\)=∇\(∑k=1KGk\)​\(θ∗\)=0,\\sum\_\{k=1\}^\{K\}\\nabla G\_\{k\}\(\\theta^\{\*\}\)=\\nabla\\left\(\\sum\_\{k=1\}^\{K\}G\_\{k\}\\right\)\(\\theta^\{\*\}\)=0,the first\-order perturbations cancel when the local optima are arithmetically averaged:

θ¯∗−θ∗=O⁡\(rα2\),Δ=O⁡\(rα4\)\.\\bar\{\\theta\}^\{\*\}\-\\theta^\{\*\}=O\(r\_\{\\alpha\}^\{2\}\),\\qquad\\Delta=O\(r\_\{\\alpha\}^\{4\}\)\.\(35\)Before Sinkhorn balancing, a component of a symmetric Dirichlet vector has variance

Var⁡\[Ms,k\]=K−1K2​\(K​α\+1\)\.\\operatorname\{Var\}\[M\_\{s,k\}\]=\\frac\{K\-1\}\{K^\{2\}\(K\\alpha\+1\)\}\.\(36\)The final allocation probabilities after Sinkhorn balancing are not themselves Dirichlet distributed\. If the balancing map is smooth in a neighborhood of the uniform allocation, the balanced perturbations retain the same asymptotic scale, giving

rα2=O⁡\(\(K​α\+1\)−1\)as​α→∞\.r\_\{\\alpha\}^\{2\}=O\\\!\\left\(\(K\\alpha\+1\)^\{\-1\}\\right\)\\qquad\\text\{as \}\\alpha\\to\\infty\.\(37\)Together with \([35](https://arxiv.org/html/2610.10931#A1.E35)\), this explains the rapid attenuation of the selection\-relevant population effect near homogeneous partitioning\. The result is local to the regimeα→∞\\alpha\\to\\inftyand does not imply monotonicity ofΔ\\Deltainα\\alphaor identify a universal optimal heterogeneity level\.

*Finite\-sample connection\.*Assume fixedKKand standard joint M\-estimation regularity conditions for the balanced partition design\. Each local empirical minimizer is root\-NNconsistent for its environment\-specific population optimum, and therefore

θ¯−θ¯∗=Op\(N−1/2\)\.\\bar\{\\theta\}\-\\bar\{\\theta\}^\{\*\}=O\_\{p\}\(N^\{\-1/2\}\)\.\(38\)WhenΔ\>0\\Delta\>0, smoothness of the empirical and population risks then gives

ΔN=Δ\+Op\(N−1/2\)\.\\Delta\_\{N\}=\\Delta\+O\_\{p\}\(N^\{\-1/2\}\)\.\(39\)Thus the empirical excess fluctuates around a non\-vanishing population quantity\. Under correct specification, all environments instead share the same population optimumθ∗\\theta^\{\*\}\. Let

Sk:=∇Lk​\(θ∗\),Hk:=∇2Λk​\(θ∗\),H¯:=1K​∑k=1KHk\.S\_\{k\}:=\\nabla L\_\{k\}\(\\theta^\{\*\}\),\\quad H\_\{k\}:=\\nabla^\{2\}\\Lambda\_\{k\}\(\\theta^\{\*\}\),\\quad\\bar\{H\}:=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}H\_\{k\}\.A first\-order M\-estimation expansion gives

θ¯−θ^N=−1K∑k=1K\(Hk−1−H¯−1\)Sk\+op\(N−1/2\)=Op\(N−1/2\)\.\\begin\{array\}\[\]\{ll\}\\bar\{\\theta\}\-\\hat\{\\theta\}\_\{N\}&=\-\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\left\(H\_\{k\}^\{\-1\}\-\\bar\{H\}^\{\-1\}\\right\)S\_\{k\}\+o\_\{p\}\(N^\{\-1/2\}\)\\\\\[10\.0pt\] &=O\_\{p\}\(N^\{\-1/2\}\)\.\\end\{array\}\(40\)SinceLNL\_\{N\}is stationary atθ^N\\hat\{\\theta\}\_\{N\}, the corresponding difference in loss is second order:

ΔN=LN​\(θ¯\)−LN​\(θ^N\)=Op​\(N−1\)\.\\Delta\_\{N\}=L\_\{N\}\(\\bar\{\\theta\}\)\-L\_\{N\}\(\\hat\{\\theta\}\_\{N\}\)=O\_\{p\}\(N^\{\-1\}\)\.\(41\)Equation \([40](https://arxiv.org/html/2610.10931#A1.E40)\) also shows why the finite\-sample discrepancy can be amplified when a local information matrix is poorly conditioned: the difference between averaged and centralized estimators involves the inverse local curvaturesHk−1H\_\{k\}^\{\-1\}\. This remains a finite\-sample effect and vanishes with increasingNN\. Combining the two regimes yields

ΔN=\{Op​\(N−1\),well specified, identifiable,Δ\+Op\(N−1/2\),when​Δ\>0\.\\Delta\_\{N\}=\\begin\{cases\}O\_\{p\}\(N^\{\-1\}\),&\\text\{well specified, identifiable\},\\\\\[4\.0pt\] \\Delta\+O\_\{p\}\(N^\{\-1/2\}\),&\\text\{when \}\\Delta\>0\.\\end\{cases\}\(42\)The separation is between a sampling discrepancy that vanishes with sample size and a selection\-relevant population effect that persists\.

*Surrogate informativeness\.*The conditional\-invariance argument above requires only that the partitioning score be a function ofxx\. This condition guarantees validity of the construction, but it does not guarantee that the induced environments are informative for DSCA\.

To see what makes a surrogate informative, consider a small perturbation around the pooled population law and write the environment weight as

K​πk​\(g⁡\(X\)\)=1\+ηk​\(X\)\.K\\pi\_\{k\}\(g\(X\)\)=1\+\\eta\_\{k\}\(X\)\.\(43\)Let

q⁡\(X\):=𝔼⁡\[∇θℓ​\(Z,θ∗\)∣X\]q\(X\):=\\mathbb\{E\}\\\!\\left\[\\nabla\_\{\\theta\}\\ell\(Z,\\theta^\{\*\}\)\\mid X\\right\]\(44\)denote the conditional population loss gradient at the pooled optimum\. A first\-order perturbation of the stationarity condition gives

θk∗−θ∗=−H−1​𝔼​\[ηk​\(X\)​q​\(X\)\]\+O⁡\(‖ηk‖2\),\\theta\_\{k\}^\{\*\}\-\\theta^\{\*\}=\-H^\{\-1\}\\mathbb\{E\}\\\!\\left\[\\eta\_\{k\}\(X\)q\(X\)\\right\]\+O\(\\\|\\eta\_\{k\}\\\|^\{2\}\),\(45\)whereH=∇2Λ​\(θ∗\)H=\\nabla^\{2\}\\Lambda\(\\theta^\{\*\}\)\. Hence the local optima move only to the extent that the covariate reweighting induced by the surrogate is aligned with regions in which the candidate has a nonzero conditional loss gradient\.

For a correctly specified model,q⁡\(X\)=0q\(X\)=0pointwise, consistently with Proposition[1](https://arxiv.org/html/2610.10931#Thmproposition1)\. For a misspecified model, a surrogate is informative when its ordering separates regions of covariate space in which the preferred calibration differs\. Conversely, if the surrogate\-induced allocation perturbations are approximately orthogonal toq⁡\(X\)q\(X\), the first\-order parameter displacement is small and DSCA provides little additional selection pressure\.

*Worked examples\.*LetY=X2Y=X^\{2\}withX∼Unif⁡\[0,1\]X\\sim\\mathrm\{Unif\}\[0,1\], and considerK=2K=2heterogeneous environments corresponding toX1=\[0,1/2\]X\_\{1\}=\[0,1/2\]andX2=\(1/2,1\]X\_\{2\}=\(1/2,1\]\. For the misspecified candidatef⁡\(x,θ\)=θ​xf\(x,\\theta\)=\\theta x,

θ1∗=38,θ2∗=4556,θ¯∗=3356≈0\.5893,θ∗=34\.\\theta\_\{1\}^\{\*\}=\\frac\{3\}\{8\},\\quad\\theta\_\{2\}^\{\*\}=\\frac\{45\}\{56\},\\quad\\bar\{\\theta\}^\{\*\}=\\frac\{33\}\{56\}\\approx 0\.5893,\\quad\\theta^\{\*\}=\\frac\{3\}\{4\}\.\(46\)Under squared\-error loss,

Λ⁡\(θ∗\)=180=0\.01250,Λ⁡\(θ¯∗\)=33115680≈0\.02111,Δ=273136≈0\.00861\.\\begin\{array\}\[\]\{lll\}\\Lambda\(\\theta^\{\*\}\)=\\frac\{1\}\{80\}=0\.01250,\\\\\[10\.0pt\] \\Lambda\(\\bar\{\\theta\}^\{\*\}\)=\\frac\{331\}\{15680\}\\approx 0\.02111,\\\\\[10\.0pt\] \\Delta=\\frac\{27\}\{3136\}\\approx 0\.00861\.\\end\{array\}\(47\)Thus DSCA raises the population loss assigned to this misspecified structure from0\.012500\.01250to approximately0\.021110\.02111, corresponding toΔ≈0\.00861\\Delta\\approx 0\.00861, or an increase of about69%69\\%relative to centralized calibration\. Under homogeneous environments the local and pooled optima coincide andΔ=0\\Delta=0\.

The second example shows why misspecification or disagreement alone is insufficient\. For the misspecified constant candidatef⁡\(x,θ\)=θf\(x,\\theta\)=\\thetaunder the same two environments,

θ1∗=112,θ2∗=712,θ¯∗=13=θ∗\.\\theta\_\{1\}^\{\*\}=\\frac\{1\}\{12\},\\qquad\\theta\_\{2\}^\{\*\}=\\frac\{7\}\{12\},\\qquad\\bar\{\\theta\}^\{\*\}=\\frac\{1\}\{3\}=\\theta^\{\*\}\.\(48\)Despite the substantial disagreement between the two local pseudo\-true parameters, their arithmetic mean coincides exactly with the centralized optimum, and thereforeΔ=0\\Delta=0\. HereH1=H2=2H\_\{1\}=H\_\{2\}=2, so this cancellation is precisely the behavior predicted by \([33](https://arxiv.org/html/2610.10931#A1.E33)\)\. DSCA introduces no population\-level selective penalty for this candidate\.

*Scope and limitations\.*The analysis above assumes identifiable population optima, smooth local risks, fixedKK, and an idealized calibration step based on exact empirical minimizers\. Symbolic simplification reduces exact algebraic non\-identifiability before calibration, but weak identification, near\-singularity, and numerical optimization error may remain in practice\. Strong heterogeneity can amplify these effects by reducing the information available within individual partitions, whereas the homogeneous limit suppresses the population DSCA effect\. These limitations motivate the empirical sensitivity analysis with respect toα\\alphaand the optimization controls reported in the main text\.

## Appendix BPer\-Function Results on the Synthetic Track

Figures[7](https://arxiv.org/html/2610.10931#A2.F7)and[8](https://arxiv.org/html/2610.10931#A2.F8)report the success rates on each synthetic function for the 100\-sample and 200\-sample settings, respectively\. Within each figure, the top row shows the low noise setting and the bottom row the high noise setting\. Each configuration has been tested with 10 different seeds, and an equation is considered recovered if it achievesR2\>0\.99R^\{2\}\>0\.99on the noiseless test set\. These results underlie the win/tie/loss counts and the aggregate counts of the main text\.

All evaluated methods struggle significantly withf3f\_\{3\}, frequently failing to identify the ground\-truth model\. Success rates remain at or near 0% across almost all configurations\. The only exceptions are DSCA and BFGS reaching 10% under BIC in the 200\-sample high\-noise regime\. This difficulty likely stems from the sine term, which usually represents a significant challenge for SR algorithms\.

Forf1f\_\{1\}andf2f\_\{2\}, DSCA generally outperforms BFGS and LM, particularly under the NLL and AIC selection criteria\. For instance, onf1f\_\{1\}in the 200\-sample low\-noise regime, DSCA achieves an 80% success rate using NLL, compared to 50% for BFGS and 30% for LM\. Performance is instead comparable across the three optimizers onf4f\_\{4\}andf5f\_\{5\}\. There, success rates typically cluster between 60% and 70% under low noise with 200 samples, and the few losses of DSCA reported in Table IV occur\.

These comparative trends are consistent across both noise regimes\. The higher noise setting yields lower overall recall, a degradation most pronounced forf4f\_\{4\}andf5f\_\{5\}\. For example, in the 200\-sample regime, high noise drops the success rates onf4f\_\{4\}from approximately 70% to between 0% and 20% for most criteria\.

Fig\. 7:Success rates on synthetic datasets \(N=100N=100\) under low and high noise regimes\.Fig\. 8:Success rates on synthetic datasets \(N=200N=200\) under low and high noise regimes\.