PaperlayerView exact revision on arXiv

Source-linked AI summary

Integration of Neural Network-Based Symbolic Regression in Deep Learning for Scientific Discovery

Samuel Kim, Peter Y. Lu, Srijon Mukherjee, Michael Gilbert, Li Jing, Vladimir Čeperić, Marin Soljačić

arXiv:1912.04825v2cs.LGcs.NEphysics.data-anstat.ML

TL;DR

The paper addresses the need for interpretable and generalizable scientific models beyond black-box neural networks. It integrates the Equation Learner network with deep learning architectures and trains the combined systems end-to-end through backpropagation. Across symbolic-regression, MNIST arithmetic, and dynamical-system tasks, the EQL-based systems recover symbolic structure and extrapolate beyond training data better than the compared ReLU architecture in reported experiments.

  • Problem

    Scientific applications need interpretable and generalizable models, while conventional neural networks are difficult to interpret and often extrapolate poorly.

  • Method

    The paper integrates the Equation Learner symbolic-regression network with other deep learning architectures so the entire system can be trained end-to-end through backpropagation.

  • Results

    Across symbolic-regression, MNIST arithmetic, and dynamical-system tasks, EQL-based architectures learn symbolic relationships and extrapolate beyond training data better than the compared ReLU architecture.

  • Takeaways & Limitations

    The results support combining symbolic regression with deep learning for interpretable and generalizable scientific modeling.

Abstract

from arXiv · show

Symbolic regression is a powerful technique that can discover analytical equations that describe data, which can lead to explainable models and generalizability outside of the training data set. In contrast, neural networks have achieved amazing levels of accuracy on image recognition and natural language processing tasks, but are often seen as black-box models that are difficult to interpret and typically extrapolate poorly. Here we use a neural network-based architecture for symbolic regression called the Equation Learner (EQL) network and integrate it with other deep learning architectures such that the whole system can be trained end-to-end through backpropagation. To demonstrate the power of such systems, we study their performance on several substantially different tasks. First, we show that the neural network can perform symbolic regression and learn the form of several functions. Next, we present an MNIST arithmetic task where a separate part of the neural network extracts the digits. Finally, we demonstrate prediction of dynamical systems where an unknown parameter is extracted through an encoder. We find that the EQL-based architecture can extrapolate quite well outside of the training data set compared to a standard neural network-based architecture, paving the way for deep learning to be applied in scientific exploration and discovery.

I. INTRODUCTION

The paper motivates interpretable, generalizable scientific models and introduces an EQL-based neural architecture that integrates symbolic regression with deep learning and end-to-end backpropagation.

  • Scientific and engineering phenomena can often be described by relatively simple mathematical equations that support broad predictions.
  • Symbolic regression searches mathematical-expression space without predefined features, producing interpretable models that may extrapolate better than black-box neural networks.Genetic-programming approaches can recover physical laws but scale poorly to large systems and may overfit.
  • The proposed EQL architecture integrates symbolic regression with other deep learning architectures while retaining interpretable and generalizable outputs.
  • Because EQL symbolic regression is trainable through backpropagation, the integrated system can be optimized end-to-end without multiple training steps.
  • EQL layers use scientific primitive functions, including component-wise sine or square operations and multi-input multiplication, instead of conventional activations.Stacking layers enables compositions of primitive functions, while depth limits expression complexity.
  • Duplicating activation functions reduces sensitivity to initialization, smooths optimization, and broadens the functions the network can fit.

A. Sparsity

The paper uses sparsity regularization to simplify learned symbolic expressions, emphasizing smoothed L0.5 regularization as a gradient-compatible alternative that strongly drives weights toward zero.

  • Sparsity regularization aims to set inactive weight parameters to zero so they can be removed from the final symbolic expression.
  • L0 regularization promotes sparsity by penalizing nonzero weights, but its NP-hard combinatorial form is incompatible with standard gradient descent.
  • L0.5 regularization remains compatible with gradient descent and is reported to enforce sparsity more strongly than L1 without penalizing weight magnitude as much.
  • Experiments with L0.3 and L0.7 show no significant overall improvement over L0.5, while L0.5 produces simpler expressions than L1.
  • The smoothed L∗ 0.5 regularizer replaces the gradient singularity near zero with a piecewise smoothing function to improve gradient-based training convergence.Experiments use threshold a = 0.01, whereas the accompanying plot uses a = 0.1 for visualization.

III. EXPERIMENTS

Experiments test EQL on analytical functions and MNIST arithmetic, showing recovery of symbolic structure under noise and end-to-end learning that combines digit representation with addition.

  • A. Symbolic Regression: The symbolic-regression evaluation tests analytical functions over x_i ∈ [−1, 1] across 20 random-initialization trials, counting recovery of equivalent correct expressions.
  • A. Symbolic Regression: The network recovers correct expressions with 10% data noise and, with three hidden layers, represents more complicated expressions such as (x_1 + x_2x_3)^3.
  • B. MNIST Arithmetic: The MNIST task learns y = ψ_1 + ψ_2 directly from two 28 × 28 handwritten-digit images, combining image recognition and symbolic regression.
  • B. MNIST Arithmetic: Shared-weight encoders map each digit image to a scalar latent variable, which the EQL network transforms into the predicted sum.
  • B. MNIST Arithmetic: The encoder and EQL network are trained jointly end-to-end with mean-squared error rather than being trained separately.

C. Dynamical System Analysis

The dynamical-system architecture encodes an unknown system parameter and uses recurrent EQL cells to learn propagation equations for kinematics and simple harmonic oscillators.

  • C. Dynamical System Analysis: Each recurrent EQL cell contains a separate EQL network for every predicted state dimension; for y = {u, v}, it uses two EQL networks.
  • C. Dynamical System Analysis: The dynamics encoder produces a single latent parameter from an input time series, while the propagating decoder predicts future states from that parameter and an initial condition.
  • C. Dynamical System Analysis: The full dynamics architecture is trained end-to-end with mean-squared error, without separately training the dynamics encoder and propagating decoder.
  • C. Dynamical System Analysis: The experiments use datasets from kinematics and simple harmonic oscillator systems, with position u and velocity v as state variables.
  • C. Dynamical System Analysis: In kinematics, the relevant varying parameter is acceleration a = F/m, and the encoder output is expected to correlate with it.

2) Simple Harmonic Oscillator (SHO):

The simple harmonic oscillator is modeled as a coupled first-order ODE, numerically integrated with finite differences, while an encoder extracts the resonant-frequency parameter for dynamical prediction.

  • The SHO describes position and velocity dynamics governed by resonant frequency ω, with ω = √(k/m) for a spring.
  • The system uses an Euler finite-difference method to integrate the SHO equations numerically.
  • State variables u and v are sampled at Δt = 0.1, with propagation length increased during training to avoid learning the identity function or unstable trajectories.
  • The initial conditions u0 and v0 feed the propagator, while the dynamics encoder is trained to correlate with ω².

D. Training

The MNIST arithmetic system is trained end-to-end with regularization that sparsifies the EQL component, learning digit representations and an additive equation while testing extrapolation beyond the training range.

  • D. Training: The network is trained by backpropagation with RMSProp using a loss that balances mean-squared error against regularization through λ.
  • D. Training: Training optionally begins with small λ, then increases λ to sparsify EQL weights before thresholded weights are zeroed, frozen, and fine-tuned.
  • A. MNIST Arithmetic: The encoder learns a strong linear mapping from MNIST images to the latent digit variable despite never receiving digit labels.
  • A. MNIST Arithmetic: The extracted EQL equation captures the additive structure, with mean absolute errors of 0.307 on training digits and 0.315 on test digits.
  • A. MNIST Arithmetic: Rounded regression outputs yield accuracies of 89.7% on training digit pairs and 90.2% on test digit pairs.
  • A. MNIST Arithmetic: For pairs with y ≥ 15, the EQL architecture loses only a few percentage points, whereas the ReLU replacement falls below 1% accuracy.

B. Kinematics

The EQL-based architectures extract latent dynamical parameters and discover propagation equations, while extrapolating beyond training regimes better than conventional ReLU networks. In the SHO task, EQL also captures higher-order corrections beyond Euler integration.

  • Kinematics: Both EQL and ReLU networks match the true kinematics solution in the training regime, but ReLU quickly diverges outside it while EQL extrapolates for several more time steps.The ReLU comparator uses two fully connected hidden layers with 50 neurons each and ReLU activations.
  • Simple Harmonic Oscillator: The SHO dynamics encoder yields a strong linear z–ω2 relation, with correlation coefficient −0.995, because ω2 requires fewer propagation operations than ω.Using ω would require an additional squaring function in the propagating equations.
  • Simple Harmonic Oscillator: A residual term 0.0133z2 in the extracted SHO expression may arise from cross-terms or incomplete convergence and could disappear after further thresholding.The authors identify this term as likely removable rather than part of the expected equation.
  • Simple Harmonic Oscillator: Within the SHO training regime, all methods fit the true solution reasonably well, but ReLU fails to extrapolate and regresses to noise around 0.Euler and EQL extrapolate reasonably well, while Runge–Kutta almost exactly fits the analytical solution.
  • Simple Harmonic Oscillator: The EQL SHO solution has smaller error than Euler, indicating that it learned higher-order corrections to the first-order Euler method.Its extracted equations match the Euler approximation in their leading terms and include an additional term corresponding to a second-order Taylor contribution.

V. CONCLUSION

The paper integrates symbolic regression with deep learning architectures in end-to-end systems spanning MNIST arithmetic and dynamical-system prediction. The authors report that SHO results suggest possible discovery of improved ODE integration techniques, while noting optimization and architecture-extension challenges.

  • Conclusion: The integrated architecture trains symbolic regression and other deep learning components end-to-end through backpropagation.The demonstrated systems include MNIST digit arithmetic and dynamical-system prediction with an extracted unknown parameter.
  • Conclusion: In the SHO system, the results suggest that symbolic regression could discover new ODE integration techniques, including improved integrators for stiff systems.The paper frames this as a potential direction rather than an established result.
  • Limitations and Future Work: The EQL network is not always able to find the correct mathematical expression because gradient-based optimization can become trapped in local minima.The authors mitigate this through increased stochasticity, duplicated activation functions, and multiple trials, but seek alternatives to repeated trials.
  • Limitations and Future Work: Future work includes integrating EQL with spatio-temporal architectures for PDE discovery and extending encoders to videos, audio, and text.Other proposed extensions include time-dependent unknown parameters and known or learnable spatial-derivative operators.

APPENDIX A EQL NETWORK DETAILS

The EQL network uses configurable activation functions, staged training, and sparsity methods to discover compact symbolic expressions. Its computational cost is comparable to conventional fully connected networks, while symbolic regression remains combinatorially hard.

  • EQL architecture: EQL hidden layers combine primitive activation functions whose duplication improves optimization robustness and broadens the representable function range.The exact duplication count is arbitrary and has little effect on performance.
  • Benchmarking: 20 trials benchmark whether the network reaches each correct trial function at least once, selecting solutions using equation simplicity and test error.Equivalent expressions and slight coefficient variations count as correct.
  • Limitations: The method does not provide asymptotic speedup because finding the correct expression remains an NP-hard combinatorial search.Its stated advantage is integrating symbolic regression with deep learning through gradient descent.
  • Computational efficiency: 37 seconds for 2 hidden layers and 51 seconds for 3 hidden layers on an Nvidia GTX 1080 Ti.These timings correspond to 20,000 training epochs.
  • Computational efficiency: The EQL network and conventional fully connected networks share O(n^2) dominant weight-matrix complexity, despite EQL activation evaluation taking O(n).The activation-function iteration is therefore not the dominant computational cost.
  • Relaxed L0 regularization: Relaxed L0 regularization uses stochastic gates and an analytical expected penalty to encourage sparse expressions during training.Unlike the basic procedure, this approach avoids a final stage that freezes small weights at zero.

APPENDIX C EXPERIMENT DETAILS

The MNIST arithmetic experiment uses a convolutional encoder and normalized scalar output to extract digit information and predict the arithmetic label, with a ReLU network as comparison.

  • Encoder: The encoder uses two convolutional-max-pooling stages followed by fully connected ReLU layers and a one-unit batch-normalized output.The convolutional layers use 32 and 64 filters, each of size 5 × 5.
  • Output normalization: The predicted output is scaled as ˆy = 9ˆy∗ + 9 before loss computation to match the expected label range.This is equivalent to normalizing y to [−1, 1].
  • Comparison model: The comparison ReLU network has two hidden layers with 50 units each and uses two-stage training with regularization only in the first phase.Small weights are frozen between phases at threshold α = 0.01.

B. Kinematics

The dynamics experiments generate time-series data for sampled parameters and use convolutional encoders with EQL or ReLU predictors under staged training schedules.

  • Kinematics: The kinematics dataset samples 100 values of a and propagates each input series for Tx = 100 time steps.The passages specify the sampling count and propagation length but not the explicit kinematic equation.
  • Kinematics: The kinematics dynamics encoder uses two 1D convolutional layers, a 16-node ReLU hidden layer, one output unit, and batch normalization.Each convolutional layer has 16 filters of length 5.
  • Kinematics: Kinematics training uses an initial phase of 5000 iterations followed by a second phase of 10000 iterations without regularization.The first phase uses learning rate 10^-2 and λ = 10^-3; small weights are frozen at α = 0.1.
  • Simple harmonic oscillator: The SHO dataset samples 1000 values of ω2, propagates inputs for Tx = 500 steps, and predicts Ty = 25 future steps at ∆t = 0.1.Input and output series use the same time step.
  • Simple harmonic oscillator: The SHO dynamics encoder retains the kinematics encoder architecture, while the EQL network reduces activation duplication because it propagates more time steps.The ReLU comparison uses four hidden layers with 50 units each.
  • Simple harmonic oscillator: SHO training progressively increases the prediction horizon from Ty = 1 to Ty = 7 before using Ty = 25 in the second phase.The supplied schedule also changes learning rate and regularization across phases.

APPENDIX D ADDITIONAL MNIST ARITHMETIC DATA

Additional MNIST trials show that the encoder can produce different latent representations while the EQL component still recovers the digit-sum relationship. Exact equations vary with initialization, but one-to-one latent mappings can remain sufficient for prediction.

  • Trial variability: Different random initializations produce variable EQL equations, motivating additional trials to assess robustness.The reported Figure 6 and Table I results come from one of several trials.
  • Trial 2: Trial 2 yields a linear relationship between the true digit φ and latent variable z, with positive rather than negative correlation, while accurately predicting the sum y.The sign of the correlation is arbitrary because no bias is imposed on z.
  • Trial 3: Trial 3 preserves a one-to-one mapping between φ and z even though the relationship is nonlinear, and the EQL network still predicts the correct sum.The result is reported for both training and test datasets in Figure 10.
  • Interpreting the learned equation: The Trial 3 equation first inverts an inverse-sine-like latent transformation and then adds the two digits.Analyzing the latent variable and resulting equation provides insight even when the exact expected equation is not recovered.
Loading 1912.04825v2…