Efficient Constant Optimization for Symbolic Regression with GPU-Accelerated Tree-Based Genetic Programming

Symbolic regression discovers closed-form expressions from data by jointly searching over both the structure and numerical parameters of candidate models. Tree-based genetic programming is the dominant approach, evolving a population of expression trees through fitness-based selection and genetic operators. However, genetic operators are ineffective at fine-tuning real-valued constants, leading to poor fitness for structurally correct expressions with suboptimal coefficients. Constant optimization addresses this by solving a continuous nonlinear least-squares subproblem for each fixed tree, typically using the Levenberg-Marquardt method.

The computational cost of constant optimization often leads modern symbolic regression systems to omit or restrict it. While GPU-accelerated frameworks improve evaluation throughput, they lack support for constant optimization, leaving a gap between structural search and the precision needed to recover governing equations.

This paper designs a GPU-resident batched Levenberg-Marquardt solver purpose-built for the heterogeneous constant optimization workload in TGP-based symbolic regression, integrating it directly into EvoGP so that constant optimization runs entirely on the GPU within each generation. The use of reverse-mode automatic differentiation and a double-precision delivery guard ensures the precision and efficiency of optimization results.

• �� Design a GPU-resident batched LM solver for heterogeneous population constant optimization.
• �� Use reverse-mode automatic differentiation to construct the Jacobian matrix for each tree, independent of constant count.
• �� Implement each iteration with a fixed number of CUDA launches, using a double-precision delivery guard to ensure result precision.

Experiments were conducted on an NVIDIA A100 using synthetic populations calibrated to snapshots from real EvoGP runs. Benchmark tests used synthetic targets with known near-global optima, allowing direct assessment of solution quality and fair throughput comparisons. Populations ranged from 10^3 to 2.56×10^5 trees, with 100 to 10,000 data points per tree. All methods received the same trees, data, and initial constants.

The solver sustains up to 5.1×10^5 trees per second on an NVIDIA A100 for early-generation populations. In a GPU-saturated benchmark configuration, it delivers roughly 9.9× the throughput of Operon running on a 64-core EPYC 7763. Integrated into EvoGP, the solver recovers governing equations on 10 of 18 constructed problems, compared to 0 for stock EvoGP.

The method can be applied to symbolic regression tasks in scientific discovery and interpretable machine learning, particularly where complex models need optimization. Its efficient constant optimization capability makes it suitable for applications requiring rapid iteration and high precision.

While the method performs excellently on GPUs, GPU memory may become a bottleneck when handling very large datasets. Additionally, the algorithm's reliance on GPU hardware may limit its application on other hardware. Future research directions include optimizing the algorithm to support larger datasets and exploring its potential application on other hardware architectures.