Ex-HiDeNN: Hierarchical Neural Symbolic Regression

• Ex-HiDeNN is a two-stage framework that first learns a smooth, separable neural surrogate using C-HiDeNN-TD and then employs symbolic regression to extract closed-form expressions.
• It incorporates Hessian-based analysis to quantify additive and multiplicative separability, guiding the regression process to reduce complexity and enhance accuracy.
• The framework demonstrates superior performance in noise robustness and RMSE reduction compared to traditional symbolic regression methods, especially for multiplicative and moderately separable functions.

Ex-HiDeNN, short for Explainable Hierarchical Deep Learning Neural Networks, is a two-stage framework for discovering closed-form, human-interpretable mathematical expressions from data by combining a structured interpolating neural surrogate with symbolic regression. Its core components are C-HiDeNN-TD, a convolutional Hierarchical Deep Neural Network with Tensor Decomposition, and PySR, a genetic-programming-based symbolic regression engine. The method is designed to be accurate, frugal, fast, separable, and scalable, with an explicit emphasis on recovering symbolic surrogates from limited observation while exploiting additive and multiplicative separability in the learned function (Batley et al., 7 Jul 2025).

1. Conceptual definition and lineage

Ex-HiDeNN is formulated as a bridge between hierarchical neural interpolation and symbolic discovery. In the first stage, it learns a continuously differentiable, structured surrogate of the target function from data. In the second stage, it uses derivatives of that surrogate to assess separability and then adapts the symbolic regression strategy accordingly. The stated goals are interpretability, accuracy, frugality or parsimony, speed and scalability, and separability-aware modeling (Batley et al., 7 Jul 2025).

The framework is built upon the authors’ earlier work on HiDeNN, Convolutional HiDeNN, and C-HiDeNN-TD. HiDeNN embeds computational-mechanics-style interpolation into neural networks. Convolutional HiDeNN and C-HiDeNN-TD extend that basis through patch functions, convolutional ideas, and tensor decomposition, yielding scalable and separable interpolating networks. Within Ex-HiDeNN, these prior developments serve as the structured substrate on which symbolic regression operates.

A central feature of Ex-HiDeNN is that symbolic regression is not applied directly to raw data. Instead, the surrogate first interpolates and smooths the data, and only then provides samples and structural information for symbolic regression. This arrangement is presented as the main reason the method can discover closed-form expressions that are more accurate and robust than traditional symbolic regression alone, often with orders-of-magnitude improvements in RMSE on the reported benchmarks.

2. Surrogate architecture and mathematical representation

The surrogate used in Ex-HiDeNN is C-HiDeNN-TD, a separable, hierarchical neural network for approximating a scalar function $u(\mathbf{p})$ with $I$ input parameters,

$\mathbf{p} = [p_1, p_2, \dots, p_I].$

Its learned mapping is written as

$$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$

Here, $i=1,\dots,I$ indexes input dimensions, $m=1,\dots,M$ indexes modes, $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$ is a one-dimensional interpolating neural basis along dimension $i$, and $u_{pi,k}^m$ are trainable coefficients (Batley et al., 7 Jul 2025).

Conceptually, the functional space is split into multiple one-dimensional subspaces, one for each input dimension. Along each axis, C-HiDeNN-TD constructs a meshfree-like interpolant with trainable patch functions. The full function is then represented as a sum over modes of products over dimensions of these one-dimensional interpolants, in a form described as structurally akin to a tensor decomposition, specifically canonical polyadic decomposition.

The one-dimensional interpolant for a parameter $p$ in a segment $I$0 is expressed as

$I$1

In this formulation, $I$2 are finite-element-like shape functions at nodes $I$3 of segment $I$4, and $I$5 are patch enrichment functions over a neighborhood of size $I$6, with dilation parameter $I$7 and polynomial order $I$8 as hyperparameters. Because these interpolants are built from specific basis and patch functions, the surrogate is described as small, structured, and continuously differentiable.

Training is posed as least-squares fitting to observed data through minimization of MSE. The loss is written as

$I$9

This representation is important because it embeds separable structure directly into the neural parametrization rather than attempting to infer it only after fitting.

3. Separability analysis and the two-step algorithm

Ex-HiDeNN relies on classical notions of additive and multiplicative separability. For additive separability,

$\mathbf{p} = [p_1, p_2, \dots, p_I].$0

and the Hessian becomes diagonal. For multiplicative separability,

$\mathbf{p} = [p_1, p_2, \dots, p_I].$1

Taking $\mathbf{p} = [p_1, p_2, \dots, p_I].$2 converts multiplicative separability into additive separability,

$\mathbf{p} = [p_1, p_2, \dots, p_I].$3

This is the basis for analyzing the Hessian of $\mathbf{p} = [p_1, p_2, \dots, p_I].$4 to detect multiplicative structure (Batley et al., 7 Jul 2025).

Given a scalar function $\mathbf{p} = [p_1, p_2, \dots, p_I].$5, Ex-HiDeNN defines the Hessian of $\mathbf{p} = [p_1, p_2, \dots, p_I].$6 as

$\mathbf{p} = [p_1, p_2, \dots, p_I].$7

with diagonal matrix $\mathbf{p} = [p_1, p_2, \dots, p_I].$8. The multiplicativity score is

$\mathbf{p} = [p_1, p_2, \dots, p_I].$9

The score satisfies $$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$0, with $$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$1 corresponding to perfect multiplicative separability and $$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$2 corresponding to purely coupled behavior. The score is computed pointwise via automatic differentiation of the trained C-HiDeNN-TD surrogate and then averaged over sampled points to obtain a domain-level separability measure.

The full pipeline is organized as a two-step algorithm. First, C-HiDeNN-TD is trained on the data to obtain an interpolant $$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$3. Second, the algorithm computes the relevant Hessian-derived separability score and branches the symbolic regression procedure according to thresholded separability regimes. The formal regimes are summarized below.

Regime	Condition on $$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$4	Symbolic form
High separability	$$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$5	$$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$6
Moderate separability	$$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$7	$$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$8
Low separability	$$u^h(\mathbf{p}_k)=u^h(p_{1,k},p_{2,k},...,p_{I,k})= \sum_{k=1}^{K}\prod_{i=1}^I\sum_{m=1}^{M}\tilde N_{pi,k}(p_{i,k};\mathbf{\theta})\cdot u_{pi,k}^m.$$9	Full multivariate symbolic regression

In the high-separability regime, Ex-HiDeNN performs one-dimensional symbolic regression independently for each variable and multiplies the resulting expressions. In the moderate regime, it uses a multi-mode sum of products of one-dimensional symbolic pieces. In the low-separability regime, it performs global symbolic regression on multivariate samples from the surrogate. This adaptive branching is presented as the mechanism that reduces combinatorial blow-up by constraining symbolic search to forms compatible with the structural information extracted from the surrogate.

4. Symbolic regression engine, implementation, and complexity

The symbolic regression component in Ex-HiDeNN is PySR, a multi-population evolutionary symbolic regression system built atop the Julia library SymbolicRegression.jl. Its evolutionary loop consists of initializing a population of candidate expressions, selecting a random subset for tournament selection using a fitness function such as MSE with optional complexity penalty, accepting a parent with a simulated-annealing-based probability, mutating the parent through operators such as $i=1,\dots,I$0, $i=1,\dots,I$1, $i=1,\dots,I$2, $i=1,\dots,I$3, $i=1,\dots,I$4, $i=1,\dots,I$5, and $i=1,\dots,I$6, and replacing an existing individual, by default the oldest (Batley et al., 7 Jul 2025).

Two PySR design features are highlighted. The first is age-regularized evolution, in which the oldest individual is replaced rather than the least fit, with the stated purpose of helping avoid premature convergence. The second is simulated annealing acceptance,

$i=1,\dots,I$7

where $i=1,\dots,I$8 and $i=1,\dots,I$9 are pre- and post-mutation fitness and $m=1,\dots,M$0 and $m=1,\dots,M$1 are hyperparameters. Within Ex-HiDeNN, aside from population size, maximum iterations, and complexity, most PySR hyperparameters are reported to use defaults.

Ex-HiDeNN also specifies standard performance metrics: Mean Square Error, Root Mean Square Error, Maximum Absolute Error, Relative Error, and coefficient of determination $m=1,\dots,M$2. For sampling, it uses Latin Hypercube Sampling for global coverage and Local Perturbation Sampling for local jitter around training points. The paper states that this combination helps symbolic regression access both global structure and local detail.

Representative training configurations are case-dependent. The demonstrative example uses 42 parameters, 100 epochs, and attains RMSE $m=1,\dots,M$3. In the fatigue example, the surrogate uses 600 parameters and 2000 epochs. In the hardness example, it uses 630 parameters and 4000 epochs. For the Lorenz system, the configuration is 581 parameters per output and 1708 epochs, with normalized RMSE $m=1,\dots,M$4. Optimization is described as standard neural training with gradient-based methods, while the exact optimizer is not specified.

A central implementation claim is that symbolic regression is NP-hard and its search space grows exponentially with dimension. Ex-HiDeNN is presented as mitigating this through dimensionality reduction in symbolic space, the tensor-decomposition architecture of C-HiDeNN-TD, and parallelism across independent one-dimensional regression problems. This suggests that its computational advantages are strongest when the target function exhibits substantial separable structure.

5. Benchmarks and dynamical system identification

The paper first presents a demonstrative example on

$m=1,\dots,M$5

Using 10,000 Latin Hypercube samples with a 70/15/15 split, C-HiDeNN-TD with 42 parameters trained for 100 epochs achieves RMSE $m=1,\dots,M$6. PySR with 1000 iterations, maximum complexity 15, and 20 populations recovers one-dimensional expressions in $m=1,\dots,M$7 and $m=1,\dots,M$8 that combine into $m=1,\dots,M$9 with RMSE $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$0. The separability score is described as essentially 1, and the example is used to illustrate near-exact multiplicative separability (Batley et al., 7 Jul 2025).

A broader evaluation is conducted on the Vladislavleva benchmark suite $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$1 through $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$2, with comparisons against original genetic programming, KAN, and PySR alone. Median RMSE over five runs is reported. Ex-HiDeNN recovers the exact functional form up to constants on multiple problems, including $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$3, $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$4, $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$5, and $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$6 to rounding. On $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$7, Ex-HiDeNN attains RMSE $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$8 compared with $\tilde N_{pi,k}(p_{i,k};\boldsymbol{\theta})$9 for PySR, $i$0 for KAN, and $i$1 for the reference method. On $i$2, PySR alone is slightly better in RMSE, though the paper notes it is more expensive. The overall pattern reported is that Ex-HiDeNN performs especially well on multiplicative and moderately to highly separable functions.

Noise robustness is examined by adding Gaussian noise $i$3 to $i$4 for $i$5 between $i$6 and $i$7. Across these noise levels, Ex-HiDeNN achieves a $i$8–$i$9 reduction in RMSE relative to PySR alone. The interpretation given is that the interpolating surrogate filters noise and captures the underlying signal before symbolic regression is performed, whereas direct symbolic regression over noisy data overfits more easily.

For dynamical system identification, the target is the Lorenz system,

$u_{pi,k}^m$0

with $u_{pi,k}^m$1, $u_{pi,k}^m$2, and $u_{pi,k}^m$3. The trajectories are integrated numerically from $u_{pi,k}^m$4 to $u_{pi,k}^m$5 with timestep $u_{pi,k}^m$6, starting from $u_{pi,k}^m$7. A C-HiDeNN-TD surrogate with one input, three outputs, 581 parameters per output, and 1708 epochs is trained, then differentiated by automatic differentiation to obtain $u_{pi,k}^m$8, $u_{pi,k}^m$9, and $p$0. PySR is applied separately to each derivative component using the operator dictionary $p$1. The recovered equations are

$p$2

These are reported to match the true Lorenz parameters to two significant digits.

6. Engineering applications, comparative position, and limitations

Three engineering applications are presented. The first concerns discovery of a closed-form fatigue equation from the NIMS rotation-bending fatigue dataset from Agrawal et al. (2014). The target is fatigue life in units of $p$3 cycles, with 25 inputs spanning alloy composition, heat treatment parameters, and processing variables. The dataset contains 437 experiments, corresponding to about $p$4 samples per dimension. The surrogate uses 600 parameters, specifically 6 modes and 4 nodes per each of the 25 inputs, trained for 2000 epochs. Its training performance is reported as RMSE $p$5 and RE $p$6. Sampling the surrogate at 520 points and applying PySR yields

$p$7

For this symbolic expression, training RMSE is 34, RE is $p$8, and $p$9; test RMSE is 20, RE is $I$00, and $I$01 (Batley et al., 7 Jul 2025).

The second application concerns prediction of Vickers hardness from micro- and nano-indentation-related descriptors using data from Tantardini et al. (2024). The six inputs are $I$02, $I$03, $I$04, $I$05, $I$06, and $I$07. The C-HiDeNN-TD surrogate uses 630 parameters, specifically 15 modes and 7 nodes per feature, trained for 4000 epochs, with normalized training RMSE $I$08. Sampling at 621 points and running PySR gives

$I$09

Its performance is reported as training RMSE $I$10 GPa, MaxAE $I$11 GPa, RE $I$12, and $I$13; test RMSE $I$14 GPa, RE $I$15, and $I$16; and full-dataset RMSE $I$17 GPa with $I$18. The paper compares this with a SISSO expression whose full-dataset performance is RMSE $I$19 GPa and $I$20, describing the Ex-HiDeNN result as about a 25-fold reduction in RMSE.

The third application targets the yield surface of a Matsuoka–Nakai-type pressure-sensitive granular material using data from Bahmani et al. (2024). Data are generated over pressure $I$21, Lode angle $I$22, and Lode radius $I$23, then converted into principal stresses and invariants $I$24, $I$25, and $I$26, producing 13,200 nine-dimensional samples with a 70/15/15 split. The discovered expression is

$I$27

The reported performance is training RMSE $I$28 MPa, RE $I$29, and $I$30, and test RMSE $I$31 MPa, RE $I$32, and $I$33. The accompanying interpretation states that the prefactor $I$34 is consistent with pressure-sensitive yield criteria and that the expression is reminiscent of the classical Matsuoka–Nakai invariant condition.

Relative to pure symbolic regression methods such as PySR, SISSO, and earlier genetic programming, Ex-HiDeNN is described as benefiting from surrogate-based denoising, derivative information, and separability-aware dimensionality reduction. Relative to KAN and related grey-box models such as INNs and KHRONOS, it is described as retaining neural expressivity while yielding fully symbolic, parsimonious models. Relative to SINDy for dynamical systems, it differs by learning a neural surrogate in time, differentiating it, and then applying symbolic regression rather than fitting a pre-defined polynomial library.

The explicit limitations identified are dependence on surrogate quality, the NP-hardness and operator-dictionary dependence of symbolic regression, reduced advantage on strongly coupled high-dimensional data, and difficulty with discontinuities or sharp kinks because the method relies on smooth interpolants and Hessians. The stated future directions are alternative surrogates such as KHRONOS, improved symbolic regression engines beyond PySR, and deeper integration of gradient and Hessian information into the symbolic discovery process. A plausible implication is that the framework is most effective when accurate smooth surrogates can be learned and when the target law possesses at least approximate separable structure.