NeuronSeek: On Stability and Expressivity of Task-Driven Neurons

1. Introduction

Deep learning has achieved significant success in a wide range of fields [1,2]. Credits for these accomplishments largely go to the so-called scaling law, wherein fundamental computational units such as convolutions and attention mechanisms are designed and scaled per an optimized architecture. Prominent examples include ResNet [3], Transformer [4], and Mamba [5]. Recently, motivated by neuroscience’s paramount contributions to AI, an emerging field called NeuroAI has garnered significant attention, which draws upon the principles of biological circuits in the human brain to catalyze the next revolution in AI [6]. The ambitious goal of NeuroAI is rooted in the belief that the human brain, as one of the most intelligent systems, inherently possesses the capacity to address complex challenges in AI development [7]. Hence, it can always serve as a valuable source of inspiration for guiding practitioners, despite the correspondence between biological and artificial neural networks sometimes being implicit.

Following NeuroAI, let us analyze the scaling law through the lens of brain computation. It is seen that the human brain generates complex intellectual behaviors through the cooperative activity of billions of mutually connected neurons with varied morphologies and functionalities [8]. This suggests that the brain simultaneously benefits from scale at the macroscopic level and neuronal diversity at the microscopic level. The latter is a natural consequence of stem cells’ directed programming in order to facilitate the task-specific information processing. Inspired by this observation, extensive research has shown that incorporating neuronal diversity into artificial neural networks can markedly enhance their capabilities [7]. For instance, researchers have transitioned beyond traditional inner-product neurons by integrating diverse nonlinear aggregation functions, such as quadratic neurons [9], polynomial neurons [10], and dendritic neurons [11]. Such networks have achieved superior performance in tasks like image recognition and generation, underscoring the promise and practicality of innovating a neural network at the neuronal level.

Along this direction, driven by the task-specific nature of neurons in our brain, recent research [12] introduced a systematic framework for task-driven neurons. Hereafter, we refer to prototyping task-driven neurons as NeuronSeek. This approach assumes that there is no single type of neurons that can perform universally well. Therefore, it enables neurons to be tailored for specific tasks through a two-stage process. As shown in Figure 1, in the first stage, the vectorized symbolic regression (SR)—an extension of symbolic regression [13,14] that enforces uniform operations for all variables—identifies the optimal homogeneous neuronal formula from data. Unlike traditional linear regression and symbolic regression, the vectorized SR searches for coefficients and homogeneous formulas at the same time, which can facilitate parallel computing. In the second stage, the derived formula serves as the neuron’s aggregation function, with learnable coefficients. Notably, activation functions remain unmodified in task-driven neurons. At the neuronal level, task-driven neurons are supposed to have superior performance than those one-size-fits-all units do, as it imbues the prior information pertinent to the task. More favorably, task-driven neurons can still be connected into a network to fully leverage the power of connectionism. Hereafter, we refer to task-driven neurons using symbolic regression as NeuronSeek-SR. NeuronSeek considers in tandem neuronal importance and scale. From design to implementation, it systematically integrates three cornerstone AI paradigms: behaviorism (the genetic algorithm runs symbolic regression), symbolism (symbolic regression distills a formula from data), and connectionism.

Despite this initial success, the following two major challenges remain unresolved in prototyping task-driven neurons.

• Non-deterministic and unstable convergence. While symbolic regression employs modified genetic programming (GP) to discover polynomial expressions [15,16], its tree-based evolutionary approach struggles to converge in high-dimensional space. The resulting formulas demonstrate acute sensitivity to certain hyperparameters, such as initialization conditions, population size, and mutation rate. This may confuse users about whether the identified formula is truly suitable for the task or due to the random perturbation, therefore hurting the legibility of task-driven neurons.
• Necessity of task-driven neurons. In the realm of deep learning theory, it was shown [17,18,19] that certain activation functions can empower a neural network using a fixed number of neurons to approximate any continuous function with an arbitrarily small error. These functions are termed “super-expressive" activation functions. Such a unique and desirable property allows a network to achieve a precise approximation without increasing structural complexity. Given the tremendous theoretical advantages of adjusting activation functions, it is necessary to address the following issue: Can the super-expressive property be achieved via revising the aggregation function while retaining common activation functions?

In this study, we address the above two issues satisfactorily. On one hand, in response to the instability of NeuronSeek-SR, we propose a tensor decomposition (TD) method referred to as NeuronSeek-TD. This approach begins by assuming that the optimal representation of data is encapsulated by a high-order polynomial coupled with trigonometric functions. We construct the basic formulation and apply TD to optimize its low-rank structure and coefficients. Therefore, we reformulate the unstable formula search problem (NeuronSeek-SR) into a stable low-rank coefficient optimization task. To enhance robustness, we introduce the sparsity regularization in the decomposition process, automatically eliminating insignificant terms while ensuring the framework to identify the optimal formula. Additionally, rank regularization is employed to derive the simplest possible representation. The improved stability of our method stems from two key factors: i) Unlike conventional symbolic regression which suffers from complex hyperparameter tuning in genetic programming [20], NeuronSeek-TD renders significantly fewer tunable parameters. ii) It has been proven that tensor decomposition with rank regularization tends to yield a unique solution [21], thereby effectively addressing the inconsistency issues inherent in symbolic regression methods.

On the other hand, we close the theoretical deficit by showing that task-driven neurons that use common activation such as ReLU can also achieve the super-expressive property. Earlier theories like [17,18,19] first turn the approximation problem into the point-fitting problem, and then uses the dense trajectory to realize the super-expressive property. The key message is the existence of a one-dimensional dense trajectory to cover the space of interest. In this study, we highlight that dense trajectory can also be achieved by adjusting aggregation functions. Specifically, we integrate task-driven neurons into a discrete dynamical system defined by the layer-wise transformation ${\mathbf{x}}_N=T\left({\mathbf{x}}_{N-1}\right)$, where T is the mapping performed by one layer of the task-driven network. There exists some initial point ${\mathbf{x}}_0$ that can be taken by T to the neighborhood of any target point $\mathbf{x}$. To approximate $\mathbf{x}$, our construction does not need to increase the values of parameters in T and the network, and we just compose the transformation module T different times. This means that we realize “super-super-expressiveness": not only the number of parameters but also the magnitudes remain fixed. In summary, our contributions are threefold:

• We introduce tensor decomposition to prototype stable task-driven neurons by enhancing the stability of the process of finding a suitable formula for the given task. Our study is a major modification to NeuronSeek.
• We theoretically show that task-driven neurons with common activation functions also enjoy the “super-super-expressive” property, which puts the task-driven neuron on a solid theoretical foundation. The novelty of our construction lies in that it fixes both the number and values of parameters, while the previous construction only fixes the number of parameters.
• Systematic experiments on public datasets and real-world applications demonstrate that NeuronSeek-TD achieves not only superior performance but also stability compared to NeuronSeek-SR. Moreover, its performance is also competitive compared to other state-of-the-art baselines.

2. Related Works

2.1. Neuronal Diversity and Task-Driven Neurons

Many groundbreaking AI advancements were inspired by the computing of biological neural systems. For instance, the neocognitron [22], a precursor to convolutional neural networks, is exactly a miniature of cortical visual circuits. The human brain is renowned for its extensive array of neurons, which vary significantly in morphology and functionality. This neuronal diversity underpins the brain’s intelligent behavior [23,24]. However, traditional artificial neural networks often prioritize the scaling law that designs a fundamental unit and replicates it across large architectures. This approach pays little attention to neuronal diversity. Recently, the idea of introducing neuronal diversity has naturally emerged, which mainly modifies the aggregation function instead of the activation function in constructing new neurons. This operation is well-grounded: Once an activation function is monotonic, a neuron’s decision boundary is exclusively shaped by its aggregation function. In particular, it replaces the inner product of a neuron with a nonlinear function such as a polynomial. Table 1 summarizes the recently proposed non-linear neurons. Notably, the complexity of neurons in [25,26,27] is of $\mathcal{O}\left(n^2\right)$, which is much larger than the conventional neuron, while neurons in [9,28,29,30] favorably enjoy the linear parametric complexity.

Following the way paved for task-specific neuronal design [12], this work constructs neurons by tensor decomposition. It simultaneously resolves theoretical deficits and improves the empirical stability of the original framework, representing a substantial advance in handling high-dimensional learning tasks like image recognition. While Chrysos et al. [31,32,33] focus on designing polynomial network architectures, our work differs by concentrating on task-specific neuronal construction and more diverse functional forms such as trigonometric functions ($sin\left(\mathbf{x}\right)$, $cos\left(\mathbf{x}\right)$).

2.2. Symbolic Regression and Deep Learning

Symbolic regression (SR) is an algorithm that searches for formulas from data [13]. Recent research in SR includes learning partial differential equations (PDEs) from data [34,35], improving the validity of SR on high-dimensional data [36], and invertible SR [37]. Besides, the Kolmogorov-Arnold Network (KAN) is an exemplary work that combines symbolic regression and deep learning [38]. KAN is inspired by the celebrated Kolmogorov-Arnold representation theorem, which states that any continuous multivariate function can be represented as a superposition of several continuous univariate functions. Such a decomposability lays KAN a solid theoretical foundation for function approximation.

Specifically, Kolmogorov’s theorem asserts that for any vector $\mathbf{x}={\left[x_1,x_2,\dots ,x_d\right]}^T\in {[0,1]}^d$, a continuous function $f\left(\mathbf{x}\right)$ can be expressed as

$$f(x_1,x_2,\dots ,x_d)=\sum _{i=0}^{2d}{g}_i\left(\sum _{j=1}^d{h}_{i,j}\left(x_j\right)\right),$$

(1)

where $g_i$ and $h_{i,j}$ are continuous univariate functions. Building upon this theorem, KAN is a neural network to exploit this decomposition by learning $g_i$ and $h_{i,j}$, which can be regarded as a special symbolic regression. Currently, this KAN-type line has been successfully used in a myraid of domains, such as medical image analysis [39] and scientific discovery [40].

Moreover, deep learning is also applied to strengthen the performance of SR algorithms. Kamienny et al. [41] combined the Transformer and symbolic regression to enhance its searching ability. Kim et al. [42] integrated the neural network with symbolic regression, and proposed a framework similar to the equation learner network. SymbolicGPT, a Transformer-based language model for symbolic regression [43], is also shown to be a competent model compared with the existing models in terms of accuracy and run time.

To summarize, symbolism is one of the earliest and most influential paradigms in the history of AI [44,45]. It is rooted in the idea that intelligence can be modeled through the manipulation of symbols—abstract representations of objects, concepts, or relationships in the real world. This paradigm assumes that reasoning and problem-solving can be achieved through symbolic manipulation. Symbolism and connectionism are highly complementary [46,47]. Incorporating techniques in symbolism can effectively solve the intrinsic problems in connectionism such as efficiency, robustness, and forgetting. Our work synergizes symbolic regression and deep learning at the neuronal level, which is a novel practice in fusing symbolism and connectionism.

2.3. Super-Expressive Activation Functions

These days, super-expressive activation functions gain lots of traction in deep learning approximation theory [48]. Unlike the traditional activation functions, super-expressive activation functions can empower a network to achieve an arbitrarily small error with a fixed network width and depth. The idea of super-expressive activation function dates back to [49], while the activation function therein was not in closed form. Recently, several studies proposed explicit super-expressive activation functions [17,18,19], such as $sin\left(x\right),arcsin\left(x\right)$, and $x-\mathtt{floor}\left(x\right)$ with $x>0$. The key to achieving super-expressiveness lies in leveraging the dense orbit to solve the point-fitting problem $(n,f_n)$, where $f_n\in [0,1],n=1,2,\dots ,N$. The fact that there exists $\theta$ such that $f_n=\theta /(\pi +n)-\mathtt{floor}(\theta /(\pi +n))$ was used in [18].

Our work shows that modifying the aggregation function can even endow a network with the super-super-expressive property. Particularly, our theoretical construction is much more efficient in memory. Previously, though the number of parameters like $\theta$ is fixed in [18], $\theta$ grows large when the approximation error goes low, and thus the essential network memory is still larger. In contrast, our construction maintains both a fixed number and magnitude of parameters, which brings a significant memory saving. Our work provides a strong theoretical foundation for task-driven neurons. Our proof is inspired by [50] that also composes a fixed block to do approximation. [50] uses the bit extract technique, while ours uses the chaotic theory.

3. Task-Driven Neuron Based on Tensor Decomposition (NeuronSeek-TD)

3.1. NeuronSeek-SR

In this section, we briefly illustrate how NeuronSeek-SR creates novel neurons. It follows a two-stage process, as illustrated in Figure 1. First, the vectorized SR identifies the optimal homogeneous neuronal formula by fitting task-specific data. Second, the derived formula is implemented as the neuron’s aggregation function with learnable coefficients.

3.1.1. Vectorized SR

Conventional SR encodes a mathematical expression as a tree-structured directed graph and utilizes GP to navigate the solution space [15]. The GP optimization process employs biologically-inspired evolutionary operators, including: i) crossover operations that exchange subtrees between candidate solutions to maintain population diversity, and ii) mutation operations that introduce random modifications to explore novel regions of the function space. While this approach offers notable flexibility for symbolic optimization, it faces significant challenges such as computational complexity, formula heterogeneity, and overfitting when applied to high-dimensional problems. All these together are not ideal for the task of prototyping new neuronal formulations.

Unlike the conventional SR, which discovers heterogeneous formulas for each variable, the vectorized SR regularizes all input variables to learn the same formula [13,14]. This single change brings huge advantages: i) Lower complexity. Homogeneous formulas generated by the vectorized SR have significantly fewer parameters, scaling linearly with the input dimensionality. ii) Faster computation and minor overfitting risks. By limiting the search space, the vectorized SR reduces the risk of overfitting and also makes the regression process faster and more computationally efficient, especially for high-dimensional data.

3.1.2. Parameterization

The learned formulas from the vectorized SR are parameterized, i.e., making the coefficients trainable. This allows neurons to fine-tune their behavior in the context of a network. While task-based neurons alone may have limited approximation power, connecting them into a neural network framework will further amplify the power of neurons, thereby enabling efficient optimization and improving performance in handling complex tasks.

Despite the success of the vectorized SR, its built-in tree-based evolutionary algorithm still suffers from poor convergence in high-dimensional space. Consequently, the derived formulas exhibit significant sensitivity to key hyperparameters—including initialization conditions, population size, and mutation rate—often producing inconsistent solutions across different runs. This inherent instability in the NeuronSeek-SR framework stems from the intrinsic issues of genetic algorithms.

3.2. NeuronSeek-TD

In this work, we propose NeuronSeek-TD, a faster and more stable alternative to conventional vectorized SR for neuronal formula discovery. We reasonably assert that the maximum required polynomial order for any given task is bounded by a constant N. By applying tensor decomposition to optimize coefficients, we reformulate the unstable formula search problem into a stable low-rank coefficient optimization task. Furthermore, we augment the polynomial structure through elementary functions, providing greater design flexibility and improved performance in high-dimensional pattern recognition tasks. The overall framework is illustrated in Figure 2, which comprises two stages: search for the neuronal formula and construct the task-driven neural network.

3.2.1. Neuronal Formula Discovery

In the first stage, we adopt a polynomial model to identify the expressions that best fit the task-specific data. For consistent processing across modalities, both tabular and image inputs are flattened into a vector representation $\mathbf{z}\in {\mathbb{R}}^d$. Then, given the input data $\mathbf{z}$, the initial neuronal formula is defined as follows:

$$\begin{array}{cc}\hfill f\left(\mathbf{z}\right)& =C_1·{\beta }_j+C_2·{\mathbf{z}}^{\top }{\mathbf{W}}^{\left[1\right]}+C_3·{\mathbf{z}}^{\top }{\mathbf{W}}^{\left[2\right]}\mathbf{z}\hfill \\ \hfill & +C_4·{\mathbf{W}}^{\left[3\right]}{\times }_1\mathbf{z}{\times }_2\mathbf{z}{\times }_3\mathbf{z}+\dots +C_N·{\mathbf{W}}^{\left[N\right]}\prod _{n=1}^N{\times }_n\mathbf{z}\hfill \\ \hfill & +D·sin\left(\mathbf{z}\right),\hfill \end{array}$$

(2)

where ${\mathbf{W}}^{\left[n\right]}$ are weights and $C_j$ are coefficients. To put it in a more compact way, we have

$$\begin{array}{cc}\hfill f\left(\mathbf{z}\right)& =\sum _{n=1}^N\left(C_n·{\mathcal{W}}^{\left[n\right]}\prod _{m=1}^n\mathbf{z}\right)+D·sin\left(\mathbf{z}\right),\hfill \end{array}$$

(3)

where the order N and the rank are two hyperparameters to be set appropriately. The main idea of this method is to outer-product a series of rank-1 tensors to represent high-order tensors so that the computation cost is reduced significantly. The tensor decomposition process could be formulated as follows:

$${\mathbf{W}}^{\left[n\right]}=\sum _{r=1}^R{\mathbf{a}}_r^{\left(1\right)}\circ {\mathbf{a}}_r^{\left(2\right)}\circ \dots \circ {\mathbf{a}}_r^{\left(n\right)},$$

(4)

where R is decomposition rank, and n is the order of the decomposed tensor $\mathbf{W}$. Here, we use COMP/PARAFAC (CP) decomposition to execute this process. CP decomposition has been widely adopted for modeling multilinear relationships in machine learning applications [31,51]. In contrast to Tucker decomposition, CP decomposition provides both computational and structural advantages for our neuronal formula discovery framework. Computationally, CP decomposition requires only $O\left(nRd\right)$ parameters for an n-th order tensor of size $d^n$, compared to Tucker’s $O(R^n+nRd)$ complexity, enabling efficient handling of high-order polynomial terms. An example of CP decomposition is illustrated in Figure 4. While other tensor decomposition approaches could be explored (e.g., tensor ring decomposition, hierarchical Tucker decomposition) [52,53], the efficiency and structural benefits of CP decomposition have already aligned well with the requirements of our polynomial modeling framework.

Figure 3. A rank-R CP decomposition of a third-order tensor.

Figure 3. A rank-R CP decomposition of a third-order tensor.

Moreover, the incorporation of the trigonometric term $Dsin\left(\mathbf{z}\right)$ serves a pivotal purpose. It extends the model’s capacity to capture periodic patterns and high-frequency components that intrinsically appear in many real-world datasets [54]. This is particularly crucial as pure polynomial expansions exhibit limited spectral sensitivity.

We then impose sparsity constraints through $L_1$ regularization on both the coefficients $C_i$ and weight tensors ${\mathbf{W}}^{\left[i\right]}$ to promote computational efficiency. Specifically, since each tensor ${\mathbf{W}}^{\left[i\right]}$ is decomposed using CP decomposition as ${\mathbf{W}}^{\left[i\right]}={\sum }_{r=1}^R{\mathbf{a}}_r^{\left(1\right)}\circ {\mathbf{a}}_r^{\left(2\right)}\circ \dots \circ {\mathbf{a}}_r^{\left(i\right)}$, we apply sparsity regularization directly to the factor matrices. Thus, the complete optimization objective becomes

$$\underset{\left\{{\mathbf{a}}_r^{\left(j\right)}\right\},C_i}{min}\sum _{k=1}^K\mathcal{L}(f\left({\mathbf{x}}_k\right),y_k)+{\lambda }_1{\parallel \mathbf{C}\parallel }_1+{\lambda }_2\sum _{i=1}^N\sum _{r=1}^R\sum _{j=1}^i{\parallel {\mathbf{a}}_r^{\left(j\right)}\parallel }_1,$$

(5)

where $\mathcal{L}(·)$ is the task-specific loss function, set to mean squared error (MSE) for regression and cross-entropy (CE) loss for classification. This structured sparsity approach encourages the discovery of interpretable and simplified polynomial formulas while maintaining computational efficiency through the low-rank tensor representation.

To further enhance the compactness of the discovered neuronal formula beyond the $L_1$ regularization’s sparsity induction, we introduce a post-training thresholding mechanism for explicit term selection. While the $L_1$ regularization effectively drives many coefficients toward zero, empirical observations reveal that a non-negligible portion of coefficients $C_i$ retain small but non-zero values, which cumulatively increase model complexity without substantial predictive contribution. The thresholding mechanism is a statistically grounded pruning strategy: After model convergence, we compute the mean ${\mu }_{\left|C\right|}$ and the standard deviation ${\sigma }_{\left|C\right|}$ of the absolute values of all coefficients $\left\{\right|C_i{\left|\right\}}_{i=1}^N$, then construct the final neuronal formula by retaining only terms whose coefficients satisfy $|C_i|\ge \tau$, where $\tau ={\mu }_{\left|C\right|}-{\sigma }_{\left|C\right|}$. This can further increase the stability of identified formulas.

Remark 1. Three key features of the proposed method warrant emphasis. i) The higher-order polynomial architecture systematically enables the learning of the polynomials. It is noticed that Eq. (2) is rather inclusive, e.g., it identifies multi-order variable interactions. This structural flexibility facilitates sophisticated nonlinear approximations, aligning the model more closely with the underlying data distribution. ii) The integration of trainable coefficients $\mathbf{C}$ with the $L_1$ regularization imposes selective sparsity constraints on the polynomial expansion. Specifically, the regularization term $\mathbf{C}$ induces sparsity by driving coefficients of less predictive terms toward zero, thereby automatically pruning nonessential higher-order components and yielding a parsimonious architecture. iii) The complementary $L_1$ regularization on the coefficients $\mathbf{W}$ enforces structured sparsity in the parameter space. This dual regularization strategy enhances feature selection while simultaneously improving computational efficiency—a critical advantage in high-dimensional contexts. Collectively, these design choices balance expressive capacity with compactness, ensuring robust performance without overparameterization.

3.2.2. Task-Driven Network Construction

In the second stage, the discovered neuronal formula is first parameterized into neuron-like structures. The selected polynomial and trigonometric terms are assigned trainable weights and integrated into neural network architectures. For tabular data processing, we employ fully connected neural networks. A polynomial with the maximum polynomial order of 3 is expressed as

$$y=\sigma ({\mathbf{W}}_1^{\top }\mathbf{z}+{\mathbf{W}}_2^{\top }\mathbf{z}\odot \mathbf{z}+{\mathbf{W}}_3^{\top }\mathbf{z}\odot \mathbf{z}\odot \mathbf{z}+b),$$

(6)

with $\sigma$ being the ReLU activation function. If a specific order term is absent in the neuronal formula, its corresponding weight and bias parameters are initialized to zero.

Moreover, we extend our framework to convolutional neurons for images. As mentioned, the trigonometric term is incorporated to model periodic patterns and high-frequency components inherent in high-dimensional images. Consequently, the complete neuronal formula is structured as follows:

$$y=\sigma ({\mathbf{W}}_1^{\top }*\mathbf{z}+{\mathbf{W}}_2^{\top }*\mathbf{z}\odot \mathbf{z}+{\mathbf{W}}_3^{\top }*\mathbf{z}\odot \mathbf{z}\odot \mathbf{z}+{\mathbf{W}}_4*sin\left(\mathbf{z}\right)+b),$$

(7)

where * denotes the convolution operation.

Finally, we introduce a simple initialization strategy to enhance training stability. Followed by [33], the weight matrix ${\mathbf{W}}_1$ is initialized as ${\mathbf{W}}_1\sim \mathcal{N}(0,{\sigma }^2)$, where $\sigma =\sqrt{1/32\left(k_n\right)}$, while all other weights are initialized to a tiny number such as 0.001. This strategy ensures that the network initially behaves like linear neurons and gradually learns higher-order terms during training. This strategy has been demonstrated the effectiveness in training quadratic neurons [55].

4. Super-Expressive Property Theory of Task-Driven Neuron

In this section, we formally prove that task-driven neurons that modify the aggregation function can also achieve the super-expressive property. Our theory is established based on the chaotic theory. Mathematically, we have the following theorem:

Theorem 1.

Let $f\in C\left(\right[0,1\left]\right)$ be a continuous function. Then, for any $\epsilon >0$, there exists a function h generated by a task-driven network with a fixed number of parameters, such that

$$\left|f\right(x)-h(x\left)\right|<\epsilon \mathrm{for}\mathrm{any}x\in [0,1].$$

(8)

Proof sketch: Our proof leverages heavily the framework of [18] that turns via the interval partition the function approximation into the point-fitting problem. Hence, for consistency and readability, we inherit their notations and framework. The proof idea is divided into three steps:

Step 1. As Figure 4 shows, in equal distance, divide $[0,1)$ into small intervals ${\mathcal{J}}_k=\left[{\displaystyle \frac{k-1}{K}},{\displaystyle \frac{k}{K}}\right)$ for $k\in \{1,2,\dots ,K\}$, where K is the number of intervals. The left endpoint of ${\mathcal{J}}_k$ is $x_k$. A higher K leads to a smaller approximation error. We construct a piecewise constant function h to approximate f. The error can be arbitrarily small as long as the divided interval goes tiny. Based on the interval division, we have

$$h\left(x\right)\approx f\left({\tilde{x}}_k\right),x\in {\mathcal{J}}_k,$$

(9)

where ${\tilde{x}}_k$ is a point from ${\mathcal{J}}_k$.

Step 2. As Figure 4 shows, use the function ${\varphi }_1\left(x\right)=\lfloor Kx\rfloor$ mapping all x in the interval ${\mathcal{J}}_k$ to k, where $\lfloor ·\rfloor$ is the flooring function. There exists a one-to-one correspondence between k and $f\left({\tilde{x}}_k\right)$. Thus, via dividing intervals and applying the floor function, the problem of approximation is simplified into a point-fitting problem. Therefore, we just need to construct a point-fitting function to map k to $f\left({\tilde{x}}_k\right)$.

Step 3. In [18], the point-fitting function is constructed as $\tau (\theta /(\pi +n\left)\right)$ based on Proposition 1. In contrast, we solve it with a discrete dynamic system that is constructed by a network made of task-driven neurons. We design a sub-network to generate a function ${\varphi }_2$ mapping k approximately to $f\left({\tilde{x}}_k\right)$ for each k. Then ${\varphi }_2\left({\varphi }_1\left(x\right)\right)={\varphi }_2\left(k\right)\approx f\left({\tilde{x}}_k\right)\approx f\left(x\right)$ for any $x\in {\mathcal{J}}_k$ and $k\in \{1,2,\dots ,K\}$, which implies ${\varphi }_2\circ {\varphi }_1\approx f$ on $[0,1)$. ${\varphi }_2$ is $T^{m\left(k\right)}\left(\xi \right)$ based on Lemma 1, where T is a unimodal function from $[0,1]\to [0,1]$ that can also induce the dense trajectory.

Proposition 1

([56]). For any given set ${\left\{w_n\right\}}_{n=1}^N\subseteq [0,1]$ and $ϵ>0$, there exists a value ${\theta }^{*}\in [0,+\infty )$, such that

$$|w_n-\tau \left({\displaystyle \frac{{\theta }^{*}}{\pi +n}}\right)|<ϵ,n=1,\dots ,N,$$

(10)

where $\tau \left(z\right)=z-\lfloor z\rfloor$.

Definition 1

([57]). A map $T:[0,1]\to [0,1]$ is said to be unimodal if there exists a turning point $\eta \in [0,1]$ such that the map f can be expressed as

$$T\left(x\right)=min\left\{T_{\mathrm{L}}\left(x\right),T_{\mathrm{R}}\left(x\right)\right\}=\left\{\begin{array}{c}{T}_{\mathrm{L}}\left(x\right),0⩽x⩽\eta ,\hfill \\ T_{\mathrm{R}}\left(x\right),\eta ⩽x⩽1,\hfill \end{array}\right.$$

(11)

where $T_{\mathrm{L}}:[0,\eta ]\to [0,1]$ and $T_{\mathrm{R}}:[\eta ,1]\to [0,1]$ are continuous, differentiable except possibly at finite points, monotonically increasing and decreasing, respectively, and onto the unit-interval in the sense that $T_{\mathrm{L}}\left(0\right)=T_{\mathrm{R}}\left(1\right)=0$ and $T_{\mathrm{L}}\left(\eta \right)=T_{\mathrm{R}}\left(\eta \right)=1$.

Lemma 1

(Dense Trajectory). There exists a measure-preserving transformation T, generated from a network with task-driven neurons, such that there exists ξ with $\{T^n\left(\xi \right):n=\mathbb{N}\}$ dense in $[0,1]$. In other words, for any x and $ϵ>0$, there exist ξ and the corresponding composition times m, such that

$$|x-T^m\left(\xi \right)|<ϵ.$$

(12)

Particularly, we refer to ξ as a transitive point.

Proof. 

Since we extensively apply polynomials to prototyping task-driven neurons, we prove that a polynomial of any order can fulfill the condition. The proof can be easily extended to other functions with minor twists.

i) Given a polynomial of any order, as shown in Figure 5, through cutting, translation, flipping, and scaling, we can construct a unimodal map $T\left(x\right):[0,1]\to [0,1]$.

ii) For a unimodal map T, we can show that for any two non-empty sets $U,V\in X$ with $\mathfrak{m}\left(U\right),\mathfrak{m}\left(V\right)>0$, there exists an m such that $T^m\left(U\right)\cap V\ne ⌀$. As Figure 6 shows, this is because a set U with $\mathfrak{m}\left(U\right)>0$ will always contain a tiny subset ${\mathcal{J}}_{mi}$ such that $T^m\left({\mathcal{J}}_{mi}\right)=[0,1]$. Then, naturally $T^m\left(U\right)\cap V\ne ⌀$ due to ${\mathcal{J}}_{mi}\subset U$.

iii) The above property can lead to a dense orbit, combining that $[0,1]$ is with a countable basis, according to Proposition 2 of [58]. For the sake of self-sufficiency, we include their proof here.

Let ${\left(V_i\right)}_{i\in I}$ be a countable base for X. For $i\in I$, the set $W_i={\bigcup }_{n\geqq 0}{f}^{-n}\left(V_i\right)$ is open by continuity of f. This set is also dense in X. To see this let U be any non-empty open set. Because of topological transitivity there exists a $k>0$ with $f^k\left(U\right)\cap V_i\ne \varnothing$. This gives $f^{-k}\left(V_i\right)\cap U\ne \varnothing$ and $W_i\cap U\ne \varnothing$. Thus $W_i$ is dense. By the Baire category theorem, the set $B={\cap }_{i\in I}{W}_i$ is dense in X. Now, the orbit of any point $\xi \in B$ is dense in X.

Because, given any non-empty open $U\subset X$, there is an $i\in I$ with $V_i\subset U$ and $k>0$ with $\xi \in f^{-k}\left(V_i\right)$ . This means $f^k\left(\xi \right)\in V_i\subset U$. Thus the orbit of $\xi$ enters any U.

□

Proof of Theorem 1. Now, we are ready to prove Theorem 1. First, the input x is mapped into an integer k, whereas k corresponds to the composition time $m\left(k\right)$. Next, a network module representing T is copied $m\left(k\right)$ times and composed $m\left(k\right)$ times such that $T^{m\left(k\right)}\left(\xi \right)\approx f\left({\tilde{x}}_k\right)$. Since $f\left({\tilde{x}}_k\right)$ can approximate $f\left(x\right)$ well, $T^{m\left(k\right)}\left(\xi \right)$ also approximates $f\left(x\right)$ well.

Theorem 1 can be easily generalized into multivariate functions by using the Kolmogorov-Arnold theorem.

Proposition 2.

For any $\mathbf{x}={\left[x_1,x_2,\dots ,x_d\right]}^T\in {[0,1]}^d$, there exist continuous functions $h_{i,j}\in C\left([0,1]\right)$ for $i=0,1,\dots ,2d$ and $j=1,2,\dots ,d$ such that any continuous function $f\in C\left({[0,1]}^d\right)$ can be represented as

$$f\left(\mathbf{x}\right)=\sum _{i=0}^{2d}{\Psi }_i\left(\sum _{j=1}^d{\Phi }_{i,j}\left(x_j\right)\right).$$

(13)

where ${\Psi }_i:\mathbb{R}\to \mathbb{R}$ is a continuous function for each $i\in \{0,1,\dots ,2d\}$.

Remark 2. Our construction saves the parametric complexity. In [18], though the number of parameters remains unchanged, the values of parameters ${\theta }^{*}$ increase as long as n increases. Therefore, the entire parametric complexity goes up. In contrast, in our construction, the parameter value of T is also fixed for different $(k,f\left({\tilde{x}}_k\right))$. Yet, T is copied and composed $m\left(k\right)$ times. The overall parametric complexity remains unchanged. What increases herein is computational complexity. Therefore, we call for brevity our construction “super-super-expressive”.

5. Comparative Experiments on NeuronSeek

This section presents a comprehensive experimental evaluation of our proposed method and the previous NeuronSeek-SR [12]. First, we systematically demonstrate the superior performance and enhanced stability of NeuronSeek-TD compared to the previous NeuronSeek-SR. The results demonstrate that pursuing stability in NeuronSeek-TD does not hurt the representation power of the resultant task-driven neurons and even present better results.

5.1. Experiment on Synthetic Data

Here, we adopt the synthetic data to compare the performance of NeuronSeek-TD and NeuronSeek-SR. The synthetic data are generated by adding noise perturbations to the known ground-truth functions. We then measure the regression error (MSE) of the formulas of NeuronSeek-SR and NeuronSeek-TD on the synthetic data.

5.1.1. Experimental Settings

Eight multivariate polynomials are constructed as shown in Table 2. For simplicity and conciseness, we prescribe:

$$\sum _{i=1}^{d}r{x}_i^s=r·\underset{s\mathrm{times}}{\underbrace{{(\mathbf{x}\odot \dots \odot \mathbf{x})}^{\top }}}=r·{\left(\mathbf{x}{\odot }^s\mathbf{x}\right)}^{\top },$$

where $\mathbf{r}=[r,\dots ,r]\in {\mathbb{R}}^d$ denotes a constant vector, $\mathbf{x}=[x_1,\dots ,x_d]$ represents the input vector, and ⊙ denotes the Hadamard product. The elements of $\mathbf{x}$ are uniformly sampled from the interval $[-1,1]$. We set $d=30$ to simulate multi-dimensional inputs.

To evaluate the noise robustness of the proposed methods, we introduce three noise levels through

$$\overline{y}=y+\delta ·\mathrm{RMS}\left(y\right)·\mathit{\eta },$$

where y denotes the noise-free polynomial output, $\overline{y}$ is the noisy observation, $\mathit{\eta }\sim \mathcal{N}(\mathbf{0},\mathit{I})$ is a standard normal random vector, $\delta$ controls the relative noise level, and $\mathrm{RMS}\left(y\right)$ represents the root mean square (RMS) of the data:

$$\mathrm{RMS}\left(y\right)=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{y}_i^2}.$$

Here, $N=1,000$ specifies the number of data points. For NeuronSeek-SR, we configure a population size of 1,000 with the maximum of 50 generations. NeuronSeek-TD is optimized using Adam for 50 epochs. The noise levels are set to $\delta \in \{0\%,10\%,20\%\}$ respectively. We quantify the regression performance using the mean squared error (MSE) of the discovered neuronal formulae. Lower MSE values indicate better preservation of the underlying data characteristics.

5.1.2. Experimental Results

The performance comparison on the synthetic data is presented in Table 2. Our proposed NeuronSeek-TD demonstrates consistent superiority over NeuronSeek-SR across all generated data, with particularly notable improvements on the data derived from the function $\left(\mathbf{x}{\odot }^2\mathbf{x}\right)-1.5(\mathbf{x}\odot \mathbf{x})+2$. This performance advantage remains robust under varying noise levels.

To assess computational efficiency, we measure the average execution time for both methods. NeuronSeek-TD achieves substantial computational savings, requiring less than half the execution time of the SR-based approach while maintaining comparable solution quality. These results collectively demonstrate that our method achieves both superior accuracy and greater computational efficiency.

5.2. Experiments on Uniqueness and Stability

Here, we demonstrate that NeuronSeek-SR exhibits instability in terms of resulting in divergent formulas under minor perturbations in initialization, while NeuronSeek-TD maintains good consistency.

5.2.1. Experimental Settings

Our experimental protocol consists of two benchmark datasets: the phoneme dataset (a classification task classifying nasal and oral sounds) and the airfoil self-noise dataset (a regression task modeling self-generated noise of airfoil). These public datasets are chosen from the OpenML website to represent both classification and regression tasks, ensuring a comprehensive evaluation across different problem types. For each dataset, we introduce four levels of Gaussian noise ($\sigma \in \{0.01,0.025,0.05,0.1\}$) to simulate real-world perturbations. We train 10 independent instances of each method under identical initialization. For comparison, NeuronSeek-SR based on GP uses a population of 3,000 candidates evolved over 20 generations, whereas NeuronSeek-TD methods—employing CP decomposition and Tucker decomposition—rely on a model optimized via Adam over 20 epochs.

To quantify the uniqueness of SR and TD methods, we track two diversity metrics:

• Epoch-wise diversity represents the number of different formulas within each epoch (or generation) across all initializations. For instance, $x^3+2x^2$ and $2x^3+x^2$ are considered identical, whereas $x^3+x^2$ and $x^3+x$ are distinct.
• Cumulative diversity measures the average number of unique formulas identified per initialization as the optimization goes on.

5.2.2. Experimental Results

The results of phoneme and airfoil self-noise datasets are depicted in Figure 8 and Figure 7, respectively. From the epoch-wise diversity (Figure 8 (a) and Figure 7 (a)), we have several important observations: Firstly, GP-based SR demonstrates a rapid increase in the number of discovered formulas within each epoch. After training, it produces more unique formulas compared to TD-based methods (CP and Tucker). This indicates that SR is sensitive to initialization conditions. Secondly, as noise levels increase, the GP method discovers nearly 10 distinct formulas across 10 initializations. With higher noise, every initialization produces a completely different formula structure. These behaviors highlight the instability of GP-based SR when dealing with noisy data, as minor perturbations in initialization or data noise can lead to substantially different symbolic expressions.

Furthermore, regarding cumulative diversity (Figure 8(b) and Figure 7(b)), we can observe that GP-based SR constantly discovers an increasing number of different formulas, while TD methods tend to stabilize during training. This indicates that the GP algorithm struggles to converge to a consistent symbolic formula. In contrast, the asymptotic behavior exhibited by TD-based methods demonstrates their superior stability and convergence properties when identifying symbolic expressions. Moreover, the CP-based method suggests the best resilience against initialization variability.

5.3. Superiority of NeuronSeek-TD over NeuronSeek-SR

To validate the superior performance of NeuronSeek-TD over NeuronSeek-SR, we conduct experiments across 16 benchmark datasets—8 for regression and 8 for classification—selected from scikit-learn and OpenML repositories. First, we use NeuronSeek-TD and NeuronSeek-SR simultaneously to search for the optimal formula. Second, we construct fully connected networks (FCNs) with neurons taking these discovered formulas as their aggregation functions. To verify the ability of NeuronSeek-TD, the structures of neural networks are carefully set, such that NeuronSeek-TD has equal or fewer parameters than NeuronSeek-SR.

Comprehensive implementation details and test results are provided in Table 3. For regression tasks, the dataset is partitioned into training and test sets in an 8:2 ratio. The activation function is $\mathtt{ReLU}$ with MSE as the loss function and $\mathtt{RMSProp}$ as the optimizer. For classification tasks, the same train-test split is applied, using the $\mathtt{Sigmoid}$ activation function. The loss function and optimizer are set to CE and $\mathtt{Adam}$, respectively. All experiments are repeated 10 times and reported by mean and standard deviation, with MSE for regression and accuracy for classification.

The regression results in Table 3 demonstrate that NeuronSeek-SR exhibits larger fitting errors compared to NeuronSeek-TD, along with higher standard deviations in MSE values. A notable example is the airfoil self-noise dataset, where NeuronSeek-TD outperforms NeuronSeek-SR by a significant margin, achieving a 9% reduction in MSE. Furthermore, NeuronSeek-TD shows superior performance on both the California housing and abalone datasets while utilizing fewer parameters. These results collectively indicate that NeuronSeek-TD possesses stronger regression capabilities than NeuronSeek-SR.

This performance advantage extends to classification tasks, where NeuronSeek-TD consistently achieves higher accuracy than NeuronSeek-SR. Particularly in the HELOC and orange-vs-grapefruit datasets, NeuronSeek-TD attains better classification performance with greater parameter efficiency. The observed improvements across both regression and classification tasks suggest that NeuronSeek-TD offers superior modeling capabilities compared to NeuronSeek-SR.

5.4. Comparison of Tensor Decomposition Methods

Our framework employs the CP decomposition for neuronal discovery. To empirically validate this choice, we conduct a comparative analysis of some tensor decomposition methods used in NeuronSeek-TD.

5.4.1. Experimental Settings

Generally, we follow the previous experimental settings in Section 5.3: 16 tabular datasets and the same network structures as described in Table 3. We adopt CP decomposition and Tucker decomposition respectively for tensor processing in NeuronSeek-TD. We employ MSE for regression tasks, while classification accuracy is utilized for classification tasks.

5.4.2. Experimental results

First, Figure 9 illustrates that, for most datasets, the two decomposition methods exhibit comparable performance. Exceptions include airfoil Self-Noise and electricity datasets, where the two methods show slight performance differences. This observation suggests that variations in TD methods do not substantially impact the final results. Second, as shown in Figure 10, the CP-based method achieves a 42% reduction in computational time. This result verifies that the CP method offers a computational advantage over the Tucker method.

6. Competitive Experiments of NeuronSeek-TD over Other Baselines

6.1. Experiment on Tabular Data

To further illustrate the superiority of NeuronSeek-TD, we compare it with other state-of-the-art models over 4 real-world tabular tasks: Electron Collision Prediction1, Asteroid Prediction2, Heart Failure Detection3, and Stellar Classification4). We adopt 7 advanced machine learning methods, namely XGBoost [59], LightGBM [60], CatBoost [61], TabNet [62], TabTransformer [63], FT-Transformer [64], DANETs [65] and NeuronSeek-SR [12] as the baselines to compare with NeuronSeek-TD. Moreover, our proposed model and NeuronSeek-SR replace standard neurons with what are discovered by NeuronSeek-TD and NeuronSeek-SR, respectively.

All datasets are randomly split into training, validation, and test sets in a 0.8:0.1:0.1 ratio. For regression tasks (electron collision and asteroid prediction), we use the MSE as an evaluation metric, while for classification tasks with imbalanced labels (heart failure detection and stellar classification), we adopt the F1 score. Detailed experiment information is shown in Table 4.

Results are compared in Table 5. It can be seen that TabTransformer demonstrates competitive performance on the particle collision dataset. However, its effectiveness sharply declines on the asteroid prediction. In contrast, CatBoost and TabNet exhibit more consistent yet suboptimal results across four datasets. Lastly, the proposed NeuronSeek-TD method is the best performer, achieving the lowest error rates and highest F1 scores uniformly across all four benchmarks. Notably, NeuronSeek-TD significantly outperforms FT-Transformer and DANET. This significant margin underscores the efficacy of its task-adaptive architecture in recognizing feature interactions without domain-specific tuning.

6.2. Experiments on Images

To evaluate the performance enhancement of NeuronSeek-TD in computer vision tasks, we integrate it into several established convolutional neural network architectures: ResNet [3], DenseNet [66], SeResNet [67], and GoogLeNet [68]. Notably, we replace the first convolutional layer of each network block with a TN convolutional layer employing the discovered formula $x+sinx$. It should be noted that the original NeuronSeek-SR was not designed for high-dimensional image tasks due to convergence difficulties. Our implementation focuses on NeuronSeek-TD and follows the standardized training protocol from a popular image classification repository5, incorporating several advanced training techniques: Learning rate scheduling via $\mathtt{WarmUpLR}$ and $\mathtt{MultistepLR}$, regularization through $\mathtt{LabelSmoothing}$, and optimization using $\mathtt{SGD}$. We use grid search to fine-tune the learning rate from 0.5 to 0.05. All models are trained for 200 epochs with a batch size of 128. We evaluate the results using classification accuracy averaged over 10 independent runs.

We adopt four standard image datasets to validate the model: CIFAR 10/100 [69], SVHN [70], and STL10 [71]. The training and validation sets are split with a ratio of 0.9 : 0.1. In particular, images from the STL10 dataset are resized to $32\times 32$ to fit the input of the model. The classification results are summarized in Table 6.

Our experiments demonstrate that task-driven neuronal layers consistently enhance the baseline performance across multiple datasets. Notably, NeuronSeek-augmented ResNet18 achieves a 2% accuracy improvement on both CIFAR 100 and STL10 datasets. While SeResNet exhibits median performance on CIFAR and SVHN datasets, its performance degrades on STL10; however, TN layers compensate for this degradation by boosting its accuracy by 10%. These results suggest that task-driven neuronal formula discovery can effectively unlock the latent potential of fundamental network architectures.

7. Discussions and Conclusions

In this paper, we have proposed a novel method in the framework of NeuronSeek to enhance the capabilities of a neural network. Compared to neural architecture search (NAS) [72], NeuronSeek provides a different approach. Rather than modifying the overall model structure, NeuronSeek optimizes solely at the neuronal level by searching for optimal aggregation formulas. Firstly, we have introduced tensor decomposition to prototype stable task-driven neurons. Secondly, we have provided theoretical guarantees demonstrating that these neurons can maintain the super-super-expressive property. At last, through extensive empirical evaluation, we show that our proposed framework consistently outperforms baseline models while achieving the state-of-the-art results across diverse datasets and tasks. These advancements establish task-driven neurons as an effective addition in the design of task-specific neural network components, offering both theoretical soundness and practical efficacy.