Neural Symbolic Regression that Scales

Problem

• Symbolic regression infers both equation structure and parameters from input-output pairs, making it a difficult combinatorial problem.The number of functions grows exponentially with expression length, while numerical constants further increase the difficulty.Source: §1 · p. 1

Limitations of prior methods

• Traditional search-based methods regress each equation from scratch and therefore do not improve when given data from different equations.Their fitness-driven search strategies do not accumulate experience across problems.Source: §1 · p. 1
• Traditional methods also have an opaque inductive bias that is difficult for users to steer toward a specific class of equations.This limits transparent control over the prior used during regression.Source: §1 · p. 1

Approach

• The proposed framework learns symbolic regression from data by procedurally generating unbounded, accurate training data at negligible cost.The pre-training equation distribution determines the final model’s prior and makes that prior easier to understand and control.Source: §1 · p. 1

Contributions

• The paper introduces the first reported symbolic regression approach that improves over time with data and compute.It is presented as a simple, flexible, and powerful framework.Source: §1 · p. 2
• The learned method significantly outperforms state-of-the-art approaches that rely on hand-designed strategies.The paper also releases its code and largest pre-trained model.Source: §1 · p. 2

Genetic programming

• Genetic programming evolves populations of candidate mathematical expressions through mutation and recombination.Eureqa is identified as a prominent commercial genetic-programming system for symbolic regression.Source: §2 · p. 2

Neural approaches

• Neural symbolic-regression networks replace standard activations with symbolic building blocks so formulas can be read from network architectures and weights.These models handle high-dimensional data and scale with the number of input-output pairs.Source: §2 · p. 2
• Activation-based neural methods require extensions for divisions and can encounter exploding gradients and numerical issues with exponential or logarithmic activations.The limitation is tied to the operations represented by the network.Source: §2 · p. 2

Latent-space search

• Variational autoencoder methods move symbolic-expression search from a discrete space into a continuous latent space.The autoencoder reconstructs expressions before a subsequent latent-space search for a fitting function.Source: §2 · p. 2

Hybrid approaches

• Other neural methods identify properties such as additive separability and compositionality, then recursively simplify data into subproblems for symbolic-regression techniques.These approaches combine neural networks with existing methods or hand-designed rules.Source: §2 · p. 2

Large-scale pre-training

• Large-scale pre-training has been applied to computer vision, natural language processing, symbolic mathematics, and automated theorem proving.The paper builds on Transformer results for symbolic integration and differential-equation solving, while mapping numerical values to symbolic expressions remains an added challenge.Source: §2 · pp. 2–3

Framework

• A symbolic regressor maps a set of input-output pairs to an equation whose function approximates the observed outputs.The learned model is parameterized as Sθ and trained from many examples.Source: §3 · p. 3

Data generation

• The pre-training dataset contains procedurally generated equations and input sets drawn from a chosen distribution that determines the model’s inductive bias.The distribution can vary equation families, constants, term counts, interactions, input-set sizes, and sampling dependencies.Source: §3.1 · p. 3

Pre-training

• The Transformer’s encoder converts numerical input-output data into a fixed-size latent representation, while the decoder predicts the corresponding equation sequence.Training examples are generated in potentially unbounded quantity from standard mathematical libraries.Source: §3.1 · p. 3
• The decoder predicts an equation skeleton in which numerical constants are replaced by placeholder symbols.For example, 4.2 sin(0.3x1) + x2 becomes ⋄sin(⋄x1) + x2, and the model is trained against this skeleton.Source: §3.1 · p. 3

Test-time search

• At test time, the encoder represents observed pairs, the decoder samples candidate skeletons, and a nonlinear optimizer fits each placeholder constant.Candidate equations are evaluated by minimizing loss between their predictions and the observed targets.Source: §3.2 · p. 3

• NeSymReS is the evaluated implementation of the proposed neural symbolic regression framework.Source: §4 · p. 3

4.2. Pre-training Data Generator

• The data generator samples equation skeletons, numerical constants, and support points, then evaluates equations to create input-output pairs and skeleton targets.Training targets omit numerical constants, and both inputs and targets use prefix notation.Source: Figure 1 · p. 4

4.1. The Model Sθ

• The Set Transformer encoder uses inducing points so computation scales as O(nm) rather than O(n^2), with m = 50, while the decoder is a regular Transformer.The encoder and decoder have 11 million and 13 million parameters, respectively.Source: §4.1 · p. 4

4.3. Inference

• At test time, beam-search generates candidate skeletons and BFGS fits their constants by minimizing squared loss against the observed outputs.The default configuration uses beam-size 32 and four BFGS restarts per equation.Source: §4.3 · p. 5

4.4. Evaluation

• Evaluation uses five datasets, including 52 AI-Feynman equations, 200 strictly out-of-sample equations, and the 12-equation Nguyen benchmark.The SOOSE variants differ in their treatment of constants, with full-constant equations outside the model's pre-training distribution.Source: §4.4 · p. 5

(i) Accuracy as a Function of Pre-training Data

• NeSymReS accuracy steadily increases as the pre-training dataset grows from 10K to 10M equation skeletons.All models use the same training procedure and iterations, with early stopping on held-out validation data.Source: (i) Accuracy as a Function of… · p. 6

(ii) Accuracy as a Function of Test-time Compute

• Across five datasets, NeSymReS outperforms all baselines in both time and accuracy across most compute budgets.The comparisons use matched or greater computing time for each baseline and run methods on a single CPU.Source: (ii) Accuracy as a Function of… · p. 6, (i) Accuracy as a Function of… · p. 6
• On AIF, NeSymReS reaches the same maximum accuracy as Genetic Programming more than three orders of magnitude faster on CPU.Source: (ii) Accuracy as a Function of… · p. 6
• NeSymReS performs relatively well on Nguyen at shorter runtimes and comparably with DSR at approximately 10^3 seconds.DSR was fine-tuned on two equations from the Nguyen benchmark.Source: (ii) Accuracy as a Function of… · p. 7

(iii) Performance Improves with more Points

• NeSymReS remains robust across different numbers of test points, and its accuracy increases as more input-output pairs become available.The evaluation varies the number of observed points from 1 to 1024; the model was pre-trained with no more than 500 points.Source: (iii) Performance Improves with more… · p. 7, Figure 5 · p. 7

Discussion

• Pre-training on millions of equations lets NeSymReS outperform several strong baselines and improve as the dataset grows.The authors emphasize that procedurally generated training data can be produced efficiently and indefinitely.Source: §6 · p. 8
• The pre-training equation distribution controls the model’s inductive bias, including its support, constants, terms, and interactions.This makes the prior targetable to specific applications.Source: §6 · p. 8
• The encoder uses a Set Transformer whose computation scales linearly with the number of input-output points and input dimensions.The reported comparison is O(nm) for the Set Transformer versus O(n^2) for a standard Transformer encoder.Source: §6 · p. 8, §A.1 · p. 11
• BFGS constant fitting can struggle with several local minima, and the presented model cannot test equations with more variables than seen during pre-training.The neural network also does not directly interact with the function evaluator.Source: §6 · p. 8
• The framework aims to scale with datapoints and variables, improve with experience, and target specific equation distributions.These properties are presented as desirable characteristics of a strong symbolic regressor.Source: §6 · p. 8

Baselines

• Deep Symbolic Regression uses standard open-source hyper-parameters, including final estimation of numerical constants.Its main hyper-parameters are the entropy coefficient λH and risk factor ϵ.Source: §A.2 · p. 12
• Genetic Programming is implemented with gplearn using mostly documented defaults and a specified operator set.The available symbols include arithmetic, square root, logarithm, exponential, negation, inversion, sine, and cosine.Source: §A.2 · p. 12
• Gaussian Process regression is included as a non-symbolic baseline that directly learns a mapping from x to y.The implementation uses sklearn defaults with a product of constant and RBF kernels.Source: §A.2 · p. 12
• Missing primitive functions lower a method’s maximum attainable performance on datasets containing those primitives.The same restriction can make other equations easier to fit because the search function set is smaller.Source: §A.2 · pp. 12–13

Training Dataset Generation

• Expression trees contain at most five non-leaf nodes, with operator choices sampled from the experiment’s weighted distribution.Leaf nodes are sampled as independent variables with probability 0.8 and integers with probability 0.2.Source: §B.1 · p. 13, Table 6 · p. 13
• The training corpus contains 10M equations compiled into evaluable functions so support points and constants can be resampled during mini-batch training.The functions are compiled with SymPy’s lambdify.Source: §B.1 · p. 13

Training

• Mini-batches contain B = 150 equations, with up to min(3, Nc) constants sampled uniformly from U(1, 5).The number and locations of constants are randomized during training.Source: §B.1 · p. 13

Training Dataset Distribution

• The 10M-equation dataset contains approximately 1.2M unique skeletons, with shorter expressions occurring more frequently.Some skeletons are repeated or mathematically equivalent.Source: §B.1 · p. 13

Addition of Numerical Constants

• Constant placeholders are added to generated or predicted skeletons by multiplying unary operators and adding constants to independent variables.Longer expressions consequently tend to receive more placeholders.Source: §B.1 · pp. 13–14

Evaluation Metrics

• Reported results are accuracies over all equations, with plot error bars representing the standard error of the mean.This convention applies across methods and datasets.Source: §B.2 · p. 14
• The evaluation distinguishes in-support accuracy Aiid from out-of-support accuracy Aood by extending each variable’s support beyond the in-sample interval.The extension length equals hi − lo for each variable present in the equation.Source: §B.2 · p. 14

SOOBE Dataset

• SOOBE uses equation skeletons that do not overlap numerically or symbolically with the pre-training dataset.Its inputs are sampled from U(−10, 10), using 500 support points per independent variable.Source: §B.2 · p. 14

Benchmarks

• The model is evaluated on AI-Feynman, SOOBE-WC, SOOBE-NC, SOOBE-FC, and Nguyen.The appendix lists the AI-Feynman equations and samples 50 of 200 SOOBE equations.Source: §B.2 · p. 14, Table 8 · p. 15, Table 9 · p. 16

• A1-metric conclusions are consistent with those drawn using the A2 metric.Source: §C.1 · p. 17
• Accuracy is plotted against pre-training dataset size under a fixed test-time computational budget of approximately 100 seconds.Source: Figure 7 · p. 17
• In-distribution and out-of-distribution accuracy are each evaluated as functions of time on a single CPU per equation.Source: Figure 8 · p. 17, Figure 9 · p. 17
• Accuracy is also evaluated as a function of the number of input-output pairs observed at test time.Source: Figure 10 · p. 17
• NeSymReS extrapolates out of distribution from relatively few support points by retrieving underlying symbolic expressions, sometimes finding valid trigonometric alternatives.The examples involve functions of two independent variables, x1 and x2; alternative forms use identities such as sin(x) = cos(x −π/2) = −cos(x + π/2).Source: §C.2 · p. 18, Figure 11 · p. 18