SMILE: Bridging Continuous Optimization and Discrete Symbolic Recovery

Conference: NeurIPS2026
arXiv: 2609.04639
Area: Interpretability
Keywords: symbolic regression, structural decomposition, symbolic activations, gated pruning, constant recovery

TL;DR

SMILE identifies decomposable structure in data, fits a network with fixed symbolic activations, and recovers compact expressions through pruning, constant refitting, and gradient-based rounding, achieving strong symbolic recovery under high noise in SRBench without universally leading noise-free recovery or prediction accuracy.

Background & Motivation

Symbolic regression (SR) seeks not only predictions but also readable, executable mathematical expressions. Genetic programming searches expression trees directly, but faces a combinatorial space involving operators, connections, and constants; neural networks fit data readily through gradient-based optimization, yet successful fitting does not imply recovery of a short, correct equation. When several simple factors are multiplied or added, their combined function can be difficult to learn, and the resulting expression can retain many floating-point coefficients.

Embedding mathematical operations such as sine and exponential into a network provides one bridge between these approaches. Continuous training can nevertheless retain too many paths, yielding a network that can be expanded into a formula rather than an informative scientific law. SMILE therefore treats data as more than a supervision signal: it first examines which variables approximately follow power laws and which disrupt that pattern, then selects direct fitting, reciprocal-target fitting, or separate fitting of two subexpressions. This reduces the compositional burden of an individual fit rather than merely reducing input dimensionality.

The paper separates the roles of training and recovery: training finds a numerically useful representation, while recovery compresses its structure and checks which constants can be simplified. Core Idea: use compositional cues in the data to constrain the fitting task of a symbolic network, then progressively convert continuous parameters into a compact discrete expression under constraints on prediction changes.

Method

Overall Architecture

The input consists of numerical sample pairs, and the output is a closed-form expression rather than a predictor that retains all network parameters. The pipeline proceeds through Variable Analysis and Decomposition, Symbolic Activation Network, Gated Pruning and Refitting, and Gradient-Based Constant Rounding. In the decomposition case, it separately fits a subexpression involving the remaining variables and a residual relationship involving the selected variable, then combines and recovers them.

%%{init: {'flowchart': {'rankSpacing': 24, 'nodeSpacing': 28, 'padding': 6, 'wrappingWidth': 400}}}%%
flowchart TD
    A["Numerical samples"] --> B["Variable Analysis<br/>and Decomposition"]
    B -->|Direct, reciprocal, or subproblems| C["Symbolic Activation<br/>Network"]
    A -.->|Fitting supervision| C
    C --> D["Gated Pruning<br/>and Refitting"]
    D --> E["Gradient-Based<br/>Constant Rounding"]
    E --> F["Final expression evaluation"]
    X["New input"] --> F

The dashed edge denotes supervision from the training data; at inference time, the recovered expression receives new inputs without rerunning structural analysis or network training. The reciprocal branch must invert the recovered result, while the decomposition branch must retain its additive or multiplicative combination.

Key Designs

1. Variable Analysis and Decomposition: identify the variable that complicates an otherwise simple relationship

For each input variable, the method independently performs log-linear regression between the absolute target and the absolute variable, then compares goodness of fit. If all variables exhibit similarly high goodness of fit, the pipeline favors direct training. It also compares low-degree polynomial fits to the target and its reciprocal; if the reciprocal is easier to fit, it learns that form first to reduce numerical instability associated with denominators. If one variable has substantially lower goodness of fit, the method treats it as a possible additional compositional factor rather than immediately increasing network depth.

Decomposition cannot obtain many records with an exactly equal variable value from continuously sampled data. The authors use windows around two fixation points, restrict the other variables to narrow bands, and estimate target sensitivity to the selected variable through local linear regression. The following constancy score is the operational criterion in the appendix: the local regression slope estimates the partial derivative, window width determines total drift, and the denominator is the target standard deviation within the windowed slab:

\[ C=\frac{|\beta|W}{\mathrm{std}(y_{\mathrm{slab}})}<\tau_{\mathrm{window}}. \]

Windows are tested from largest to smallest, preferring one that contains more samples and passes the threshold. The defaults are 6 candidate relative half-widths from 20% to 2%, with a threshold of 0.05. If all fail, the method still falls back to the smallest window: fixing a variable remains approximate and does not receive an error guarantee for every problem. It first recovers the relationship among the remaining variables within a window, then checks whether multiplicative residuals around the two fixation points exhibit an approximately constant ratio. If so, it fits division residuals; otherwise, it defaults to subtraction residuals, and finally multiplies or adds the recovered components.

This design suits cases where one variable disrupts an otherwise approximate power-law pattern, but low log-linear goodness of fit does not prove separability. Zero inputs or targets also require numerical handling in the logarithmic scan. Multiple coupled non-power-law variables, near-zero residual denominators, and too few usable samples within a window can all weaken decomposition decisions.

2. Symbolic Activation Network: let continuous parameters control directly expandable mathematical operations

Each hidden layer contains five fixed activations: sine, multiplication, identity, logarithm, and exponential. Sine, identity, logarithm, and exponential operate on linear combinations of inputs; the multiplication neuron learns an exponent for each input in logarithmic space. A single unit can therefore represent powers, roots, and ratios without searching operators individually. Fully dense residual connections link all layers, and an appended constant 1 supplies a bias path.

\[ m=\exp\!\left(\sum_j w_j\log z_j\right)=\prod_j z_j^{w_j}. \]

This real-valued identity requires positive logarithm inputs; it does not establish safe computation for arbitrary signed inputs and arbitrary real exponents. The implementation clamps logarithm inputs from below and exponential inputs from above, while penalizing inputs in clamped regions. This controls numerical behavior during training but changes the function within those regions, and the clipping boundaries are not everywhere smooth. Interpreting the network as an ordinary mathematical expression still requires checking its valid domain and whether clamping affected the fit.

The network jointly learns edge gates and activation gates: each edge weight is multiplied by its gate's sigmoid value, and each neuron output is further scaled by its activation gate. All gate logits start at 0, corresponding to an opening of 0.5. A sparsity penalty favors closing unnecessary paths, while another gate penalty discourages intermediate openings in preparation for discrete removal.

The authors construct sigmoid functions with identity, exponential, and logarithm, then emulate sufficiently wide MLPs by stacking and replicating blocks in parallel to argue universal approximation for the network family. This is an expressivity argument with sufficient capacity, not a guarantee that the experimental shallow networks with five neurons per layer and 1–2 layers represent arbitrary functions, nor a guarantee that gradient-based training discovers the true symbolic equation.

3. Gated Pruning and Refitting: remove structure before recalibrating the retained coefficients

Learned gates are not already a completed discrete structure selection. During recovery, the method tentatively closes individual gates, selects the candidate with the smallest effect on training goodness of fit, and accepts closure if the drop remains below a tolerance. It repeats until no further removal is possible. Edge gates remove connections, activation gates remove entire neurons, and only the surviving subgraph is composed layer by layer into an expression.

Transferring training weights unchanged into the reduced expression is unreliable: removing one path changes the optimal values of other coefficients. The method therefore refits numerical constants in the extracted expression by least squares, providing a better starting point for rounding. Pruning uses local removal costs on the training set, and the resulting subgraph is a greedy outcome rather than a globally minimal expression. Removal order can matter when paths compensate for one another.

4. Gradient-Based Constant Rounding: simplify according to output sensitivity rather than decimal places alone

A coefficient close to an integer is not necessarily safe to round: it may occur inside an exponential or a high-frequency sine, strongly affecting the output. Conversely, a larger coefficient change may be acceptable on an insensitive path. SMILE substitutes a candidate simple constant, estimates the resulting output change, and requires it to remain below a threshold across all checked samples:

\[ |c_i-r_i|\left|\frac{\partial\hat f}{\partial c_i}(r_i,\mathbf{x})\right|<\tau_{\mathrm{round}},\qquad \forall\mathbf{x}\in S. \]

The appendix first uses the mean value theorem to derive an error bound controlled by the maximum derivative along the parameter-change interval. To reduce computation, the implemented procedure evaluates the derivative only at the candidate rounded value as an approximation to that maximum. These are different claims: an endpoint-based first-order local sensitivity check does not guarantee the interval-wide error bound, nor prove that the chosen integer, rational, or familiar constant is the true constant that generated the data. It establishes passage of an approximate check, not a proof of scientific discovery.

A Worked Example

For II.15.4 in Appendix Table 3, the target is \(-\mu B\cos(\theta)\) and the recovered expression is \(-\mu B\sin(\pi/2-\theta)\). Combining multiplicative variables with an angle variable illustrates why requiring every variable to follow the same power law is insufficient. A decomposition route can first fit the multiplicative part at a fixed angle, then learn the angle-dependent residual.

This example illustrates the mechanism rather than claiming that the paper records this problem's actual internal route. Recovery permits a sine-based equivalent instead of the original cosine form, so a correct expression need not have the same syntax tree. This problem succeeds at all four tested noise levels in the table, but that does not establish recovery of arbitrary trigonometric equations.

Loss & Training

The training loss combines data mean squared error, a clamping-region penalty, and a gate penalty. The first supplies fitting supervision; the latter two respectively discourage reliance on invalid numerical regions and encourage sparse structure.

\[ \mathcal L_{\mathrm{total}}=\mathcal L_{\mathrm{MSE}}+\lambda_{\mathrm{clamp}}\mathcal L_{\mathrm{clamp}}+\lambda_{\mathrm{gate}}\mathcal L_{\mathrm{gate}}. \]

The appendix uses Adam with learning rate 0.1, 1000 epochs, batch size 512, and early-stopping patience 300, with 1 or 2 hidden layers. The direct route first tries 1 layer and uses 2 if needed; decomposed subproblems use 1 layer. The logarithm lower threshold is 0.005, the exponential upper threshold is 4, the pruning tolerance is 0.01, and the rounding tolerance is 0.001.

Consequently, not requiring users to select operators or tune structure per problem does not mean having no hyperparameters. The method still depends on training settings and window, pruning, rounding, and clamping thresholds; loss weights also enter the objective. Appendix Table 1 does not specify the values of the two loss weights.

Key Experimental Results

Main Results

Evaluation uses 119 Feynman equations, 14 Strogatz ODE problems, and 57 filtered black-box regression problems with continuous features and input dimension no greater than 10. Each Feynman problem supplies 1 million points, from which 10,000 training points are sampled randomly. Results are reported as averages over 3 independent trials.

The symbolic solution rate (SSR) is the fraction of equations recovered exactly; the accuracy solution rate is the fraction of problems exceeding \(R^2>0.999\) in test goodness of fit. Complexity counts expression-tree nodes, with each operator, variable, and constant counting as one node. Black-box tasks have no ground-truth expressions, so only median goodness of fit and complexity can be evaluated; these results do not establish discovery of true equations.

Noise is zero-mean Gaussian perturbation whose standard deviation equals the target root mean square multiplied by the noise coefficient. The tested levels are \(\sigma\in\{0,0.001,0.01,0.1\}\). The table retains conclusions supported by the main text and appendix rather than guessing absolute SSR values from figures whose numerical coordinates are absent from the cache.

Dataset / Condition	Metric	Reported SMILE result	Comparison boundary / Evidence
Feynman, noise-free	SSR	Third place	Behind PySR and ParFam; §4, Figure 3
Feynman, highest noise	SSR	First place; decreases by at most 8% as noise rises	PySR and ParFam decrease by more than 30%; §4
Feynman	Accuracy solution rate	Lower than the best baselines	Shallow networks favor recoverability over maximum fitting accuracy; §4
Strogatz, noise-free	SSR	Second place	Behind PySR; §4, Figure 5
Strogatz, highest noise	Accuracy solution rate	Best; decreases by at most 20% from the noise-free setting	PySR and ParFam decrease by approximately 70%; Appendix E.2
57 black-box problems	Median goodness of fit / Complexity	On the Pareto front; complexity second only to DSR	Median goodness of fit remains below several baselines; Appendix E.3

The 8%, 30%, 20%, and 70% entries retain the source's wording for decreases. The paper does not specify whether these are relative changes or percentage-point changes, so the units should not be reinterpreted. Being on the Pareto front also does not mean highest accuracy: it means no compared alternative improves both accuracy and complexity simultaneously.

Ablation Study

The following table comes from Appendix E.4, which removes one component while keeping the rest unchanged. Ablations introduce a per-stage time budget and count budget violations as timeouts, so reduced recovery includes cases where the pipeline does not finish within the budget.

Config	Reported change from the full pipeline	Time / Failures	Mechanistic interpretation
Without variable analysis and decomposition	SSR decreases by approximately 15%; accuracy solution rate by approximately 30%	Training time roughly doubles; no timeouts	Fitting the entire target increases expression length and structural burden
Without gradient-based constant rounding	SSR decreases by approximately 15%; accuracy solution rate by approximately 30%	Training time roughly doubles; no timeouts	Floating-point coefficients and exponents cannot collapse into exact constants
Without gated pruning	SSR decreases by approximately 35%	A significant fraction of equations time out	Redundant active coefficients make later optimization and rounding costly
Without constant refitting	SSR decreases by approximately 35%	A significant fraction of equations time out	Coefficients are not recalibrated after pruning, impeding simplification by rounding

Ablation percentages likewise retain the original units rather than being relabeled as percentage points. The text provides neither a verifiable exact timeout rate nor the ablation-stage budget in Table 1, so neither is assigned a number here.

Key Findings

• Pruning and constant refitting do more than shorten formulas: they determine whether later stages finish within the time budget, making them controls on the computational bottleneck of recoverability.
• In Appendix Table 3, II.15.4 and II.6.15b succeed at every tested noise level, whereas I.6.2a and I.40.1 succeed only without noise. An average robustness advantage does not imply stability for every equation.
• Greater depth can improve prediction fit, but the authors report that it does not improve SSR: additional degrees of freedom increase expression complexity and make pruning harder. Fitting observations and discovering formulas must be distinguished.
• Timing scope needs clarification: the main text describes SMILE as taking minutes and competitors hours, while Appendix Table 2 reports 43.9 s per problem on average. Hardware is a Xeon Gold 5420+ and RTX A5000 with 24 GB. Without a stated scope, these descriptions cannot be merged into a single precise speedup factor.

Highlights & Insights

• Moving structural cues ahead of fitting is more targeted than merely regularizing a symbolic network. The method uses a variable's disruption of regularity to choose the fitting task, reducing the compositional burden on a shallow model.
• Output sensitivity controls constant rounding rather than a fixed number of decimal places. This can inform post-processing of interpretable models, but actual post-rounding errors and valid domains still need independent checks.
• Exact symbolic recovery allows mathematical equivalence without syntactic identity. Replacing cosine with a sine-based expression reminds evaluators that different trees can represent the same relationship.

Limitations & Future Work

• Power-law cues and single-variable additive or multiplicative decomposition suit physical laws better than arbitrary black-box relationships. More general separability tests and explicit reporting of decomposition-failure fallbacks would broaden the evidence.
• Shallow structure makes pruning tractable but limits fitting of complex equations. Adaptive depth must control recovery cost rather than selecting models solely through higher goodness of fit.
• Real-valued domains of logarithms and powers, behavior in clamped regions, and post-rounding evaluation across domains need separate checks. Low in-data error does not ensure extrapolation free of singularities or undefined inputs.
• Gradient-based rounding uses a local approximation. Checking sensitivity across the parameter interval or directly verifying pre- and post-rounding output changes would strengthen it; generation of simple-constant candidates also needs further specification.
• The source contains inconsistent counts: Appendix C states 252 total problems, but its groups of 119, 14, and 122 sum to 255. The main text describes 15 original baselines plus 4 additions, while Appendix D describes 14 original baselines plus 4 additions yet claims 19 overall. This note retains the explicitly stated evaluated groups without correcting the full collection or baseline count on the authors' behalf.
• Averages over three trials alone do not establish statistical significance. The cached text does not explain a consistent confidence-interval protocol, baseline rerun budgets, or timing scope; speed and recovery claims need checking against complete experimental artifacts.

Related Work & Insights

• vs PySR: PySR evolves discrete expression trees and optimizes constants, whereas SMILE fits a symbolic network continuously before recovering structure. PySR still has higher noise-free Feynman SSR; SMILE's advantages primarily concern robustness and the simplicity / time trade-off.
• vs AI Feynman: Both exploit separable structure, but SMILE checks residual relationships using recovered symbolic subexpressions rather than only neural approximations. Its structural analysis remains heuristic and should not be treated as an entirely new separability principle.
• vs ParFam / EQL: ParFam requires specified function families and polynomial degrees, while EQL shares the starting point of symbolic activations. SMILE connects decomposition, gated pruning, and constant recovery into a pipeline, but does not directly compare against EQL on the full benchmark, so comprehensive empirical superiority over EQL is not established.
• Research direction: Checking expression compactness, correct recovery, and domain validity separately helps prevent attractive formulas from being mistaken for true laws. Stress tests with signed inputs, near-singular regions, and strong variable coupling would be useful extensions.

Rating

• Novelty: 4/5 — The integration of structural analysis, continuous symbolic networks, and discrete recovery is distinctive, with foundations in related approaches.
• Experimental Thoroughness: 3/5 — Ground-truth equations, black-box tasks, and ablations are covered, but counts, timing scope, and statistical reporting need clarification.
• Writing Quality: 3/5 — The pipeline is understandable, but boundaries between universal approximation and shallow practical networks, and between theoretical bounds and endpoint approximations, need sharper treatment.
• Value: 4/5 — The method offers useful ideas for compact formula recovery from noisy scientific data, but should not replace maximum-accuracy regression or validation of scientific laws.