SymMatika: Structure-Aware Symbolic Discovery

## Abstract

Abstract Symbolic regression (SR) seeks to recover closed-form mathematical expressions that describe observed data. While existing methods have advanced the discovery of either explicit mappings (i.e., y = f ​ ( 𝐱 ) y=f(\mathbf{x}) ) or discovering implicit relations (i.e., F ​ ( 𝐱 , y ) = 0 F(\mathbf{x},y)=0 ), few modern and accessible frameworks support both. Moreover, most approaches treat each expression candidate in isolation, without reusing recurring structural patterns that could accelerate search. We introduce SymMatika , a hybrid SR algorithm that combines multi-island genetic programming (GP) with a reusable motif library inspired by biological sequence analysis. SymMatika identifies high-impact substructures in top-performing candidates and reintroduces them to guide future generations. Additionally, it incorporates a feedback-driven evolutionary engine and supports both explicit and implicit relation discovery using implicit-derivative metrics. Across benchmarks, SymMatika achieves state-of-the-art recovery rates on the Nguyen and Feynman benchmark suites, an impressive recovery rate of 61% on Nguyen-12 compared to the next best 2%, and strong placement on the error-complexity Pareto fronts on the Feynman equations and on a subset of 57 SRBench Black-box problems. Our results demonstrate the power of structure-aware evolutionary search for scientific discovery. To support broader research in interpretable modeling and symbolic discovery, we have open-sourced the full SymMatika framework.

## 1 Introduction

In the late 16^th century, Tycho Brahe meticulously recorded the positions of celestial bodies. His data laid the foundation for Johannes Kepler, who in 1601 derived analytical expressions of the motion of these planets. These expressions launched a scientific revolution in their discovery that Mars’ orbit was in fact an ellipse. This was an early instance of symbolic regression (SR), where governing laws are distilled from observational data as interpretable mathematical expressions. The majority of SR methods similarly aim to discover symbolic expressions $f$ that relate variables in data $(\mathbf{X},y)$, with $X_{i}\in\mathbb{R}^{n}$, $y_{i}\in\mathbb{R}$, in the form of $y=f(\mathbf{x})$. Some others seek to uncover implicit relations $F(\mathbf{x},y)=0$ to discover invariants such as conservation laws.

Symbolic regression poses a computationally hard challenge, as the space of symbolic expressions grows exponentially with expression length Lu et al. (2016); Virgolin and Pissis (2022). The dominant approach, genetic programming (GP) Koza (1990), stochastically evolves populations of candidate expressions using mutation and crossover. Tools such as PySR Cranmer (2023), GPlearn Stephens et al. (2016), Operon Burlacu et al. (2020), and AFP Schmidt and Lipson (2010) apply GP to flexibly search across diverse mathematical structures without assuming parametric forms. However, their convergence is often slow due to uninformed exploration, generating many unpromising candidates across generations.

To improve efficiency, neural-guided SR approaches integrate deep learning to bias the search toward promising expressions. Methods such as DSR Petersen et al. (2019), AI Feynman Mundhenk et al. (2021), and NGGPPS Udrescu and Tegmark (2020) use RNNs or Transformers to guide expression generation or predict symbolic transformations. Deep learning and GP-based approaches are alternative methods of symbolic regression. We make comparisons between GP-based and deep-learning algorithms on benchmark performances later in Sec. [4](https://arxiv.org/html/2507.03110v2#S4).

Recent algorithms Petersen et al. (2019); Mundhenk et al. (2021) start to use recurrent neural networks (RNNs) to learn trends in expression coefficients. However, they do not capture higher-order syntactic structures, such as the recurring mathematical motifs found in top-performing individuals. Therefore, these approaches miss the opportunities to build complex expressions from learned substructures. Identifying and reusing these structures could accelerate discovery by enabling recombination of partial solutions into globally correct forms. Moreover, most existing SR systems specialize in either explicit mappings or implicit relations, but not both. The only known system supporting both, Eureqa Schmidt and Lipson (2009), has remained closed sourced and now is integrated into a commercialized platform.

In this work, we introduce SymMatika, a unified symbolic regression framework that discovers both explicit and implicit mathematical relations by combining multi-population feedback-based genetic programming with learned structural pattern reuse. Given data $\mathcal{D}_{E}=\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$ for explicit tasks or $\mathcal{D}_{I}=\{(\mathbf{x}_{i})\}_{i=1}^{N}$ for implicit ones, SymMatika iterates over three key phases: (1) feedback-based GP, which adaptively adjusts operation weights based on evolutionary context; (2) structural analysis, which extracts frequent syntactic substructures from high-performing expressions; and (3) population-specific updates to generation parameters, enabling structurally guided tree construction. The algorithm returns the fittest expression across all populations. In experiments on the Nguyen benchmark and Feynman equations, SymMatika outperforms state-of-the-art methods in accuracy and convergence. It also recovers implicit governing equations from experimental data in the Eureqa dataset with up to $100\times$ faster, demonstrating both versatility and performance.

## 2 Related Work

#### Genetic Programming for Symbolic Regression.

GP has long been a cornerstone of symbolic regression, originating with Koza’s foundational work Koza (1990, 1994), which introduced the idea of evolving expression trees via biologically inspired operations such as mutation and crossover. Modern GP-based SR systems such as Operon Burlacu et al. (2020) and PySR Cranmer (2023) inherit this lineage. Operon uses a steady-state GP model with tournament selection, while PySR employs a multi-population island model and simulated annealing to balance exploration and exploitation during evolution. These frameworks rely on fixed operator probabilities, which can hinder adaptive search in complex landscapes. Feedback-based adaptive crossover-rate in evolutionary computation Guan et al. (2024) proposes a modifiable crossover distribution to optimize crossover points, however it does not control the probability of crossover being selected over other genetic operations (i.e. mutation or single-node crossover). Self-adjusting mutation rates with provably optimal success rules Doerr et al. (2019) proposes updating mutation rates based on fitness-success of offspring, yet it does not consider the overall evolutionary progress of the algorithm.

SymMatika builds on this GP foundation but introduces two key innovations. First, it uses feedback-based operator scheduling that dynamically adjusts mutation, crossover, and selection rates based on each population’s recent evolutionary progress (Sec. [3.2](https://arxiv.org/html/2507.03110v2#S3.SS2)). Second, it augments the search with a reusable structural motif library, extracted from high-performing expressions across populations and generations (Sec. [3.3](https://arxiv.org/html/2507.03110v2#S3.SS3)). This enables the recombination of semantically meaningful substructures, promoting convergence toward globally correct solutions. Unlike previous GP models that operate solely at the token or subtree level with fixed heuristics, SymMatika evolves both structural and operator-level strategies in tandem.

#### Neural-Guided Symbolic Regression.

Neural-guided SR emerged to improve search efficiency by integrating neural networks into symbolic discovery Martius and Lampert (2016); Alaa and Van der Schaar (2019); Kamienny et al. (2022); Biggio et al. (2021); Champion et al. (2019). Among these, AI Feynman Udrescu and Tegmark (2020) is notable for using neural networks to detect symmetries, units, and separability in the data, enabling recursive decomposition of complex equations. More directly comparable to our method is DSO Mundhenk et al. (2021), which combines an RNN-based generator with stateless, random-restart GP loops. Their model generates $N$ candidate expressions per iteration, refines them over $S$ GP steps, and uses the best $M$ expressions to update the RNN via policy gradients.

While both our SymMatika and DSO use learning-based components to guide symbolic search, our approach diverges in two critical ways. First, instead of training a monolithic RNN to guide sampling, we perform population-specific frequency analysis over high-performing expressions to update generation parameters, enabling interpretable and efficient adaptation. Second, and more fundamentally, DSO learns only at the token level (e.g., which operators and constants are promising), whereas SymMatika extracts and reuses high-impact structural patterns – subtrees that recur in successful individuals. This is especially useful to recombine partial symbolic solutions into novel candidates for full solutions, facilitating both exploration and repair. Furthermore, our focus is orthogonal to recent work using deep neural networks to uncover latent variables Chen et al. (2022), equations Brunton et al. (2016), or structures Huang et al. (2024). These approaches aim to extract informative features from data. In fact, these approaches can integrate symbolic regression into their process of discovering governing principles.

#### Eureqa and Implicit Symbolic Regression.

Most modern SR systems focus on discovering explicit functional mappings $y=f(\mathbf{x})$, while implicit relations $F(\mathbf{x},y)=0$, commonly seen in physical systems governed by conservation laws or symmetries, remain underexplored. Eureqa Schmidt and Lipson (2009) is one of the few frameworks capable of discovering both explicit and implicit expressions. It uses a Pareto-optimized GP framework that perturbs, recombines, and simplifies expressions, scoring candidates via error and complexity. Importantly, Eureqa incorporates implicit-derivative metrics to recover nontrivial invariant equations.

Despite its versatility, Eureqa is over 16 years old and struggles on modern benchmarks due to fixed operator schedules and heuristics tailored for low-dimensional problems (originally $\leq 4$ variables). SymMatika preserves Eureqa’s implicit-derivative loss but augments it with improved GP techniques including feedback-driven operator tuning, structural motif reuse, and multi-population coordination. These enhancements enable SymMatika to recover implicit equations up to $100\times$ faster and achieves superior performance on explicit SR benchmarks such as Feynman and Nguyen.

## 3 SymMatika

SymMatika is composed of two core components: (1) a multi-population, feedback-driven genetic programming engine, and (2) a library of high-impact symbolic motifs that capture reusable substructures. In this section, we introduce the formal setup and describe each component, followed by their integration for discovering both explicit and implicit expressions in data.

### 3.1 Problem Setup

We represent symbolic expressions as algebraic trees $\tau$, with internal nodes as unary or binary operators (e.g., $+$, $\times$, $\cos$, $\log$) and leaf nodes as constants or variables. A pre-order traversal of $\tau$ yields the symbolic expression $f$, whose quality is assessed using task-specific fitness metrics.

#### Explicit relations.

Given data $\mathcal{D}_{E}=\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$ with $\mathbf{x}_{i}\in\mathbb{R}^{d}$, the goal is to find $f$ such that $y\approx f(\mathbf{x})$. We define fitness as mean-log-error (MLE):

$$ $\mathcal{L}_{\mathcal{D}_{E}}(f)=-\frac{1}{N}\sum_{i=1}^{N}\log\left(1+\left|y_{i}-f(\mathbf{x}_{i})\right|\right)$ $$

#### Implicit relations.

When target variables are not explicitly labeled, as is common in physical systems governed by invariants, we aim to discover expressions $f(\mathbf{x})=0$ that characterize underlying constraints or symmetries. Given a candidate expression $f(x_{1},x_{2},\dots,x_{n})=0$, to quantify fitness, we apply the implicit function theorem. Treating $x_{i}$ as an implicit function of $x_{j}$, the derivative is given by:

$$ $\frac{\partial x_{i}}{\partial x_{j}}=-\frac{\frac{\partial f}{\partial x_{j}}}{\frac{\partial f}{\partial x_{i}}}$ $$

We compute symbolic partial derivatives across all variable pairs and compare them against finite-difference numerical estimates. To account for variable interdependencies, which frequently occur in coupled dynamical systems, we generalize the paired partial derivative as:

$$ $\frac{\partial x_{i}}{\partial x_{j}}=\frac{\partial x_{i}+\partial x_{p}\cdot\frac{\Delta x_{p}}{\Delta x_{i}}}{\partial x_{j}+\partial x_{q}\cdot\frac{\Delta x_{q}}{\Delta x_{j}}}$ $$

where $x_{p}$ and $x_{q}$ are variables interdependent with $x_{i}$ and $x_{j}$, respectively. $\Delta$ represents finite differences in numerical analysis (including time-series data), or the numerical approximation of the derivative, and $\partial$ represents the calculus partial derivative of a symbolic expression. We evaluate all possible such pairings and take the worst-case pairing for evaluation to penalize expressions that perform well only under selective dependencies.

Let $M_{s}$ and $M_{n}$ denote the symbolic and numerical paired-partial derivative matrices (shape $\binom{d}{2}\times N$). The implicit fitness is:

$$ $\mathcal{L}_{\mathcal{D}_{I}}(f)=-\frac{1}{N}\sum_{i=1}^{N}\log\left(1+\left\|M_{n}(\mathbf{x}_{i})-M_{s}(\mathbf{x}_{i})\right\|\right)$ $$

### 3.2 Feedback-Based Genetic programming

Tree-based genetic programming uses two core operators: crossover, which swaps subtrees at selected nodes between parents, and mutation, which perturbs nodes within an individual. These are typically coupled with a selection mechanism (e.g., tournament selection) that favors high-fitness candidates. While effective, traditional GP applies static operator rates and rigid selection rules, often overlooking valuable substructures in lower-fitness individuals and failing to adapt to population dynamics. To address these limitations, we extend GP with four key mechanisms:

#### Single-node crossover.

Traditional subtree crossover introduces large structural changes. To enable finer control during late-stage optimization, we introduce single-node crossover, which swaps only one same-type node between trees (i.e. binary operation $\leftrightarrow$ binary operation, unary operation $\leftrightarrow$ unary operation, variable $\leftrightarrow$ variable, etc.).

#### Temperature-guided selection and mutation.

Simulated annealing is a well-established strategy for balancing exploration and exploitation in evolutionary algorithms. PySR Cranmer (2023) incorporates this by rejecting a mutation with probability $p=\exp\left(\frac{L_{F}-L_{E}}{\alpha T}\right)$, where $L_{E}$ and $L_{F}$ are the fitness scores of an expression before and after mutation, $T\in[0,1]$ is the annealing temperature, and $\alpha$ is a scaling hyperparameter. While effective for mutation acceptance, PySR continues to rely on tournament selection, which prioritizes only the fittest individuals, favoring exploitation over exploration.

To promote further exploration of potentially promising but lower-fitness individuals, we extend simulated annealing to the selection mechanism itself. During selection, we randomly sample a subset of individuals from a population and assign each a Boltzmann probability: $p_{f}=\exp(\frac{\mathcal{L}_{\mathcal{D}}(f)}{T})$ where $\mathcal{L}_{\mathcal{D}}(f)$ is the fitness of candidate $f$. Selection is then performed via roulette sampling over normalized probabilities. At high temperatures ($T\to 1$), selection approximates uniform sampling, promoting exploration. As $T$ decreases, the probability mass concentrates on high-fitness candidates, effectively converging toward tournament-style selection. This adaptive strategy allows the selection mechanism to gradually shift from exploration to exploitation as the evolutionary process progresses.

We apply a similar temperature-dependent strategy to mutation. At high temperatures, coarse mutations such as subtree replacement are favored to encourage diversity. As $T$ lowers, finer mutations, such as constant perturbations or operator swaps, are more likely. Each population $P_{i}$ maintains a temperature-adjusted distribution over mutation types, enabling population-specific tuning of structural granularity. Together, these temperature-guided mechanisms provide principled control over both the scope of variation and the selective pressure during evolution.

#### Feedback-based operator scheduling.

Traditional GP systems use fixed rates for genetic operators such as mutation and crossover Schmidt and Lipson (2009). While simple, static rates are suboptimal: coarse-grained changes like subtree crossover are useful early in evolution, whereas fine-grained adjustments like coefficient tuning are better suited for later stages. To enable adaptive behavior, we introduce feedback-based scheduling that dynamically adjusts operator probabilities based on the evolutionary context of each population.

Let $g_{0}$ be the average fitness of the top-$M$ individuals in a population prior to GP loop, and let $g_{n}$ be the same statistic after $n$ generations. Let $h$ denote the number of consecutive generations with negligible fitness improvement (plateau), defined as $|g_{n}-g_{n-1}|<\epsilon$ for some threshold $\epsilon=1e^{-6}$. We define the operator probability function:

$$ $\mathbb{P}(g_{n},h)=\begin{cases}m_{i}\pm(\frac{|g_{0}-g_{n}|}{g_{0}}\cdot|m_{f}-m_{i}|)\pm 0.1h&h\geq 2\\ m_{i}\pm(\frac{|g_{0}-g_{n}|}{g_{0}}\cdot|m_{f}-m_{i}|)&\text{else}\end{cases}$ $$

Here, $m_{i}$ and $m_{f}$ are the initial and final operator probabilities for each genetic operation (e.g. crossover, mutation), with signs determined by whether each operator frequency should increase or decrease over time. Crossover starts with probability 60% and decreases to 5%, single-node crossover starts with probability 10% and increases to 15%, mutation starts with probability 30% and increases to 80%. This function is safe-guarded by these probability bounds and our probabilistic selection model, so it is resilient to abnormal degradations in average fitness.

#### Island populations with migration.

We evolve a set of populations or “islands” in parallel. The island model for GP is a powerful tool for simultaneously evolving multiple distinct evolutionary paths Duarte et al. (2017); Whitley et al. (1999). We allow for small migrations between islands by swapping small subsets of individuals between populations. Initially, 1% of island populations are swapped every 20 iterations. These subsets increase to 2% of island populations according to a growing control rate $\alpha_{M}$. We keep migrations small and infrequent to promote population diversity and prevent convergent evolution, which we tend to observe when migration rates exceed 3%.

### 3.3 Motif Library and Structure Reuse

Figure: Figure 2: Motif-based recombination across evolving populations with shared top performers.
Refer to caption: Diagrams/MotifLibrary.jpg

After each iteration of GP, we generate new candidates to replace subsets in each population $P$ of worst-performing individuals. Candidate generation in SymMatika is modeled as a parameterized distribution over expression trees, denoted $p(\tau\mid\theta)$, where $\tau$ is a symbolic expression and $\theta$ encodes node-level generation probabilities (e.g., for selecting constants, variables, and operators) Mundhenk et al. (2021); Petersen et al. (2019). These parameters are updated based on frequency feedback from high-fitness individuals: for each operator type $\theta_{i}$, we compute its frequency $f_{\theta_{i}}$ across top-performing expressions, and apply the update:

$$ $\theta_{i}\leftarrow\theta_{i}+\beta_{\theta_{i}}f_{\theta_{i}}$ $$

Here, $\beta_{\theta_{i}}$ is a scaling factor that accounts for natural biases, e.g., binary operators ($\beta_{\theta_{i}}=0.2$) such as $+$ and $\times$ may occur more often than unary functions ($\beta_{\theta_{i}}=1.0$) like $\log$ or $\exp$, and are thus adjusted more conservatively. These biases were selected from observations and tuning of our parameterized candidate generator during initial testing. This feedback mechanism enables population-specific adaptation of the generative distribution, which can bias expression construction toward promising token-level patterns.

However, frequency-based feedback alone is limited: it captures which tokens are useful, but cannot be used to explicitly discover structural patterns or partial solutions in individuals. For instance, in a symbolic regression problem with correct solution $f=x^{2}+y^{2}+z^{2}$, a partial solution $f^{\prime}=x^{2}+y^{2}-\cos(z)$ would not be recognized as having some correct terms through frequency analysis. To better identify and reuse partial solutions, we introduce a structural motif discovery mechanism inspired by biological sequence motifs.

In molecular biology, a sequence motif is a recurring nucleotide or amino-acid pattern with functional significance Tateno et al. (1997); Liu et al. (2002); Chou and Schwartz (2011); Grant and Bailey (2021). Analogously, we define a symbolic motif as a high-impact subexpression, a subtree within an expression that contributes strongly to its overall fitness. We use this intuition to construct a reusable library of symbolic motifs and employ them for guided recombination.

#### Initialization and Data Structure.

To implement structural reuse, we maintain a central motif library $\mathbb{M}_{f}$ that stores high-impact symbolic subexpressions across populations. This library is structured as a $d\times k$ table, where $d$ corresponds to the number of data variables and each row holds up to $k$ motifs associated with a particular variable, where $k$ is typically a small integer ($\leq 20$). Each entry in the table consists of a motif $\mathcal{M}_{v}=(\tau\textquoteright,I(\tau\textquoteright))$, where $\tau\textquoteright$ is a symbolic subtree and $I(\tau\textquoteright))$ is its estimated impact score on fitness. We calculate impact score with the following equation $I(\tau\textquoteright))=\mathcal{L}(\tau)-\mathcal{L}(\tau-\tau\textquoteright)$, where $\mathcal{L}(\tau)$ is the loss of the original tree expression and $\mathcal{L}(\tau-\tau\textquoteright)$ is the loss of the original tree without the subtree $\tau\textquoteright$. This quantifies the importance of the substructure to the overall quality of the expression.

#### Motif Generation.

After each generation, we examine the top $M$ individuals from each island population to identify candidate motifs. For a given expression $\tau$ (e.g., $\tau=3x^{2}-\cos(y)$), we extract the left and right subtrees of each internal node (e.g., $\tau_{L}=3x^{2}$, $\tau_{R}=-\cos(y)$. Note: $\cos(y)$ is a negative term in the example, so we extract it as $-\cos(y)$) and assess their fitness contribution.

Each motif is then associated with a specific input variable $v$ by scanning the subtree in pre-order and assigning it to the first variable encountered. This heuristic is based on the assumption that symbolic components are typically rooted in a dominant variable, and indexing motifs by variable helps maintain coverage and diversity during recombination. The newly extracted motifs are compared to the existing entries in row $v$ of $\mathbb{M}_{f}$. If a new motif has a higher impact than the lowest-ranked entry, it replaces it. If not and the row has not been filled, the new motif will be added. The motifs in each row are then re-sorted by their impact scores to retain the most promising candidates.

#### Structure-Aware Expression Synthesis.

Once the motif library is updated, new expressions are synthesized by sampling one or more motifs from each row of $\mathbb{M}_{f}$ until all variables are represented. These motifs are then assembled into full expressions using randomly selected binary algebraic operators such as $+$, $-$, $\times$, or $\div$, resulting in composite symbolic expressions $\tau\textquoteright$. The resulting expressions are inserted into a dedicated motif population $P_{\mathcal{M}}$, which co-evolves alongside the main island populations.

#### Co-Evolve Motif-Population and Main Populations.

To propagate promising structural components throughout the evolutionary process, high-performing expressions from $P_{\mathcal{M}}$ are periodically injected back into the main populations. This allows strong partial solutions—discovered independently across islands and generations—to be recombined and reused in novel ways, ultimately improving both convergence speed and expression quality.

This structure-aware mechanism complements token-level parameter updates by identifying semantically meaningful subtrees for reuse and recombination. It accelerates symbolic discovery by enabling partial solution reuse and coordinated variable coverage. Motif recombination also improves robustness to local optima by synthesizing expressions from independently validated high-impact components. We visualize this pipeline in Fig. [2](https://arxiv.org/html/2507.03110v2#S3.F2), outline the full algorithm in Alg. [1](https://arxiv.org/html/2507.03110v2#alg1), visualize the full model in Fig. [1](https://arxiv.org/html/2507.03110v2#S0.F1), and list hyperparameters in Appendix Tab. [A.3](https://arxiv.org/html/2507.03110v2#Ax1.SS3).

### 3.4 Implementation Details

SymMatika is a multi-core C++23 library, compiled with clang++ and optimized using -O3 and -march=native flags. It leverages SymEngine for symbolic computation, Eigen for linear algebra and paired-partial derivative operations, and OpenMP for parallelization. The full framework is open-sourced to support future work in interpretable modeling and symbolic discovery.

## 4 Experiments

For all experiments, SymMatika is instantiated with population sizes of 10,000 individuals per island, followed by truncation to the top 400 for the GP-loop. The number of islands $I$ is set proportional to expression complexity as $I\propto 2d$, where $d$ is the number of variables (capped at $I=8$ for large SR-problems with five or more variables due to computational considerations). We run each experiment for 1500 generations of GP.

We report: (1) recovery rates on the Nguyen benchmark suite (with an emphasis on Nguyen-12 recoveries), (2) error-complexity plots of the Feynman equations (mean proportion of $R^{2}>0.99$ vs. equation complexity) and 57 SRBench Black-box problems (median $R^{2}$ vs. equation complexity), and (3) implicit-relation discovery in physical systems using Eureqa datasets. We also performed ablation studies to assess the contribution of core components of SymMatika. Our goals are to demonstrate the effectiveness of SymMatika on both explicit and implicit SR tasks, highlight cases where it outperforms prior work, and understand how different algorithm modules contribute to performance.

All experiments are run on a 2023 Apple M3 Max MacBook Pro with a 14-core CPU, 30-core GPU, and 36GB unified memory. Unlike other neural-guided approaches, SymMatika does not require GPUs nor access to large language models (LLMs). Since SymMatika uses multi-core execution, and many baselines are either single-core or multi-core with unknown runtimes, we report recovery rates (not wall-clock time) for Nguyen and Feynman.

### 4.1 Nguyen Benchmark

**Table 1: Exact recovery rates (%) on the Nguyen benchmark over 100 runs. SymMatika achieves the highest average performance (96.5%) and recovers Nguyen-12 at significantly higher rates than existing algorithms (61% success).**
| Dataset | Expression | SymMatika | NGGPPS | DSR | PQT | Eureqa | Operon | PySR |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Nguyen-1 | $x^{3}+x^{2}+x$ | 100 | 100 | 100 | 100 | 100 | — | — |
| Nguyen-2 | $x^{4}+x^{3}+x^{2}+x$ | 100 | 100 | 100 | 99 | 100 | — | — |
| Nguyen-3 | $x^{5}+x^{4}+x^{3}+x^{2}+x$ | 100 | 100 | 100 | 86 | 95 | — | — |
| Nguyen-4 | $x^{6}+x^{5}+x^{4}+x^{3}+x^{2}+x$ | 100 | 100 | 100 | 93 | 70 | — | — |
| Nguyen-5 | $\sin(x^{2})\cos(x)-1$ | 100 | 100 | 72 | 73 | 73 | — | — |
| Nguyen-6 | $\sin(x)+\sin(x+x^{2})$ | 100 | 100 | 100 | 98 | 100 | — | — |
| Nguyen-7 | $\log(x+1)+\log(x^{2}+1)$ | 98 | 97 | 35 | 41 | 85 | — | — |
| Nguyen-8 | $\sqrt{x}$ | 100 | 100 | 96 | 21 | 0 | — | — |
| Nguyen-9 | $\sin(x)+\sin(y^{2})$ | 100 | 100 | 100 | 100 | 100 | — | — |
| Nguyen-10 | $2\,\sin(x)\cos(y)$ | 100 | 100 | 100 | 91 | 64 | — | — |
| Nguyen-11 | $x^{y}$ | 100 | 100 | 100 | 100 | 100 | — | — |
| Nguyen-12 | $x^{4}-x^{3}+\tfrac{1}{2}y^{2}-y$ | 61 | 0 | 0 | 0 | 0 | 2 | 0 |
| Average |  | 96.5 | 91.4 | 83.6 | 75.2 | 73.9 | — | — |

Figure: Figure 3: Error–complexity (rank) trade-off on the Black-box (left) and Feynman (right) datasets.
Refer to caption: Diagrams/feynman_blackbox.jpg

Tab. [1](https://arxiv.org/html/2507.03110v2#S4.T1) reports recovery rates over 100 independent runs per task. We take reported results from Petersen et al. (2019); Mundhenk et al. (2021) on NGGPPS (Note: PQT is one method of training NGGPPS’ RNN), DSR, and Eureqa. SymMatika records an average recovery rate of 96.5%, outperforming all other methods by statistically significant margins (maximum $p<10^{-5}$ from $z$-score testing). We ran 100 runs with random seeds on Nguyen-12 using PySR and Operon with the exact experiment settings (details in Appendix [A.3](https://arxiv.org/html/2507.03110v2#Ax1.SS3)) and reported 0/100 and 2/100 recoveries respectively; given the simpler nature of Nguyen 1-11, we felt it was not necessary to perform testing on problems proven to be easily solvable. Notably NGGPPS Mundhenk et al. (2021) reports a 12% success rate on a relaxed variant with expanded data range (Nguyen-$12^{\star}$), our performance remains significantly stronger on the original Nguyen-12 dataset. These results highlight the advantage of combining structural motif reuse with feedback-driven parameter tuning, especially for long and composite expressions.

### 4.2 Feynman Equations Benchmark

**Table 2: Exact recovery counts on the Feynman benchmark (100 tasks). SymMatika recovers 73 expressions, outperforming all other algorithms.**
| GPlearn | AFP-FE | DSR | uDSR | AIFeynman | PySR | LaSR | SymMatika |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 20/100 | 26/100 | 23/100 | 40/100 | 38/100 | 59/100 | 72/100 | 73/100 |

Next, we evaluate on the Feynman benchmark, a widely adopted dataset of 100 symbolic physics expressions derived from the Feynman Lectures. As in prior work, expressions are sampled to form datasets $\mathcal{D}=\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$, and success is defined as recovery of the symbolic target. We performed two sets of experiments: (1) recovery rates on 100 Feynman problems outlined in LaSR Grayeli et al. (2024), and (2) error-complexity plot (mean $R^{2}>0.99$ proportion vs. equation complexity) against SRBench La Cava et al. (2021) results. We use the same experimental settings as for the Nguyen benchmark (Appendix [A.3](https://arxiv.org/html/2507.03110v2#Ax1.SS3)).

As shown in Tab. [2](https://arxiv.org/html/2507.03110v2#S4.T2), SymMatika recovers 73/100 expressions outperforming all methods, including GPlearn Stephens et al. (2016), DSR Petersen et al. (2019), uDSR Landajuela et al. (2022), AIFeynman Udrescu and Tegmark (2020), PySR Cranmer (2023), and LaSR Grayeli et al. (2024). LaSR Grayeli et al. (2024) extends PySR Cranmer (2023) with natural language priors from LLMs via a ConceptAbstraction function, enabling semantic guidance during search. While this approach is orthogonal to ours, our results show that structurally guided recombination and operator adaptation alone are competitive with language-informed strategies.

To showcase some sample output, we run SymMatika on a sample equation from the Feynman Lectures on Physics Feynman et al. (2015), the relativistic Doppler shift formula $\omega=\frac{1+\frac{v}{c}}{\sqrt{1-\frac{v^{2}}{c^{2}}}}\omega_{0}$. After running SymMatika on its dataset, it returns the expression $\omega=\sqrt{\frac{c+v}{\left(c-v\right)}}\omega_{0}$ with an MLE loss of $-3.10862\times 10^{-16}$, and the discovered expression is mathematically equivalent to the ground truth expression, albeit syntactically slightly different. We provide outputs on a subset of sampled Feynman Equations in Appendix Sec. [A.1](https://arxiv.org/html/2507.03110v2#Ax1.SS1).

For the error-complexity plot, SymMatika places impressively on the Pareto front. We report mean $R^{2}>0.99$ proportion of 0.942 and complexity of 13.294 across 6 trials (Fig. [3](https://arxiv.org/html/2507.03110v2#S4.F3)). We surpass all algorithms on SRBench – including Operon, AIFeynman and SBG-GP to name a few — with maximum $p=0.044$ — except for MRGP Arnaldo et al. (2014), which reports mean proportion 0.931 yet complexity 26316, significantly larger than SymMatika. On the Black-box subset, however, we report much better results in both error and complexity then MRGP. Additionally, two versions of TPSR Shojaee et al. (2023) slightly outperform our mean ($R^{2}>0.99$) accuracy with 0.949 (TPSR[$\lambda=0.1$]) and 0.952 (TPSR[$\lambda=0$]) to our 0.942, however our complexity is significantly improved with our 13.294 compared to their 84.42 (TPSR[$\lambda=0$]) and 57.22 (TPSR[$\lambda=0.1$]). Overall, SymMatika places strongly on the Pareto front and leads in model complexity.

### 4.3 Black-box Benchmark

We evaluated SymMatika on 57 Black-box problems gathered from SRBench used by TPSR Shojaee et al. (2023) which include synthetic and real-world datasets and we recorded median $R^{2}$ and mean equation complexity for the error-complexity plot (Fig. [3](https://arxiv.org/html/2507.03110v2#S4.F3)). We ran 5 full trials and recorded median $R^{2}$ of 0.931 and equation complexity of 28.059. We outperform nearly all benchmarked algorithms in median $R^{2}$ and match TPSR [$\lambda=0.5$] (complexity 82.58), where $\lambda$ is a hyperparameter balancing fitting and complexity reward), and falling slightly behind TPSR [$\lambda=0$] at 0.938 (complexity 129.85), Operon at 0.937 (complexity 61.74) and TPSR [$\lambda=0.1$] at 0.945 (complexity 95.71). Our complexity significantly improves on these algorithms by a factor of 3-4.5 for TPSR and a factor of 2 for Operon, and MRGP which reports a much-lower median $R^{2}$ of 0.502 larger and complexity 9878.172. TPSR records one trial so we cannot report meaningful $p$-testing given limited experimental data, yet we can report that the Operon’s outperformance of us is not statistically significant ($p=0.71138$).

### 4.4 Eureqa Physical Systems

To evaluate implicit symbolic regression, we compared against Eureqa on four time-series physical systems from its supplementary materials, including circular, pendulum, spherical, and double pendulum motions. Unlike the explicit SR benchmarks, these systems do not provide labeled outputs $y$, and recovery involves discovering hidden algebraic constraints $f(\mathbf{x})=0$ from the raw dynamics.

**Table 3: Runtime (s) on Eureqa dataset. SymMatika converges $10\times$–$100\times$ faster.**
| System | SymMatika (s) | Eureqa (s) |
| --- | --- | --- |
| circle_1 | 0.15 | 25 |
| pendulum_h_1 | 0.98 | 31 |
| sphere_1 | 3.00 | 320 |
| dbl_pendulum_h_1 | 900.00 | 27000 |

SymMatika used the same experimental setup as with previous benchmarks. As Eureqa supports 32-core execution with distributed infrastructure, we report wall-clock convergence times to compare implicit SR performance (Tab. [3](https://arxiv.org/html/2507.03110v2#S4.T3)). SymMatika recovers the ground-truth implicit expressions outlined in the Eureqa S.O.M. at very improved speeds. We observe $\sim$100$\times$ speedups for simpler systems like circle_1, and 10$\times$ for more complex systems like double_pendulum_h_1. We believe that these gains stem not only from modern hardware but also from algorithmic improvements, including adaptive operator scheduling, and motif-driven structural reuse.

### 4.5 Ablation Study

**Table 4: Recovery rates for individual benchmark problems with model ablations across 20 independent trials per problem. 95% confidence intervals are obtained from the recovery rates across all 12 Nguyen problems for each ablation.**
|  | Recovery rate (%) |  |  |  |
| --- | --- | --- | --- | --- |
| Problem | Baseline | $\theta$ Generation | Motif Library | Full Model |
| Nguyen-1 | 100 | 100 | 100 | 100 |
| Nguyen-2 | 100 | 100 | 100 | 100 |
| Nguyen-3 | 100 | 100 | 100 | 100 |
| Nguyen-4 | 90 | 95 | 100 | 100 |
| Nguyen-5 | 95 | 100 | 100 | 100 |
| Nguyen-6 | 85 | 90 | 95 | 100 |
| Nguyen-7 | 60 | 80 | 85 | 95 |
| Nguyen-8 | 100 | 100 | 100 | 100 |
| Nguyen-9 | 100 | 100 | 100 | 100 |
| Nguyen-10 | 100 | 100 | 100 | 100 |
| Nguyen-11 | 100 | 100 | 100 | 100 |
| Nguyen-12 | 5 | 20 | 40 | 65 |
| Nguyen average | 86.3 | 90.4 | 93.3 | 96.7 |
| 95% confidence intervals | $\pm$15.92 | $\pm$13.02 | $\pm$9.82 | $\pm$5.7 |

To assess the contribution of core components of SymMatika, we perform ablations on the Nguyen benchmark under four configurations: (1) a baseline with standard GP (no parameter tuning or motif reuse), (2) baseline + $\theta$-based parameterized candidate generation, (3) baseline + structural motif library, and (4) the full model with both mechanisms enabled. Each configuration is evaluated across 20 independent runs per task as shown in Tab. [4](https://arxiv.org/html/2507.03110v2#S4.T4).

The full model consistently achieves high recovery rates, with improvements observed on nearly every task. Parameterized generation provides consistent boosts across all benchmarks, while motif reuse is especially helpful on longer expressions (e.g., the challenging Nguyen-12 Sun et al. (2025)). In isolation, the motif mechanism raises the Nguyen-12 recovery rate from 5% to 40%, and further improves to 65% when combined with $\theta$-based parameter adaptation. These results confirm that token-level trends and reusable substructures capture complementary inductive biases that are critical for recovering long or irregular expressions.

## 5 Conclusions

We presented SymMatika, a symbolic regression framework that combines feedback-guided genetic programming with a reusable structural motif library to recover both explicit and implicit expressions from data. Through adaptive operator scheduling, motif-based recombination, and implicit-derivative fitness evaluation, SymMatika outperforms state-of-the-art methods on the Nguyen and Feynman benchmarks and uniquely recovers complex and implicit equations, including Nguyen-12 and Eureqa’s physical systems. These results highlight the value of integrating structural priors into symbolic search. Future work will explore robustness to noise and extend the framework with a user-friendly UI to support broader adoption in scientific discovery.

## 6 Limitations

The primary limitation in our approach is our results on the Feynman Equations. Although we report results outperforming or performing on par with leading models, we are still unable to solve a subset of 27 problems. These problems generally involve long expressions with complicated structures, such as $n=\frac{n_{0}}{\exp(\mu_{m}B/(k_{b}T))+\exp(-\mu_{m}B/(k_{b}T))}$, problem $\text{II}.35.18$ in the Feynman Equations Udrescu and Tegmark (2020). Other unsolved problems include large arguments in unary operators (e.g. $x=\sqrt{x_{1}^{2}+x_{2}^{2}-2x_{1}x_{2}\cos(\theta_{1}-\theta_{2})}$). To attempt further discovery of Feynman equations, we will conduct experiments in an improved hardware setup with more CPU cores to evolve more island populations in parallel with larger populations and for longer. Additionally, we will further investigate methods of discovering implicit relations in higher-dimensional data (i.e. $\geq 5$ variables).

## 7 Acknowledgments

This work is supported by ARO under award W911NF2410405, DARPA FoundSci program under award HR00112490372, DARPA TIAMAT program under award HR00112490419, and by gift supports from BMW and OpenAI.

## References

- Alaa and Van der Schaar [2019]
A. M. Alaa and M. Van der Schaar.
Demystifying black-box models with symbolic metamodels.
*Advances in neural information processing systems*, 32, 2019.
- Arnaldo et al. [2014]
I. Arnaldo, K. Krawiec, and U.-M. O’Reilly.
Multiple regression genetic programming.
In *Proceedings of the 2014 annual conference on genetic and evolutionary computation*, pages 879–886, 2014.
- Biggio et al. [2021]
L. Biggio, T. Bendinelli, A. Neitz, A. Lucchi, and G. Parascandolo.
Neural symbolic regression that scales.
In *International Conference on Machine Learning*, pages 936–945. Pmlr, 2021.
- Brunton et al. [2016]
S. L. Brunton, J. L. Proctor, and J. N. Kutz.
Discovering governing equations from data by sparse identification of nonlinear dynamical systems.
*Proceedings of the national academy of sciences*, 113(15):3932–3937, 2016.
- Burlacu et al. [2020]
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*, pages 1562–1570, 2020.
- Champion et al. [2019]
K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton.
Data-driven discovery of coordinates and governing equations.
*Proceedings of the National Academy of Sciences*, 116(45):22445–22451, 2019.
- Chen et al. [2022]
B. Chen, K. Huang, S. Raghupathi, I. Chandratreya, Q. Du, and H. Lipson.
Automated discovery of fundamental variables hidden in experimental data.
*Nature Computational Science*, 2(7):433–442, 2022.
- Chou and Schwartz [2011]
M. F. Chou and D. Schwartz.
Biological sequence motif discovery using motif-x.
*Current protocols in bioinformatics*, 35(1):13–15, 2011.
- Cranmer [2023]
M. Cranmer.
Interpretable machine learning for science with pysr and symbolicregression. jl.
*arXiv preprint arXiv:2305.01582*, 2023.
- Doerr et al. [2019]
B. Doerr, C. Doerr, and J. Lengler.
Self-adjusting mutation rates with provably optimal success rules.
In *Proceedings of the Genetic and Evolutionary Computation Conference*, pages 1479–1487, 2019.
- Duarte et al. [2017]
G. Duarte, A. Lemonge, and L. Goliatt.
A dynamic migration policy to the island model.
In *2017 IEEE congress on evolutionary computation (CEC)*, pages 1135–1142. IEEE, 2017.
- Feynman et al. [2015]
R. P. Feynman, R. B. Leighton, and M. Sands.
*The Feynman lectures on physics, Vol. II: The new millennium edition: mainly electromagnetism and matter*, volume 2.
Basic books, 2015.
- Grant and Bailey [2021]
C. E. Grant and T. L. Bailey.
Xstreme: Comprehensive motif analysis of biological sequence datasets.
*BioRxiv*, pages 2021–09, 2021.
- Grayeli et al. [2024]
A. Grayeli, A. Sehgal, O. Costilla Reyes, M. Cranmer, and S. Chaudhuri.
Symbolic regression with a learned concept library.
*Advances in Neural Information Processing Systems*, 37:44678–44709, 2024.
- Guan et al. [2024]
X. Guan, T. Yang, C. Zhao, and Y. Zhou.
Feedback-based adaptive crossover-rate in evolutionary computation.
In *Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence*, pages 6923–6930, 2024.
- Huang et al. [2024]
K. Huang, D. H. Cho, and B. Chen.
Automated discovery of continuous dynamics from videos.
*arXiv preprint arXiv:2410.11894*, 2024.
- Kamienny et al. [2022]
P.-A. Kamienny, S. d’Ascoli, G. Lample, and F. Charton.
End-to-end symbolic regression with transformers.
*Advances in Neural Information Processing Systems*, 35:10269–10281, 2022.
- Koza [1990]
J. R. Koza.
*Genetic programming: A paradigm for genetically breeding populations of computer programs to solve problems*, volume 34.
Stanford University, Department of Computer Science Stanford, CA, 1990.
- Koza [1994]
J. R. Koza.
Genetic programming as a means for programming computers by natural selection.
*Statistics and computing*, 4:87–112, 1994.
- La Cava et al. [2021]
W. La Cava, B. Burlacu, M. Virgolin, M. Kommenda, P. Orzechowski, F. O. de França, Y. Jin, and J. H. Moore.
Contemporary symbolic regression methods and their relative performance.
*Advances in neural information processing systems*, 2021(DB1):1, 2021.
- Landajuela et al. [2022]
M. Landajuela, C. S. Lee, J. Yang, R. Glatt, C. P. Santiago, I. Aravena, T. Mundhenk, G. Mulcahy, and B. K. Petersen.
A unified framework for deep symbolic regression.
*Advances in Neural Information Processing Systems*, 35:33985–33998, 2022.
- Liu et al. [2002]
J. S. Liu, M. Gupta, X. Liu, L. Mayerhofere, and C. E. Lawrence.
Statistical models for biological sequence motif discovery.
In *Case Studies in Bayesian Statistics: Volume VI*, pages 3–32. Springer, 2002.
- Lu et al. [2016]
Q. Lu, J. Ren, and Z. Wang.
Using genetic programming with prior formula knowledge to solve symbolic regression problem.
*Computational intelligence and neuroscience*, 2016(1):1021378, 2016.
- Martius and Lampert [2016]
G. Martius and C. H. Lampert.
Extrapolation and learning equations.
*arXiv preprint arXiv:1610.02995*, 2016.
- Mundhenk et al. [2021]
T. N. Mundhenk, M. Landajuela, R. Glatt, C. P. Santiago, D. M. Faissol, and B. K. Petersen.
Symbolic regression via neural-guided genetic programming population seeding.
*arXiv preprint arXiv:2111.00053*, 2021.
- Petersen et al. [2019]
B. K. Petersen, M. Landajuela, T. N. Mundhenk, C. P. Santiago, S. K. Kim, and J. T. Kim.
Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients.
*arXiv preprint arXiv:1912.04871*, 2019.
- Schmidt and Lipson [2009]
M. Schmidt and H. Lipson.
Distilling free-form natural laws from experimental data.
*science*, 324(5923):81–85, 2009.
- Schmidt and Lipson [2010]
M. D. Schmidt and H. Lipson.
Age-fitness pareto optimization.
In *Proceedings of the 12th annual conference on Genetic and evolutionary computation*, pages 543–544, 2010.
- Shojaee et al. [2023]
P. Shojaee, K. Meidani, A. Barati Farimani, and C. Reddy.
Transformer-based planning for symbolic regression.
*Advances in Neural Information Processing Systems*, 36:45907–45919, 2023.
- Stephens et al. [2016]
T. Stephens et al.
Genetic programming in python, with a scikit-learn inspired api: gplearn.
*Documentation at https://gplearn. readthedocs. io/en/stable/intro. html*, 2016.
- Sun et al. [2025]
C. Sun, S. Shen, W. Tao, D. Xue, and Z. Zhou.
Noise-resilient symbolic regression with dynamic gating reinforcement learning.
*arXiv preprint arXiv:2501.01085*, 2025.
- Tateno et al. [1997]
Y. Tateno, K. Ikeo, T. Imanishi, H. Watanabe, T. Endo, Y. Yamaguchi, Y. Suzuki, K. Takahashi, K. Tsunoyama, M. Kawai, et al.
Evolutionary motif and its biological and structural significance.
*Journal of molecular evolution*, 44:S38–S43, 1997.
- Udrescu and Tegmark [2020]
S.-M. Udrescu and M. Tegmark.
Ai feynman: A physics-inspired method for symbolic regression.
*Science advances*, 6(16):eaay2631, 2020.
- Virgolin and Pissis [2022]
M. Virgolin and S. P. Pissis.
Symbolic regression is np-hard.
*arXiv preprint arXiv:2207.01018*, 2022.
- Whitley et al. [1999]
D. Whitley, S. Rana, and R. B. Heckendorn.
The island model genetic algorithm: On separability, population size and convergence.
*Journal of computing and information technology*, 7(1):33–47, 1999.

## Appendix

### A.1 Sample Outputs From Subset of Feynman Equations

**Table 5: Comparison of ground-truth vs. discovered equations.**
| Equation Number | Ground Truth Equation | Discovered Equation |
| --- | --- | --- |
| I.18.4 | $r=\frac{m_{1}r_{1}+m_{2}r_{2}}{m_{1}+m_{2}}$ | $r=\frac{r_{1}m_{1}}{\left(m_{1}+m_{2}\right)}+\frac{r_{2}m_{2}}{\left(m_{1}+m_{2}\right)}$ |
| II.11.20 | $P_{*}=\frac{0.333333333333333E_{f}p_{d}^{2}n_{\rho}}{Tk_{b}}$ | $P_{*}=\frac{n_{\rho}p_{d}^{2}E_{f}}{3k_{b}T}$ |
| I.37.4 | $I_{*}=I_{1}+I_{2}+2\sqrt{I_{1}I_{2}}\cos{\left(\delta\right)}$ | $I_{*}=I_{1}+I_{2}+2\sqrt{I_{1}I_{2}}\cos{\left(\delta\right)}$ |
| I.24.6 | $E_{n}=\frac{1}{4}m(\omega^{2}+\omega_{0}^{2})x^{2}$ | $E_{n}=\frac{1}{2}m(\omega^{2}+\omega_{0}^{2})\frac{1}{2}x^{2}$ |
| III.17.37 | $f=\beta(1+\alpha\cos(\theta))$ | $f=\beta+\beta\cos{\left(\theta\right)}\alpha$ |
| I.27.6 | $f_{f}=\frac{d_{2}}{\left(n+\frac{d_{2}}{d_{1}}\right)}$ | $f_{f}=\frac{1}{\frac{1}{d_{1}}+\frac{n}{d_{2}}}$ |
| I.47.23 | $c=\sqrt{\frac{\gamma pr}{\rho}}$ | $c=\sqrt{\frac{\gamma pr}{\rho}}$ |
| I.12.11 | $F=q(E_{f}+Bv\sin(\theta))$ | $F=qE_{f}+qBv\sin(\theta)$ |

We demonstrate outputs on a subset of the Feynman Equations and compare ground truth equations to discovered equations. The ground truth equations and discovered equations are syntactially different although mathematically equivalent. These discoveries were found with implementation details consistent with all other experiments. All discovered equations record an MLE loss of $\leq 10^{-16}$.

### A.2 SymMatika Algorithm

Figure: Algorithm 1 SymMatika Algorithm

### A.3 Tabulated SymMatika Parameters

**Table 6: Experimental settings of SymMatika.**
| Parameters | Explanation |
| --- | --- |
| $G=1500$ | number of generations in evolutionary algorithm. |
| $m$ | number of distinct island populations (see Section [4](https://arxiv.org/html/2507.03110v2#S4)). |
| $k=200$ | number of iterations of GP per generation (see Fig [1](https://arxiv.org/html/2507.03110v2#S0.F1)). |
| $t=0,1$ | type of SR problem, i.e. explicit or implicit relation (see Algorithm [1](https://arxiv.org/html/2507.03110v2#alg1)). |
| $\alpha_{M}=0.02$ | migration rate of individuals between island populations (see Algorithm [1](https://arxiv.org/html/2507.03110v2#alg1)). |
| $\beta_{\theta}=(\underbrace{0.2,\dots,0.2}_{5\ \text{times}},\underbrace{1.0,\dots,1.0}_{13\ \text{times}})$ | natural bias term for node type $\theta_{i}$ (see Section [3.3](https://arxiv.org/html/2507.03110v2#S3.SS3)). |
| $m_{i}=0.60,0.10,0.30$ | initial genetic operator probabilities (crossover, single-node crossover, mutation) |
| $m_{f}=0.05,0.15,0.80$ | final genetic operator probabilities (crossover, single-node crossover, mutation) |

Tab. [6](https://arxiv.org/html/2507.03110v2#Ax1.T6) is a list of SymMatika’s experimental settings. We performed each experiment with random seeds.

### A.4 SymMatika GUI

With the open-sourced release of our codebase, we also provide a user-friendly GUI (Fig. [4](https://arxiv.org/html/2507.03110v2#Ax1.F4) and Fig. [5](https://arxiv.org/html/2507.03110v2#Ax1.F5)) which allows for easy and intuitive symbolic discovery for users who prefer graphical interface interactions with SymMatika.

Figure: Figure 4: SymMatika GUI open screen
Refer to caption: Diagrams/SymMatikaOpenScreen.png

Users should navigate through each header top-to-bottom beginning with Data and add the necessary information. In each header, there is also an instruction popup button to help explain the content for each section.

Figure: Figure 5: SymMatika GUI screen after solving Nguyen-5
Refer to caption: Diagrams/SymMatikaEndScreen.png

Installation and setup instructions are available on our GitHub: [https://github.com/generalroboticslab/SymMatika](https://github.com/generalroboticslab/SymMatika).

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: What specific improvements to the algorithm or hardware setup do you believe could help tackle the subset of unsolved Feynman equations? ​ Are there plans to incorporate additional techniques, such as deep learning or domain-specific priors, to address these challenges? 

How does SYMMATIKA's structural motif reuse compare to neural-guided symbolic regression methods in terms of scalability and interpretability? Could the framework benefit from integrating neural network components for further optimization?


While SYMMATIKA demonstrates impressive speed improvements over Eureqa, how does it compare in terms of computational resource requirements with other modern SR frameworks like PySR or TPSR?
soundness: 3
strengths: 1. SYMMATIKA supports both explicit and implicit relation discovery, making it applicable to a wide range of scientific problems.

2.The introduction of a motif library enables the identification and recombination of high-impact substructures, accelerating convergence and improving robustness to local optima.
summary: This work proposed a novel symbolic regression framework-SYMMATIKA that combines feedback-driven genetic programming with a reusable structural motif library to discover both explicit and implicit mathematical expressions from data. By leveraging adaptive operator scheduling, motif-based recombination, and implicit-derivative fitness evaluation, SYMMATIKA achieves state-of-the-art recovery rates on standard benchmarks like Nguyen and Feynman equations, significantly outperforming existing methods in accuracy, convergence speed, and model complexity.
weaknesses: However, the framework has limitations, including its inability to solve certain complex Feynman equations and a lack of experiments on real-world applications. 

While the results are promising, some comparisons lack statistical significance due to limited experimental data.  

Furthermore, the novelty of SYMMATIKA is somewhat constrained as it builds upon  implicit-derivative metrics, and its structural motif reuse, while innovative, is inspired by established concepts from biological sequence analysis.

---

confidence: 5
contribution: 3
initial_rating: 6
presentation: 3
questions: Please address the questions in the weakness section, thanks.
soundness: 3
strengths: Conceptual novelty in motif library that measures impact via ablation and re-injects high-impact subexpressions.

Suitable benchmark datasets and benchmark algorithms are evaluated, with appropriate ablation study (albeit a very short one on a single benchmark, i.e., only Nguyen).

Tackles both explicit and implicit relationships, in which the latter is impactful but has not been addressed much in existing literature.

Organization and flow of paper is well-designed.
summary: The work introduces SymMatika, a symbolic regression algorithm that integrates multi-island genetic programming with a motif library for structural reuse and feedback-driven operator scheduling.
weaknesses: It is hard to see what the gap in performance for SymMatika with its adjacent algorithms. For example, in Fig. 3, it seems that Operon is better than SymMatika in terms of R^2, but maybe the difference is only less than 0.01 R^2. It is not possible for the reader to determine the difference as it is reported now. [1] has also proved that adding or removing algorithms that are not on the Pareto front can paradoxically cause the set of Pareto-optimal algorithms to change when aggregating ranks. Thus, I would recommend to analyze Pareto-optimality by using the actual metric, instead of the rank of the metric.

[1] Fong, Kei Sen, and Mehul Motani. "Pareto-Optimal Fronts for Benchmarking Symbolic Regression Algorithms.", In ICML’25

No confidence intervals or error bars, but there is an indication of z-score testing. Can the paper describe the statistical test procedure in more details? I had assumed Wilcoxon ranked sign test would be the de facto test, but I would like to know more details about the full testing procedure used. Also, can the paper include the results of SymMatika against the other methods as well, I can only find the results for Operon vs SymMatika. This is easy to resolve, so in the recommendation score, I am assuming that the confidence intervals/error bars and statistical test results for the other comparisons will be added.

I can’t seem to find SNIP [2] in the left of Figure 3. Why is this so? 

[2] Meidani, Kazem, et al. "SNIP: Bridging mathematical symbolic and numeric realms with unified pre-training.", In ICLR’24

I also cannot find the work [2] cited in the references list, despite being used as a baseline in the right of Figure 3. I have not checked for all the other algorithms, SNIP only came to mind because it was used inconsistently, as mentioned above. Can the paper make sure that all algorithms used in the paper are cited?

Please provide a clear, self-contained, reproducible definition of equation "recovery" that is used in this paper. If constant matching is involved, please explain how 2 constants are determined to be "equal".

---

confidence: 3
contribution: 3
initial_rating: 4
presentation: 2
questions: See the weaknesses.
soundness: 3
strengths: 1.Novel integration of motif-level structural reuse in symbolic regression.  

2.Feedback-driven adaptive operator scheduling for efficient search.  

3.Supports both explicit and implicit relation discovery.  

4.Strong empirical results across multiple standard benchmarks.  

5.No reliance on deep neural networks or GPUs, making it computationally efficient and accessible.

6.Unified symbolic regression framework for both explicit and implicit relations.
summary: This paper proposes SYMMATIKA, a structure-aware symbolic regression (SR) framework that can discover both explicit relations and implicit relations.  SYMMATIKA integrates two main innovations: 1) Feedback-driven multi-island genetic programming, which adaptively tunes mutation, crossover, and selection rates based on evolutionary progress. 2) A reusable motif library, inspired by biological sequence motifs, to identify and reuse high-impact symbolic subexpressions for faster convergence. The method achieves state-of-the-art recovery rates on Nguyen, Feynman, SRBench, and Eureqa benchmarks.  Notably, it recovers 96.5% of Nguyen tasks (including 61% success on Nguyen-12, compared to 2% for prior methods), and converges 10–100× faster than Eureqa on implicit physical systems.  The system is fully implemented in optimized C++ and requires only multi-core CPUs.
weaknesses: 1. **Limited theoretical analysis of convergence and motif importance.**  

   The paper does not provide a formal convergence argument for the co-evolution of motifs and populations (Sec. 3.3–3.4). Motif impact $I(\tau') = L(\tau) - L(\tau - \tau')$ is introduced heuristically (line 302 – 310) without theoretical justification that this local fitness differential leads to global improvement.  While empirically effective, there is no analysis of stability (e.g., whether motif reuse can cause premature convergence). Related to convergence theory in GP (e.g., Whitley 1999, Doerr 2019), SYMMATIKA’s dynamics remain descriptive rather than proven.

2. **Evaluation lacks noisy and high-dimensional datasets (≥5 variables).**  

   Experiments (Sec. 4, Tables 4–5) primarily involve low-dimensional benchmarks (Nguyen: ≤2 vars; Feynman: ≤4 vars).  The authors explicitly note in Appendix A.5 (L791 – L809) that performance deteriorates for equations with ≥5 variables, acknowledging scalability issues.  No robustness study under observational noise was performed, which limits the method’s applicability to real scientific data.

3. **Some Feynman equations remain unsolved; scalability is limited.**  

   Appendix A.5 reports 27 unsolved Feynman problems, particularly those with long nested unary operators or composite exponentials (e.g., problem II.35.18).  The authors attribute this to limited population size and CPU-only hardware, but it also indicates that motif reuse and feedback adaptation do not fully handle high-complexity functional structures.  This limitation suggests that the system still struggles with extremely deep symbolic trees.

4. **Results mostly deterministic; stochastic robustness not studied.**  

   All experiments (Sec. 4.1–4.4) are run with fixed hyperparameters and deterministic seeds.  There is no variance or sensitivity analysis across random seeds, mutation rates, or population sizes.  As a result, the general stability of convergence and reproducibility under stochastic perturbations remain uncertain.

---

confidence: 4
contribution: 2
initial_rating: 4
presentation: 3
questions: 1. The paper combines several established techniques (e.g., Boltzmann selection, multi-island evolution, reusable library search). What is the theoretical foundation that unifies these techniques into a novel methodology?
2. The experiments for implicit regression focus on four test functions and one baseline, while related work discussion omits recent advances in implicit regression. How does the proposed method explicitly improve upon  recent implicit regression techniques, and what justifies the choice of the baseline comparison algorithm?
soundness: 2
strengths: The proposed method can simultaneously support both explicit and implicit regression tasks, and demonstrates superior performance over existing approaches across multiple benchmark tests.
summary: This paper proposes SYMMATIKA, a symbolic regression framework that combines multi-island genetic programming and a reusable symbol library to accelerate search, supporting both explicit (y=f(x)) and implicit (F(x,y)=0) regression tasks. Experimental results demonstrate its superiority over existing methods on benchmarks, particularly in recovery rate and computational efficiency.
weaknesses: This paper appears to be more of a combination and fine-tuning of existing proposed methods.
For instance, the paper introduces a Boltzmann selection mechanism, which has already been in use in the Evolutionary Algorithms  field for over 20 years. The multi-island evolution mechanism is also a very classic approach. The paper proposes a feedback-based operator scheduling mechanism, but the design of this mechanism seems ad-hoc, lacking sufficient explanation and theoretical depth. Search mechanisms based on a reusable library have also been proposed by scholars in the Genetic Programming field long ago.
Furthermore, the review of relevant literature is insufficient. For example, the discussion on Implicit Regression only involves a single paper from 2009, overlooking new methods proposed in recent years. For implicit regression, the experimental section only analyzes the method on four simple test functions and is compared against only one other algorithm, making the results unconvincing.