Progressive Subexpression Reuse in Symbolic Regression: Insights from RL-based Search and a Genetic Programming Realization

1 Introduction

Symbolic regression (SR) aims to infer explicit mathematical expressions from data, producing compact and interpretable models that are especially valuable in scientific discovery [Schmidt and Lipson, 2009; Udrescu and Tegmark, 2020]. In most SR settings, candidate models are represented as symbolic expression trees, and solving SR amounts to searching over a discrete and hierarchical space of such trees.

* Corresponding author: Rongye Shi (shirongye@buaa.edu.cn).

However, the space grows combinatorially with expression size and operator choices, making search inefficient and unstable, particularly on complex problems where many locally plausible structures exist.

Genetic programming (GP) is a classical framework for SR [Koza, 1994], evolving expression trees through selection, mutation, and crossover. Despite decades of progress, GPbased SR can struggle on harder tasks, where success often hinges on discovering useful building blocks (subexpressions) and reliably accumulating them so that increasingly complex structures can be composed. Many extensions strengthen the within-generation search dynamics by modifying selection and variation, including ϵ-lexicase selection [La Cava et al., 2016], semantic techniques [Virgolin et al., 2019], and linkage-learning-based recombination [Virgolin, 2020].

Recently, reinforcement learning (RL) has emerged as a strong alternative for SR [Petersen et al., 2019; Kamienny et al., 2022; Holt et al., 2023]. Beyond final performance, prior analyses suggest that RL-generated solutions often improve progressively during training, with parts of the target expression emerging over time [Du et al., 2024]. Motivated by this, we analyze expression samples during search and track the empirical frequencies of ground-truth subexpressions when the target expression is available. As shown in Fig. 1, RLbased SR (DSR [Petersen et al., 2019]) exhibits a clear hierarchical pattern: simpler subexpressions tend to emerge earlier and become frequent, followed by more complex ones, whereas standard GP shows more oscillatory and less ordered dynamics under the same evaluation budget. Additional results for other RL-based solvers and GP across different tasks are provided in the Appendix.

These observations motivate a mechanism-level hypothesis: effective SR benefits from progressive subexpression reuse, a closed-loop process in which promising subexpressions are (i) retained from current search outcomes and (ii) reintroduced to shape the next-stage sampling/initialization distribution, enabling stable accumulation of increasingly complex structures. This raises a natural question: can we realize such progressive reuse within GP, without policy training, by intervening only on stage-wise population construction?

We answer this question affirmatively and propose Reinforcement Genetic Programming (RGP), a fully gradientfree SR framework that keeps the GP inner loop unchanged and instead implements progressive reuse at the process level. RGP maintains a reusable subexpression pool that is refreshed from each stage's search outcomes and directly guides next-stage population construction via pool-guided initialization. In doing so, RGP offers a simple, trainingfree instantiation of the retention and reintroduction mechanism suggested by RL dynamics, and yields competitive performance against both GP baselines and learning-guided SR methods under matched evaluation budgets.

Our contributions are as follows:

• We propose a simple, reproducible protocol to track subexpression frequency trajectories during SR search, providing process-level evidence of a typically hierarchical emergence from simpler to more complex components in RL-based SR relative to GP baselines.
• We propose RGP, a process-level retention and reintroduction method that externalizes reusable subexpressions into a pool refreshed from current-stage outcomes and uses it to bias next-stage population construction, while keeping standard GP operators unchanged.
• We demonstrate on standard SR benchmarks that, under matched operator/constraint settings and budget-aligned evaluation, RGP matches or outperforms strong baselines, highlighting the practical value of explicit processlevel retention-reintroduction.

2 RelatedWork

Symbolic regression has been studied from multiple perspectives that differ in how symbolic expressions are generated and optimized. We focus on two lines most relevant to our work: GP-based symbolic regression and learningguided symbolic regression. Other paradigms, including physics-inspired decomposition, neural symbolic regression, pretrained sequence models, and tree search, have also been explored [Udrescu and Tegmark, 2020; Kim et al., 2020; Kamienny et al., 2022; Holt et al., 2023; Sun et al., 2022].

2.1 Genetic Programming for Symbolic Regression

Symbolic regression has traditionally been studied within the framework of genetic programming (GP), following its original formulation by Koza [Koza, 1994]. In GP-based approaches, candidate expressions are represented as syntax trees and evolved through selection and variation operators, making GP naturally suited to the discrete and hierarchical structure of symbolic expressions.

A longstanding challenge in GP-based symbolic regression is the identification and reuse of building blocks, i.e., subexpressions that contribute positively to overall performance. Most existing approaches address this issue implicitly. Pareto-based evolutionary strategies, such as AgeFitness Pareto optimization (AFP) [Schmidt and Lipson, 2010], introduce auxiliary objectives to promote diversity and mitigate premature convergence, which can indirectly increase the survival of useful subexpressions. Semantic methods further exploit behavioral information of programs. For example, ϵ-lexicase selection [La Cava et al., 2016] favors individuals that perform well on difficult subsets of training cases, while semantic backpropagation techniques [Virgolin et al., 2019] guide recombination based on desired intermediate outputs, encouraging the preservation of semantically meaningful subexpressions. Dependency-aware recombination methods, such as GP-GOMEA [Virgolin, 2020], explicitly model interactions between program components to better preserve cooperating substructures during crossover. Similar reusable-component ideas also appear beyond GP-based SR, especially in program synthesis systems that learn libraries of reusable programs or abstractions, such as DreamCoder [Ellis et al., 2021] and InnateCoder [Moraes et al., 2025].

2.2 Learning-guided Symbolic Regression

Reinforcement learning (RL) has emerged as a strong paradigm for symbolic regression by casting expression construction as a sequential decision-making problem. In this formulation, an agent generates an expression token by token according to a policy and receives a reward based on the data-fitting quality of the resulting expression. A representative example is Deep Symbolic Regression (DSR) [Petersen et al., 2019], which learns a sampling policy over symbolic token sequences via policy optimization, enabling effective exploration in large discrete spaces.

Beyond purely RL-driven generators, a growing line of work combines learned generation with evolutionary or search-based refinement. A common strategy is to use a learned generator or RL-guided search agent to shape the starting distribution of candidate expressions, and then apply GP-style variation and selection for downstream refinement. For example, GEGL [Ahn et al., 2020] uses a learned generator to seed the GP population, while NGGP [Mundhenk et al., 2021] couples evolutionary search with policy learning to bias generation toward high-reward expressions. RSRM [Xu and Sun, 2024] further leverages a reinforcementguided search agent to seed and steer symbolic refinement. These methods highlight the value of learned sampling or decision biases in symbolic search, while evolutionary operators provide complementary refinement and diversity.

2.3 Position of RGP

RGP is closest to GP reuse mechanisms, linkage-learning GP, and neural/RL-guided population seeding, but it introduces the reuse bias at a different control point. Standard elitism or Hall-of-Fame mechanisms retain complete high-fitness expressions as candidate solutions, whereas RGP extracts their subexpressions and turns them into explicit reusable initialization units. Linkage-learning methods such as GP-GOMEA use dependency information among tree components to guide recombination inside the GP loop; RGP does not learn linkage structure and does not modify selection, crossover, or mutation. Instead, it keeps the GP search loop standard and applies the retained subexpressions only before each stage, by biasing the construction of the next initial population. Neuralor RL-guided seeding methods learn a generator or policy to propose expressions, while RGP derives its bias directly from the evolving Hall-of-Fame-mined subexpression pool. Thus, RGP is not a new genetic operator or a learned proposal model, but an outer-loop retention-reintroduction mechanism for stage-wise population construction.

3 Preliminaries

3.1 Symbolic Trees

In symbolic regression, a mathematical expression is commonly represented as a symbolic tree. A symbolic tree is an ordered tree whose internal nodes correspond to operators (e.g., +, ×, sin), while leaf nodes represent variables or constants; each symbolic tree uniquely defines a mathematical expression, and vice versa.

3.2 Genetic Programming for Symbolic Regression

Genetic programming (GP) is a tree-based evolutionary search paradigm for symbolic regression. It maintains a population of symbolic expression trees initialized either randomly or from predefined primitives. At each generation, new candidate expressions are produced through subtree crossover and mutation, while selection is guided by data-fitting performance. To preserve high-quality solutions, GP-based methods often employ elitism mechanisms, such as a Hall-of-Fame, which retains top-performing expressions across generations.

3.3 Reinforcement Learning-based Symbolic Regression

Reinforcement learning-based symbolic regression formulates expression construction as a sequential decision-making process by generating a symbolic expression as a preorder traversal sequence s, enabling batch generation and optimization of expressions. Given a parameterized policy πθ, the training objective is to maximize the expected reward over sampled sequences:

where R(s) is defined based on the data-fitting quality of the symbolic expression τ(s) decoded from sequence s, such as the negative prediction error on the training data.

4 Method

This section presents Reinforcement Genetic Programming (RGP), a fully gradient-free symbolic regression framework that operationalizes the progressive subexpression reuse pattern observed in RL-based SR. RGP implements this idea via a stage-level retention-reintroduction loop: at each stage, we refresh a subexpression pool and use it to bias the initialization of the next population. Importantly, the GP inner loop remains unchanged-selection, crossover, mutation, and constant optimization are all standard-and our only intervention is stage-wise population construction, providing a minimal and attribution-friendly mechanism for reuse without policy training. We refresh the pool between stages (rather than every generation) so that recombination can proceed on a fixed initialization distribution within a stage, yielding more stable genetic search before the memory signal is updated.

4.1 The Framework of RGP

The core idea of RGP is to externalize stage-wise search feedback as reusable subexpressions mined from elite expressions and to reuse them to shape the next-stage initialization distribution, as illustrated in Fig. 2 and summarized in Algorithm 1. We use stage to denote the outer closed-loop iteration in Algorithm 1, while generation refers to the inner GP evolution within a stage.

RGP consists of three components that operate in a closed loop:

• Subexpression Pool. An explicit repository that stores reusable subexpressions mined from elite solutions and is refreshed across stages to balance exploitation and exploration.
• Pool-guided Population Construction. A procedure that constructs the initial GP population by incorporating subexpressions sampled from the pool into randomly generated expression trees.
• Genetic Programming Search. A standard GP process that evolves the constructed population to produce highquality symbolic expressions.

At stage t, we (i) refresh the subexpression pool from the hall-of-fame elites obtained up to stage t-1, and (ii) use the refreshed pool to construct a pool-guided initial population via injection. A standard GP run then evolves this population, and the newly found elites are fed back to refresh the pool for stage t+1. Through repeated stages, this closed-loop process progressively biases GP toward high-potential structures while retaining exploration.

4.2 Subexpression Pool

Fig. 1 suggests that RL-based SR may implicitly realize progressive subexpression reuse through its learned sampling distribution: once useful subexpressions appear in highreward expressions, they tend to reoccur more frequently in subsequent samples and later support the emergence of more complex structures. To examine this hypothesis, we seek a non-RL instantiation that implements the same retentionreintroduction loop within GP.

RGP instantiates progressive reuse via an explicit subexpression pool that biases population construction while being refreshed by GP search outcomes. Concretely, at stage t we (i) construct the pool using elites available up to stage t-1, and (ii) use the pool to construct the initial population for stage t. The pool merges two sources:

where S util t-1 collects reusable subexpressions mined from elites (exploitation), and S exp t-1 injects novel subexpressions to expand structural coverage and maintain diversity (exploration).

Exploitation via Promising Subexpressions

The exploitation component explicitly retains subexpressions that repeatedly occur in high-quality solutions. At the beginning of stage t, we update the reusable part of the pool by mining subexpressions from the global hall-of-fame Et-1. Conceptually, it proceeds in three steps: (i) mine candidates from elites, (ii) score each candidate, and (iii) select a compact reusable set.

Step 1: Mining candidates from elites. We maintain a global hall-of-fame set Et that stores the best top-K expressions found up to stage t, together with their fitness values. Each element τ ∈ Et is a symbolic expression tree with fitness r(τ ) (e.g., negative MSE), where a larger value indicates better fit. We decompose each τ ∈ Et-1 into its subtrees and collect syntactically valid subexpressions. We represent each subexpression s as a subtree, and define its size |s| as the number of tokens (nodes), implemented as the length of its preorder traversal sequence. To avoid overly trivial fragments and overly large structures, we keep only subexpressions with sizes |s| ∈ [la, ha]. Let Sb denote the deduplicated candidate set:

where Sub(τ ; [la, ha]) returns the set of subexpressions (subtrees) of τ whose sizes |s| lie in [la, ha].

Step 2: Scoring candidates by quality and recurrence. Not all candidates in Sb are equally useful. We summarize each s ∈ Sb by (i) its average fitness among elites that contain it and (ii) its occurrence frequency within the hall-of-fame:

where s ⪯ τ indicates that s is a subtree of τ, and I[·] is the indicator function.

Step 3: Selecting a reusable exploitation set. To form a compact exploitation memory, we rank candidates with a simple hierarchical rule: longer-first, then higher quality, then higher recurrence. We prioritize longer subexpressions because they often subsume shorter ones as subtrees and thus act as more informative, compositional building blocks. Moreover, longer subexpressions occur less frequently in elites, so their quality estimates can be noisier; a quality-first rule may prematurely discard potentially useful longer structures. We

Algorithm 1 RGP (overall loop)

Require: Dataset D; operator/grammar F; stages T; population size N; injection prob Pin; early-stop threshold ϵ; pool hyperparameters; hall-of-fame size K.

Ensure: Best expression τ ⋆ .

• 1: Notation: Et = global hall-of-fame (best-so-far expressions); Ht = history of used subexpressions; Pt = subexpression pool.
• 2: Initialize E0←∅, H0←∅, τ ⋆ ←∅
• 3: for t = 1 to T do
• 4: S util t-1 ← UpdateUtil(Et-1) // mine & select from HOF (Sec. 4.2)
• 5: S exp t-1 ←UpdateExp(Ht-1) // sample unseen valid (Sec. 4.2)
• 6: Pt ←S util t-1 ∪ S exp t-1 ; Ht ←Ht-1 ∪Pt
• 7: Sample base population {τ (0) i }Ni=1 ∼ RandInit(F)
• 8: {τ init i }Ni=1 ← Inject({τ (0) i },Pt, Pin) // subtree replacement (Sec. 4.3)
• 9: ˆ τt ← GP({τ init i }Ni=1 ,D) // best-of-stage (standard GP)
• 10: Et ← SelectTop(Et-1 ∪ {ˆ τt}; K, r)
• 11: τ ⋆ ← argmaxτ∈Et∪{τ ⋆ } r(τ)
• 12: if MSE(τ ⋆ ; D) ≤ ϵ then
• 13: break
• 14: end if
• 15: end for
• 16: return τ ⋆

therefore use fitness and recurrence only to break ties among similarly sized candidates.

Concretely, we sort s ∈ Sb by decreasing |s|, then decreasing ¯ r(s), and finally decreasing p(s). To reduce redundancy, we prune any candidate that is a strict subtree of a previously retained one under the same ranking. Finally, we keep the top-k candidates, denoted as

where ↓ denotes descending-order sorting. This set serves as the exploitation set for population construction at stage t.

Exploration via Structural Diversity

To mitigate premature convergence, RGP includes an explicit exploration component that promotes structural diversity in the subexpression pool. We maintain a lightweight history set Ht to record previously used subexpressions. Exploration candidates are drawn from the set of valid subexpressions, i.e., syntactically valid trees under the operator/grammar constraints with sizes |s| ∈ [lb, hb]. Since this range is bounded, the candidate set is finite; in practice, we choose a small hb and pre-enumerate all valid candidates with |s| ∈ [lb, hb] once before search into a cache.

At the end of stage t-1, we sample exploration subexpressions while excluding any candidate already in Ht-1:

After forming the pool at stage t, we record all subexpressions used in the pool:Ht ← Ht-1 ∪Pt. (8)

This history-based sampling encourages exploration of unseen structures and helps maintain diversity in the pool.

4.3 Population Construction Guided by the Subexpression Pool

Given a subexpression pool Pt, the remaining question is how to use this search feedback to improve symbolic search.

Recent hybrid symbolic regression methods [Mundhenk et al., 2021] often rely on a trained neural generator to propose candidate expressions as an initial population, followed by GP for downstream refinement. By shaping the starting distribution of candidates, such approaches can work synergistically with GP's evolutionary operators and improve the efficiency of subsequent search.

Following this insight, we use the subexpression pool to guide GP population construction at each stage. Specifically, at stage t, we first generate a base population {τ (0) i }Ni=1 using standard random tree initialization. We then perform poolguided injection with a fixed probability Pin per individual. For each base individual τ (0) i , we sample a Bernoulli gate zi ∼ Bernoulli(Pin). If zi = 1, we sample a small set (a fixed number in all experiments; see Appendix for details) of subexpressions S (i) in ⊂ Pt and apply a subtree-replacement operator that replaces randomly selected subtrees of τ (0) i with subexpressions from S (i) in to obtain a candidate ˜ τi; otherwise we keep the original individual. We accept the replacement only if the resulting tree satisfies standard GP validity constraints (syntactic consistency and a maximum size limit). If a sampled replacement is invalid, we resample locations/subexpressions up to a small fixed number of attempts; otherwise we keep the original tree. Formally,

The constructed population {τ init i }Ni=1 is then evolved by a standard GP run at stage t.

To further mitigate stagnation, RGP employs a lightweight restart rule: if the reusable exploitation set remains unchanged for a fixed number of consecutive stages, we clear the subexpression pool (optionally retaining the single best individual) and reinitialize the process, encouraging the search to escape local optima.

5 Experiments

In this section, we conduct comprehensive experiments to evaluate RGP on standard symbolic regression benchmarks, including comparative studies, ablation analyses, and efficiency evaluation.

5.1 EvaluationMetrics

We evaluate performance using the recovery rate [Petersen et al., 2019], defined as the fraction of independent runs that exactly recover the ground-truth expression. Compared to evaluations based on fixed datasets [Xu and Sun, 2024] or datafitting metrics such as MSE [Jiang et al., 2024], this metric provides a more faithful assessment of symbolic reconstruction capability. Formally, the recovery rate is defined as:

where Dj denotes a dataset randomly sampled from the ground-truth function f, A is the symbolic regression algorithm, J is the number of independent runs, and ≡ indicates symbolic equivalence between expressions. Unless stated otherwise, we compare methods under a matched evaluation budget measured by the total number of fitness evaluations (i.e., evaluated candidate expressions).

5.2 Baselines and Benchmarks

In our comparison, we include the following baseline methods in symbolic regression: GP, the classical genetic programming algorithm for symbolic regression; DSR, a policy gradient-based method that formulates expression generation as a sequential decision process [Petersen et al., 2019]; NGGP, which employs an RNN to initialize the GP population and iteratively updates the RNN using GP search results [Mundhenk et al., 2021]; DGSR, a transformer-based pretrained generative model for symbolic regression [Holt et al., 2023]; RSRM, which combines double Q-learning with Monte Carlo Tree Search to balance exploration and exploitation [Xu and Sun, 2024]; and UDFS, which performs unbiased symbolic regression by systematically searching the space of expression DAGs [Kahlmeyer et al., 2025].

For fair comparison, we use the same operator set and expression constraints (e.g., size/depth limits) across methods whenever applicable; method-specific settings follow the original implementations.

To assess the stability and efficiency of our approach, we evaluate it on three widely used symbolic regression benchmarks. These include: Nguyen, a widely used benchmark of symbolic expressions with one or two independent variables; R*, which consists of three rational functions with polynomial terms in both numerator and denominator, designed to increase structural complexity; and Livermore, a challenging set of equations involving high-degree exponentials, trigonometric components, and highly nonlinear polynomials.

5.3 Main Results

As shown in Table 1, under our matched operator/constraint settings and evaluation budgets, RGP achieves the highest (or tied-highest) recovery rates among the evaluated methods across the considered benchmarks. These improvements are consistent with our hypothesis that stage-wise retention and reintroduction of promising subexpressions can stabilize the accumulation of building blocks during GP search. RGP makes this mechanism explicit by mining reusable subexpressions from a global Hall-of-Fame of elites into a pool, which then guides population construction across stages, providing a more persistent reuse signal than relying only on withingeneration survival.

5.4 SRBench Results

We further evaluate RGP on SRBENCH [La Cava et al., 2021], which contains both white-box and black-box symbolic regression tasks. Specifically, we consider (i) the SRBENCH Feynman Equations suite (100 physics formulas with available ground-truth expressions) and (ii) the standard SRBENCH black-box tasks where the underlying functional forms are unknown. Throughout, we follow the official SRBENCH protocol, including the prescribed train/test splits and metric definitions. For the Feynman suite, we report the solution rate under the benchmark definition, i.e., the fraction of equations for which the recovered expression exactly matches the target up to symbolic equivalence, and we additionally apply the benchmark's additive Gaussian noise injection to the training targets under multiple relative noise levels to emulate measurement noise. For black-box tasks, where exact recovery is not defined, we report the SRBENCH test metric as predictive performance together with the corresponding model size (number of nodes), reflecting the accuracysimplicity tradeoff. On the Feynman data, RGP achieves superior solution rates compared to the baselines, while on the black-box tasks, it strikes a balance between predictive accuracy and model simplicity, achieving competitive R2 performance with relatively simpler expressions. Additional implementation details are provided in the Appendix.

5.5 Ablation

We conduct ablation studies on the Livermore benchmark to isolate the contributions of key components in our framework. Table 2 reports symbolic recovery rates averaged over multiple independent runs.

Model 1 (Full) corresponds to the complete framework, which maintains a subexpression pool with both exploitation (reusing promising subexpressions) and exploration (using structurally novel components), together with pool-guided subtree replacement during population initialization and standard GP search. Model 2 (w/o Exploitation) disables the reuse of promising subexpressions mined from elites, while retaining exploration and the replacement operation, isolating the role of exploitation in enabling progressive reuse. Model 3 (w/o Exploration) restricts the pool to promising subexpressions only, removing the use of novel components to assess how much diversity contributes beyond reuse. Model 4 (Random Pool) preserves the overall framework but replaces the pool contents with randomly sampled subexpressions matched in size, testing whether gains arise from meaningful reuse rather than the replacement operation itself. Model 5 (Pure GP) removes the subexpression pool entirely, reducing the method to standard GP with random restarts.

Overall, Model 3 (w/o Exploration) remains competitive with Model 1 (Full), whereas variants that disable Exploitation degrade substantially, indicating that retaining and reintroducing promising subexpressions is particularly important for RGP's performance.

5.6 Time Efficiency

We study the time-accuracy tradeoff of RGP by varying the initialization budget Binit, defined as the number of candidate expressions constructed for stage-wise population initialization (i.e., before GP evolution), while keeping the remaining GP inner-loop computation fixed (e.g., population size, number of generations, and fitness evaluations per stage). For each budget level, we run RGP on the full Nguyen benchmark and report the recovery rate averaged across all Nguyen tasks together with the corresponding wall-clock time measured under a fixed hardware/software setup (Appendix). As shown in Fig. 4, increasing Binit consistently improves recovery but exhibits diminishing returns at larger budgets, whereas runtime grows monotonically, quantifying a practical knob to trade compute for symbolic recovery.

6 Conclusion

We analyzed the training dynamics of RL-based symbolic regression and found empirical evidence consistent with progressive subexpression reuse, where useful subexpressions emerge early, become increasingly sampled, and later support the construction of more complex structures. Motivated by this mechanism-level perspective, we proposed Reinforcement Genetic Programming (RGP), a purely GP-based approach that explicitly retains and reintroduces promising subexpressions via a dynamically maintained pool, without relying on RL. Experiments on standard benchmarks show that RGP achieves competitive or superior performance with significantly reduced reliance on expensive policy training. These results suggest that modeling how useful subexpressions are retained and reused over search can be a major factor in effective symbolic regression.

Acknowledgments

This work was supported by the National Science and Technology Major Project (2022ZD0117801), the National Natural Science Foundation of China (62306023), and the Beijing Nova Program (20240484490).