STEP: Separable Tensor Polynomial Layers For Sustainable Deep Learning

1. Introduction

Deep neural networks (DNNs) are now one of the best tools for scientific computing, uncertainty quantification (UQ), and finding the numerical solution of partial differential equations (PDEs) [3,4]. They are very good at approximating high-dimensional nonlinear functions [5]. Their empirical success has been paralleled by a growing corpus of theoretical research that offers robust approximation assurances and elucidates how deep architectures can surmount the conventional curse of dimensionality that constrains classical numerical techniques. A significant area of research focuses on the expressive capacity of deep neural networks (DNNs) for parametric and countably-parametric partial differential equations (PDEs) [6,7]. In a groundbreaking study by Yarotsky et al., authors showed that deep ReLU networks can achieve dimension-independent approximation rates for generalized polynomial chaos (gPC) representations of holomorphic solution maps [8]. This proves that neural networks can effectively approximate complex models used in engineering and the physical sciences. These concepts were subsequently expanded to Bayesian inverse problems, such as Normalizing flow, wherein deep architectures yield dimension-independent approximation rates for posterior densities [9,10,11]

Figure 1. Schematic of the STEP architecture. Each component of the input $\mathbf{x}=(x_1,\dots ,x_d)$ is expanded in the Chebyshev basis $T_0,\dots ,T_N$. For each rank component $r=1,\dots ,R$, the layer forms one univariate polynomial per input variable and multiplies these polynomials across the d dimensions. The output is the sum of the R products, weighted by the scalars ${\lambda }_r$.

Figure 1. Schematic of the STEP architecture. Each component of the input $\mathbf{x}=(x_1,\dots ,x_d)$ is expanded in the Chebyshev basis $T_0,\dots ,T_N$. For each rank component $r=1,\dots ,R$, the layer forms one univariate polynomial per input variable and multiplies these polynomials across the d dimensions. The output is the sum of the R products, weighted by the scalars ${\lambda }_r$.

From the perspective of approximation theory, it is now firmly established that depth fundamentally enhances expressive efficiency. Yarotsky [8] demonstrated that deep ReLU networks approximate smooth and continuous functions with near-optimal complexity, while the analyses in [12] established precise lower bounds and optimality results for deep and sparsely connected architectures. Later research showed that deep networks can match or beat the approximation rates of traditional affine representation systems like wavelets, curvelets, and shearlets [13,14]. This shows that they are good for structured function classes.

These advancements in approximation theory have significantly influenced the numerical resolution of high-dimensional partial differential equations (PDEs). A series of studies by Hutzenthaler, Jentzen, Grohs, and their associates [13,14], established that DNN-based solvers surmount the curse of dimensionality for semilinear Kolmogorov equations, Hamilton–Jacobi–Bellman equations, Black–Scholes PDEs, and nonlinear stochastic control problems. In these contexts, the quantity of network parameters necessary to attain a specified accuracy increases only polynomially with respect to the dimension and the inverse error tolerance, signifying a substantial enhancement over traditional discretization methods.

More theoretical work has been done, such as full training error analyses for stochastic gradient descent [15], space-time error estimates for deep PDE solvers [16], and the design of specialized architectures that take advantage of the structural properties of the problems at hand. Convolutional and tensorized networks [17,18,19] for parametric diffusion and transport equations [20] and links between high-dimensional quadrature strategies and neural approximation rates [21] are two examples. Polynomial-based spectral concepts have also impacted practical deep learning models; for example, [22] presented a Chebyshev-based convolutional neural network for lung cancer CT analysis, showcasing enhanced performance compared to conventional CNNs.

Even with these improvements, approximating functions in high dimensions remains challenging. Classical multivariate polynomial expansions require ${(N+1)}^d$ coefficients while fully connected neural layers require $\Theta \left(dW\right)$ parameters (where N stands for the polynomial degree, d stands for data dimension and W stands for the network width). As a result, both methods become impractical in moderate to high dimensions. This drives the development of structured representations that keep their expressive power while significantly reducing the number of parameters. Related efforts have revisited the Kolmogorov–Arnold representation theorem, replacing fixed activations with learnable spline- or Chebyshev-parametrized edge functions [23,24], though these follow a different architectural design than the tensorized layer proposed here.

We propose a low-rank tensorized polynomial framework that yields explicit, scalable constructions compatible with modern deep architectures. Our method uses a simple, scalable architecture that fits well with modern deep learning design while dramatically reducing parameter complexity and still has strong theoretical guarantees offering a principled alternative to both dense neural layers and classical high-dimensional polynomial expansions. Table 1 shows a detailed comparison of representative orthogonal and structured convolutional layers highlighting the differences between the current methods and the proposed approach.

To deal with the computational and theoretical problems of current high-dimensional approximation methods, we present the Separable Tensor Polynomial (STEP) layer. This framework integrates low-rank tensor decomposition with Chebyshev spectral expansions to establish a Green AI [29] paradigm: attaining high-fidelity representations with significantly diminished depth, parameters, and energy consumption. The principal contributions of this study are encapsulated as follows:

• High-Efficiency "Green" Architecture: STEP uses a mode-separable structure that makes the number of parameters much simpler, going from the exponential $\mathcal{O}\left(N^d\right)$ of classical expansions to the linear $\mathcal{O}(R·d·N)$. This lets STEP work as a "shallow-yet-deep" model, achieving accuracy that is competitive with deep networks while using up to $1000\times$ fewer parameters. STEP captures global non-linear dependencies in a single layer by directly adding Chebyshev spectral representations to the neural activation path. This gets around the need for extreme depth and effectively reduces the "Curse of Dimensionality" in high-dimensional function approximation.
• Differentiable Tensor Interactions: We create a layer that can be differentiated natively and that does rank-separable interactions between spatial and spectral modes. This makes it possible to easily add STEP to modern architectures (like CNNs and dense networks) while keeping the theoretical rigor of low-rank tensor approximations. STEP offers verifiable convergence assurances that surpass the functionalities of conventional sparse grids and dense deep neural networks (DNNs). A comparative evaluation (see Table 1) shows that STEP is the only method that combines scalability, theoretical soundness, and compatibility with optimization based on gradients.

This work goes beyond the STEP as an architecture for shallow learning but rather to create a single mathematical framework for tensorized interactions that improves existing methods like sparse grids and polynomial expansions.

The remainder of this paper is organized as follows. In Section 2, we present the mathematical foundations of STEP, including its forward and backward passes as well as the spectral constructions. Section 3 focuses on theoretical evaluation with regards to complexity and error. Section 4 and Section 5 are where we formalize the experimental methodology and the empirical evaluation occurs respectively, focusing on shallow models, high-accuracy benchmarks, and energy efficiency. Finally, Section 6 concludes the paper with additional future work.

2. Proposed Method

The fully theory with all formal definitions and background required is presented here. In Appendix A we provide Table A1 which summarizes the notation used throughout the paper and in Appendix B, we present how the model parameters are learned.

2.1. Full Multivariate Chebyshev Expansion and Coefficient Tensor

Any function $f(x_1,\dots ,x_d)\in L_w^2\left({[-1,1]}^d\right)$ can be approximated by a truncated tensor-product Chebyshev expansion of maximum degree N:

$$f(x_1,\dots ,x_d)\approx {\displaystyle \sum _{i_1=0}^N}\dots {\displaystyle \sum _{i_d=0}^N}{a}_{i_1,\dots ,i_d}{\displaystyle \prod _{j=1}^d}{T}_{i_j}\left(x_j\right).$$

(1)

We define the coefficient tensor ${\mathcal{A}}_N$ of order d with mode size $(N+1)$ as

$${\mathcal{A}}_N[i_1,\dots ,i_d]=a_{i_1,\dots ,i_d},{\mathcal{A}}_N\in {\mathbb{R}}^{{(N+1)}^d}.$$

(2)

In subsequent subsections, we propose to approximate a high-dimensional function f by first computing its multivariate Chebyshev expansion and then applying a low-rank tensor decomposition.

2.2. Low-Rank Multivariate Approximation Using Chebyshev Tensor Decompositions

We propose to approximate the full coefficient tensor ${\mathcal{A}}_N$ from (2) using a rank-R CP decomposition:

$${\mathcal{A}}_{N,R}[i_1,\dots ,i_d]={\displaystyle \sum _{r=1}^R}{\lambda }_r{c}_{i_1,1,r}\dots c_{i_d,d,r}.$$

(3)

The resulting functional approximation of f is then given by:

$$f_{N,R}\left(\mathbf{x}\right)={\displaystyle \sum _{r=1}^R}{\lambda }_r{\displaystyle \prod _{j=1}^d}\left({\displaystyle \sum _{i=0}^N}{c}_{i,j,r}{T}_i(2x_j/R_{x_j}-1)\right),$$

(4)

where $T_i$ is the Chebyshev polynomial of degree i, $R_{x_j}$ is the range of variable $x_j$, N is the same truncation order for each variable, $c_{i,j,r}$ are coefficients, ${\lambda }_r$ is a scalar weight for component r, and R is the tensor rank. The total number of parameters reduces to $\mathcal{O}(R·d·N)$, achieving significant compression when $R\ll {(N+1)}^{d-1}$.

2.3. Tensor Formulation

We define the Chebyshev basis vector $\mathbf{T}\left(x_j\right)$ and their coefficient vector ${\mathbf{c}}_{j,r}$ as the following, where from this point onwards we implicitly assume all inputs have been normalized to $[-1,1]$, absorbing the affine scaling $2x_j/R_{x_j}-1$ from (4):

$$\mathbf{T}\left(x_j\right)={\left[T_0\left(x_j\right),T_1\left(x_j\right),\dots ,T_N\left(x_j\right)\right]}^{\top }\in {\mathbb{R}}^{N+1},$$

$$\mathrm{with}{\mathbf{c}}_{j,r}={\left[c_{0,j,r},c_{1,j,r},\dots ,c_{N,j,r}\right]}^{\top }.$$

Then each inner sum becomes a dot product ${\sum }_{i=0}^N{c}_{i,j,r}{T}_i\left(x_j\right)={\mathbf{c}}_{j,r}^{\top }\mathbf{T}\left(x_j\right).$ Thus (4) becomes

$$f_{N,R}\left(\mathbf{x}\right)={\displaystyle \sum _{r=1}^R}{\lambda }_r{\displaystyle \prod _{j=1}^d}\left({\mathbf{c}}_{j,r}^{\top }\mathbf{T}\left(x_j\right)\right)$$

(5)

From the above, the full tensor product Chebyshev feature is defined as $\mathbf{T}\left(\mathbf{x}\right)=\mathbf{T}\left(x_1\right)\otimes \dots \otimes \mathbf{T}\left(x_d\right)\in {\mathbb{R}}^{{(N+1)}^d},$ with their rank-1 coefficient tensor ${\mathbf{C}}_r={\lambda }_r{\mathbf{c}}_{1,r}\otimes {\mathbf{c}}_{2,r}\otimes \dots \otimes {\mathbf{c}}_{d,r}$, then the approximation becomes

$$f_{N,R}\left(\mathbf{x}\right)={\displaystyle \sum _{r=1}^R}\langle {\mathbf{C}}_r,\mathbf{T}\left(\mathbf{x}\right)\rangle ,$$

(6)

which is precisely the CANDECOMP/PARAFAC (CP) decomposition of a Chebyshev coefficient tensor.

2.4. Algorithm and Implementation Details

We describe the pseudocode in Algorithm 1 for the proposed method as follows:

Algorithm 1 Low-Rank Functional Approximation via CP Decomposition

1: Step 1: Compute Full Coefficients   2: Compute the full set of coefficients $a_{i_1,\dots ,i_d}$ using orthogonality: $$a_{i_1,\dots ,i_d}={\displaystyle \frac{2^d}{{\pi }^d}}{\int }_{-1}^1\dots {\int }_{-1}^{1}f\left(\mathbf{x}\right){\displaystyle \prod _{j=1}^d}{T}_{i_j}\left(x_j\right){\displaystyle \prod _{j=1}^d}w\left(x_j\right)d\mathbf{x}.$$   3: Step 2: Build Coefficient Tensor   4: Store the coefficients $a_{i_1,\dots ,i_d}$ into the tensor ${\mathcal{A}}_N\in {\mathbb{R}}^{{(N+1)}^d}$.   5: Step 3: CP Decomposition of Coefficient Tensor   6: Initialize factor vectors ${\mathbf{c}}_{j,r}$ randomly for $r=1,\dots ,R$ and $j=1,\dots ,d$.   7: Repeat until convergence (via Alternating Least Squares in Appendix B.1): $$\underset{{\mathbf{c}}_{j,r},{\lambda }_r}{min}{∥{\mathcal{A}}_N-{\displaystyle \sum _{r=1}^R}{\lambda }_r{\mathbf{c}}_{1,r}\otimes \dots \otimes {\mathbf{c}}_{d,r}∥}_F^2.$$   8: Normalize factors so that $\parallel {\mathbf{c}}_{j,r}{\parallel }_2=1$, absorbing the scaling into ${\lambda }_r$.   9: Step 4: Construct Functional Approximation 10: Assemble the univariate terms $g_{j,r}\left(x_j\right)$ from the decomposed factors: $${\displaystyle f_{N,R}(x_1,\dots ,x_d)={\displaystyle \sum _{r=1}^R}{\lambda }_r{\displaystyle \prod _{j=1}^d}{g}_{j,r}\left(x_j\right),}$$ where $${\displaystyle g_{j,r}\left(x_j\right)={\displaystyle \sum _{i=0}^N}{c}_{i,j,r}{T}_i\left(x_j\right).}$$

2.5. STEP As Neural Network Layer

We propose a new layer defined by:

$$y=\sigma \left({\displaystyle \sum _{r=1}^R}{\displaystyle \prod _{j=1}^d}{\mathbf{w}}_{j,r}^{\top }{\mathbf{t}}_j\right)$$

(7)

where ${\mathbf{t}}_j=\mathbf{T}\left(x_j\right)\in {\mathbb{R}}^{N+1}$ is the expanded Chebyshev polynomial feature vector for dimension j, and ${\mathbf{w}}_{j,r}\in {\mathbb{R}}^{N+1}$ are the learnable counterparts of the theoretical coefficient vectors ${\mathbf{c}}_{j,r}$ in the earlier subsections, parameterizing the mode-j activation weight for rank r.

2.5.1. Forward/Backward Propagation

Given an input vector $\mathbf{x}=(x_1,\dots ,x_d)$, the forward pass of the STEP layer proceeds as follows. First, each scalar input $x_j$ is expanded into its Chebyshev feature vector ${\mathbf{t}}_j$. For each rank component $r=1,\dots ,R$, we compute the separable product

$$z_r={\displaystyle \prod _{j=1}^d}{\mathbf{w}}_{j,r}^{\top }{\mathbf{t}}_j,$$

where ${\mathbf{w}}_{j,r}$ is the weight vector associated with input dimension j and rank r. The total pre-activation is obtained by summing over all rank components:

$$z={\displaystyle \sum _{r=1}^R}{z}_r.$$

Finally, the layer output is produced by applying a non-linear activation function $\sigma (·)$:

$$y=\sigma \left(z\right).$$

This formulation explicitly decouples the interactions across dimensions via the low-rank structure, yielding a parameter count that scales linearly with the input dimension d, the polynomial degree N, and the rank R, rather than exponentially.

We now discuss backward propagation. Let ${\sigma }^{\prime }\left(z\right)$ denote the derivative of the activation function. The gradient with respect to ${\mathbf{w}}_{j,r}$ is:

$${\displaystyle \frac{\partial y}{\partial {\mathbf{w}}_{j,r}}}={\sigma }^{\prime }\left(z\right)\left({\displaystyle \prod _{j=1,j^{\prime }\ne j}^d}{\mathbf{w}}_{j^{\prime },r}^{\top }{\mathbf{t}}_{j^{\prime }}\right){\mathbf{t}}_j.$$

2.6. STEP Embedding Layer

The STEP-Embedding mechanism regularizes the token representation space through a spectral bottleneck. The forward mapping is defined as:

$$\mathbf{z}=\mathrm{Mish}\left(\left(\mathrm{Emb}\left(x\right)\mathbf{B}\right){\alpha }_{\mathrm{in}}\right){\alpha }_{\mathrm{out}}+\mathbf{b},$$

(8)

where $\mathrm{Emb}\left(x\right)$[30] denotes the standard token embedding lookup and $\mathbf{B}\in {\mathbb{R}}^{R\times R}$ is the fixed Chebyshev basis matrix. By forcing the embedding through an R-dimensional manifold $\alpha$, the architecture filters out high-frequency noise inherent in high-dimensional dense embeddings.

3. Theoretical Performance and Complexity

In this section we analyze the proposed method’s properties of convergence, error, and complexity. We also apply it to specific models, and try to justify its efficiency using the proposed method.

3.1. STEP Framework

3.1.1. Convergence Analysis

In Appendix C we present a hierarchical analysis of the approximation error, decomposing it into contributions from polynomial truncation and tensor rank reduction. Here we state the main Theorem which ensures convergence of the model,

Theorem 1

(Two-Stage Convergence). Let $f\in L^2({[-1,1]}^d,w)$. For the approximation $f_{N,R}$ defined by the low-rank Chebyshev model,

$$\underset{N\to \infty }{lim}\underset{R\to \infty }{lim}{\parallel f-f_{N,R}\parallel }_{L^2}=0.$$

Furthermore, if f is analytic, the convergence with respect to N is exponential.

3.1.2. Error Analysis

We now analyze the total approximation error of the low-rank tensor Chebyshev expansion by decomposing it into two distinct components: the spectral truncation error and the low-rank CP approximation error.

Theorem 2

(Error Decomposition). Let $f\in L_w^2\left({[-1,1]}^d\right)$, let $f_N$ be its N-truncated tensor Chebyshev expansion, and let $f_{N,R}$ be the rank-R CP approximation. If f is analytic in a Bernstein ellipse of radius $\rho >1$, and bounded by the best rank-R approximation of the coefficient tensor, then the total $L^2$ error is bounded by $C{\rho }^{-N}+{\parallel {\mathcal{A}}_N-{\mathcal{A}}_{N,R}\parallel }_F$.

Proof. 

By the triangle inequality:

$$\begin{array}{cc}\hfill \parallel f-f_{N,R}{\parallel }_{L^2}& =\parallel (f-f_N)+(f_N-f_{N,R}){\parallel }_{L^2}\hfill \\ \hfill & \le \underset{\mathrm{Polynomial}\mathrm{truncation}\mathrm{error}}{\underbrace{\parallel f-f_N{\parallel }_{L^2}}}+\underset{\mathrm{Tensor}\mathrm{rank}\mathrm{error}}{\underbrace{\parallel f_N-f_{N,R}{\parallel }_{L^2}}}\hfill \\ \hfill & \le C{\rho }^{-N}+{\parallel {\mathcal{A}}_N-{\mathcal{A}}_{N,R}\parallel }_F.\hfill \end{array}$$

□

Remark: If f has bounded mixed partial derivatives up to order s, or is analytic in a Bernstein ellipse of radius $\rho >1$, then we have:

$$\parallel f-f_N{\parallel }_{L^2}\le C{N}^{-s}\mathrm{or}\mathrm{equivalently}{\parallel f-f_N\parallel }_{\infty }\le C{\rho }^{-N}.$$

Furthermore, if the tensor rank error is bounded by the best rank-R approximation error of the coefficient tensor, then we have:

$$\parallel f_N-f_{N,R}{\parallel }_{L^2}={\parallel {\mathcal{A}}_N-{\mathcal{A}}_{N,R}\parallel }_F.$$

3.1.3. Complexity Analysis

Table 2 compares the computational complexity for function evaluation between the classical full tensor method and the proposed low-rank approach. The classical method exhibits exponential complexity of $\mathcal{O}\left({(N+1)}^d\right)$, whereas the proposed method reduces this to $\mathcal{O}(R·d·N)$, which is linear in both dimension d and polynomial degree N. The complexity ratio,

$${\displaystyle \frac{\mathrm{Classical}\mathrm{Complexity}}{\mathrm{Proposed}\mathrm{Complexity}}}={\displaystyle \frac{{(N+1)}^d}{R·d·N}}=\mathrm{Exponential}\mathrm{gain},$$

demonstrates an exponential gain in favor of the low-rank approximation, particularly for high-dimensional problems. This reduction in complexity makes the proposed approach well-suited for applications where d is large and fast function evaluations are required.

3.2. STEP Neural Network Layer

It is noted that the STEP layer from (7) has a parameter count:

$$R{\displaystyle \sum _{j=1}^d}\dim \left({\mathbf{t}}_j\right)=R·d·(N+1),$$

where $\dim \left({\mathbf{t}}_j\right)=N+1$ denotes the length of the expanded feature vector. The layer performs a mapping $\mathbf{y}=\sigma \left(z\right)$, where the pre-activation $z\in {\mathbb{R}}^U$ is computed through two distinct stages:

1.

Spectral Expansion: Given input $\mathbf{x}\in {\mathbb{R}}^d$, we first apply Min-Max normalization to map $\mathbf{x}$ to the domain $[-1,1]$, followed by the Chebyshev lifting:

$${\mathbf{x}}_{\mathrm{spec}}=\mathrm{Norm}\left(\mathbf{x}\right)\mathbf{V},\mathrm{where}\mathbf{V}\in {\mathbb{R}}^{d\times d}.$$

(9)

2.

Low-Rank Projection: The spectral features ${\mathbf{x}}_{\mathrm{spec}}$ are compressed through a rank-R bottleneck defined by learnable factors ${\alpha }_{\mathrm{in}}\in {\mathbb{R}}^{d\times R}$ and ${\alpha }_{\mathrm{out}}\in {\mathbb{R}}^{R\times U}$:

$$z=\left({\mathbf{x}}_{\mathrm{spec}}{\alpha }_{\mathrm{in}}\right){\alpha }_{\mathrm{out}}+\mathbf{b}.$$

(10)

This composition reduces the parameter count from $\mathcal{O}(d·U)$ in a standard dense layer to $\mathcal{O}\left(R\right(d+U\left)\right)$.

4. Experimental Methodology

We conduct a comprehensive study spanning classical high-dimensional function approximation and modern convolutional neural network architectures. The experiments are designed to highlight STEP’s ability to drastically reduce parameter count and computational cost while maintaining or improving accuracy relative to full-rank, conventional, and deep benchmark models.

For each benchmark function or network architecture, we perform the following systematic analysis:

1.

Function or Model Definition: We formally define the target function $f\left(\mathbf{x}\right):{[-1,1]}^d\to \mathbb{R}$ or the CNN architecture under study. For functions, we outline their domain and classical full-tensor polynomial expansion. For CNNs, we specify the baseline kernel size, layer configuration, and dataset such as MNIST [31].

2.

STEP Construction and Low-Rank Analysis: We construct the low-rank tensorized approximation $f_{N,R}$ or decomposed convolutional filters using CP decomposition or mode-wise separable factorization. The decay of CP weights (${\lambda }_r$) or singular modes is analyzed to justify the selected rank R.

3.

Comparative Performance Metrics:

• Accuracy Comparison: Quantitative evaluation of approximation errors for functions ($\parallel f-f_{\mathrm{approx}}{\parallel }_{L^2}$) and classification accuracy for CNNs. Comparisons include full-rank expansions, sparse grids, Tensor-Train decompositions, and baseline deep models.
• Efficiency Comparison: Metrics include parameter count, storage requirements, and computational cost (operations per evaluation, multiplications/additions per pixel or voxel). Tables summarize gains achieved through low-rank tensorization.

4.

Summary of Findings: For each experiment, we highlight the minimum rank R required to achieve target accuracy, quantify compression and speedup factors, and discuss implications for shallow STEP networks, Green AI deployment, and high-dimensional learning efficiency.

This structured approach ensures that all experiments, from smooth function approximation to CNN compression and benchmark deep networks, are analyzed consistently. The results empirically confirm the theoretical advantages of STEP in reducing computational complexity while preserving or enhancing predictive accuracy.

5. Results

This section provides a rigorous empirical validation of the Separable Tensor Polynomial (STEP) framework. Additional results, namely the evolution of validation accuracy over training epochs together with model size and CPU inference latency for each experiment, are available in Appendix D, which highlights the stability of STEP over training.

5.1. STEP For Function Approximation

In this section, we provide a unified presentation of illustrative examples demonstrating the power of tensorized polynomial approximations using the STEP framework.

We use the following functions below and try to approximate them with STEP.

1.

Function 1 (2D):$f(x_1,x_2)=exp\left(x_1{x}_2\right)$, rank $R=2$.

2.

Function 2 (3D):$f(x_1,x_2,x_3)={\displaystyle \frac{1}{1+10(x_1^2+x_2^2+x_3^2)}}$, smooth radial function; classical ${(N+1)}^3$ vs STEP $3RN$.

3.

Function 3 (4D):$f(x_1,\dots ,x_4)=exp(x_1{x}_2+x_3{x}_4)$, rank $R=2$.

For each function, we include: definition, classical approximation, tensorized approximation construction, error comparison tables, and discussion.

From the above Table 3, we can see that across all examples, the STEP framework consistently demonstrates parameter efficiency as it dramatically reduces parameters compared to the full tensor approximations. Along with this is the increased precision as for both low and moderate dimensional functions, small errors can be achieved even if the rank R is small. We also note the advantages of scalability as higher dimensions require parameters without exponential growth and the ability to control accuracy through increasing the rank R which implies that the framework works for deep learning models.

5.2. STEP vs Traditional MLP

This section evaluates STEP model using the formulation discussed in Section 3.2 on two different datasets: a synthetic function dataset, and the Fashion-MNIST [32] dataset converted into a 1D representation. These datasets represent different levels of complexity, including nonlinear functional learning, structured mathematical patterns, and real-world image classification.

The workflow is shown in Figure 2 for clarity of the experiment. For a fair comparison, the traditional MLP and the STEP-based MLP are trained under identical experimental conditions in terms of epochs and dataset splits. While the traditional MLP relies on fully connected dense layers, the STEP model replaces the integration stage with a spectral tensorized projection that significantly reduces the number of parameters while maintaining competitive accuracy. We observe throughout the experiments here that (1) we can model high dimensional nonlinear interactions with a parameter count that scales linearly with dimension d and rank R, avoiding the exponential blow-up typical of full tensors, (2) gradients are explicitly defined, ensuring STEP integrates seamlessly with standard optimization routines, even in large-scale high-dimensional settings allowing for drop in replacements to standard models, (3) STEP embodies Green AI principles by substantially reducing memory and energy demands while maintaining strong predictive accuracy, making it highly suitable for resource-constrained and edge deployment scenarios.

5.2.1. Experiment 1: Function Dataset

The first experiment evaluates both models on a nonlinear synthetic function dataset defined as:

$$y=\left(sin\left(X_0\right)cos\left(X_1\right)+{\displaystyle \frac{X_2^2}{100}}+(X_4>2)\right)>0$$

Purely from a function approximation perspective, we see that STEP better captures its shape due to higher validation accuracy and doesn’t fit to noise as much compared to the MLP in Table 4. This also hints that deep learning models easily overfit in terms of parameter size and STEP provides a concrete answer to parameter efficiency without either overfitting or underfitting.

5.2.2. Experiment 2: Fashion-MNIST (Converted to 1D)

In the second experiment, the models are evaluated on the Fashion-MNIST dataset after converting the images into a 1D representation.

Here, STEP demonstrates again its parameter efficiency in Table 5. It uses 5% of parameters of the MLP while having comparable validation accuracy. Due to the reduction in parameters, inference speed has increased slightly, in further experiments in the next section, we demonstrate the speedup is significant for larger models in the parameter sizes up to millions.

5.3. STEP vs. Traditional CNN

This section presents a controlled comparison between the (STEP) model and a traditional Convolutional Neural Network (CNN) using the MedMNIST [33], PneumoniaMNIST, and Imagenette [34] datasets. Models are trained under identical experimental conditions to ensure a fair and unbiased comparison in terms of accuracy, model complexity, and computational efficiency.

The main difference between the two models lies in the integration mechanism. While the traditional CNN relies on spatial flattening followed by dense fully connected layers, the STEP model replaces this stage with localized spectral projection using Chebyshev polynomial filters defined on the interval $[-1,1]$ as shown in Figure 3.

5.3.1. Experiment 1: MedMNIST Dataset (Epoch 15)

The first experiment evaluates the performance of the STEP model compared with a traditional CNN on the MedMNIST dataset after 15 training epochs.

Running on the MedMNIST dataset, we see in Table 6 the potential for real-time applications with large models as we can have comparable accuracy while only using 10% of the parameters. In addition to this, latency is halved, further supporting the use of STEP in embedded devices due to reduced memory and latency.

5.3.2. Experiment 2: PneumoniaMNIST Dataset (Epoch 15)

To further evaluate the robustness of the STEP architecture, we compare it with traditional deep learning models on the PneumoniaMNIST dataset.

From Table 7, the STEP model is smaller by a wide margin than every baseline against which it is compared. This is especially prevalent in the case where we compare STEP against ResNet50, where STEP is over 550× smaller while achieving higher accuracy (0.8333 vs. 0.8077); it also slightly outperforms the Traditional CNN at less than a ninth of its size. The only case where accuracy is significantly lower is for the VGG16 model, but it should be noted that STEP is 346.7x smaller so there has to be a tradeoff between efficiency and performance to some degree.

5.3.3. Experiment 3: Imagenette Dataset (Epoch 10)

To evaluate scalability, we extend the comparison to the Imagenette dataset using models trained from scratch.

From the Table 8 above, STEP’s shallow tensorized architecture attains competitive accuracy relative to deeper and heavily parameterized models like VGG16 and ResNet50. Remarkably, STEP uses orders of magnitude fewer parameters (126k vs 24.1M for ResNet50) while maintaining superior validation accuracy (0.586 vs 0.308 at Epoch 10), demonstrating that depth and parameter count are not strictly necessary to achieve high performance when expressive spectral representations are employed.

5.3.4. Experiments 4 & 5: Hybrid VGG16-STEP and ResNet50-STEP on the Imagenette Dataset

This section evaluates the STEP layer integrated into VGG16 and ResNet50 architectures on the Imagenette dataset. The goal is to assess how spectral tensor expansion reduces model complexity while maintaining classification performance across different deep architectures.

Similar to the previous subsection, incorporating STEP layers maintains accuracy as shown in Table 9 as shown in the comparison between the Hybrid and Standard VGG16 while reducing the parameter size up to half. The comparison between Hybrid and Standard ResNet50 is in contrast as it shows extreme compression with only roughly 10% of the original parameters while performing competitively, it is noted that the latency is significantly reduced due to the parameter reduction.

5.4. STEP VS LLM Tokenizer Generations

To evaluate the scalability and robustness of the STEP architecture, we benchmark the STEP-Embedding layer against standard dense embedding matrices across three sizes of Large Language Model tokenizers [35]: GPT-2 [36] size tokenizer (vocab size: 50,257), GPT-3 [37] p50k_base size tokenizer (vocab size: 50,281), and GPT-4 cl100k_base size tokenizer (vocab size: 100,277). The target embedding dimension was held constant at 256. We analyzed both parameter efficiency and geometric fidelity. This projection from (8) acts as a structural regularizer, preserving the semantic identity of tokens while achieving $>96\%$ parameter reduction. As shown in the workflow of Figure 4, we set a low rank bottleneck $R=8$ for compression. Appendix E provides a graphical understanding of how the embedding space is affected.

The traditional embedding approach scales linearly with vocabulary size, resulting in severe parameter bloat. For the GPT-4 size tokenizer, the standard dense matrix requires over 25.6 million parameters just for the input layer. By routing the token identities through an 8-dimensional Chebyshev manifold, the STEP architecture achieved a massive reduction in memory footprint across all generations, consistently exceeding $96.8\%$ compression (see Table 10). A primary concern with extreme bottleneck compression is the loss of token identity (linear collapse). To test this, we performed a reverse-decoding proof using Nearest-Neighbor Euclidean Distance ($L_2$ Norm) to map the outputted 256-dimensional continuous vectors back to the discrete vocabulary space. Across diverse grammatical structures and all model generations, the STEP Rank-8 layer achieved a $100\%$ string-match reconstruction rate as shown in Table 10. This empirically shows that the inner-manifold Mish activation [38] prevents structural amnesia and preserves absolute semantic identity.

6. Conclusions

This work introduced the Separable Tensor Polynomial (STEP) framework, a principled approach for low-rank tensorized polynomial expansions in both function approximation and deep learning. STEP overcomes the classical curse of dimensionality by employing structured tensor decompositions, reducing parameter growth from exponential to linear in the dimension ($\mathcal{O}(R·d·N)$) while maintaining high approximation fidelity. Further experiments on the Imagenette dataset confirmed that shallow STEP models achieve competitive validation accuracy relative to deeper architectures such as VGG16 and ResNet50, while using orders of magnitude fewer parameters (126k vs. 24.1M for ResNet50) and lower inference latency. Integrated directly into these architectures, STEP reduces parameter counts by up to 91.8% while preserving over 97% of classification accuracy. Notably, a single STEP layer consistently replaces multiple conventional dense or convolutional layers across experiments, with very few trainable parameters relative to traditional baselines, confirming that a shallow architecture need not sacrifice competitive performance. These results demonstrate that STEP provides a powerful mechanism for model compression, accelerated computation, and energy-efficient deployment, making it highly suitable for resource-constrained and Green AI scenarios. Future work will focus on extending STEP to larger datasets, exploring adaptive rank selection, integrating with pruning and quantization [39], and assessing its deployment efficiency in edge AI environments.

Appendix B.1. Alternating Least Squares

We aim to minimize the mean-squared-error loss between the true values $f^{\left(k\right)}$ and their approximations $f_{N,R}^{\left(k\right)}$:

$$L={\displaystyle \sum _{k=1}^M}{\left(f^{\left(k\right)}-f_{N,R}^{\left(k\right)}\right)}^2.$$

The optimization is performed via an alternating least-squares (ALS) scheme. At each step, we update the coefficient vector ${\mathbf{c}}_{j,r}$ for a single variable j and rank index r while keeping all other coefficients fixed. To isolate its contribution, we define the intermediate quantity

$$Q_{j,r}^{\left(k\right)}={\lambda }_r{\displaystyle \prod _{\begin{array}{c}{j}^{\prime }=1,\\ j^{\prime }\ne j\end{array}}^d}\left({\mathbf{c}}_{j^{\prime },r}^{\top }\mathbf{T}\left(x_{j^{\prime }}^{\left(k\right)}\right)\right),$$

which collects the product of all other terms in the rank-r component. The model approximation can then be rewritten linearly in the target coefficients as

$$f_{N,R}^{\left(k\right)}={\displaystyle \sum _{r=1}^R}{Q}_{j,r}^{\left(k\right)}\left({\mathbf{c}}_{j,r}^{\top }\mathbf{T}\left(x_j^{\left(k\right)}\right)\right).$$

Let $\mathbf{y}={(f^{\left(1\right)},\dots ,f^{\left(M\right)})}^{\top }$ be the vector of observed outputs and define the diagonal weight matrix $W_r=diag(Q_{j,r}^{\left(1\right)},\dots ,Q_{j,r}^{\left(M\right)})$ for each r. Substituting the linearized form into the loss L yields a weighted least-squares problem for ${\mathbf{c}}_{j,r}$, with the solution

$${\mathbf{c}}_{j,r}={\left({\Phi }_j^{\top }{W}_r^{\top }{W}_r{\Phi }_j\right)}^{-1}{\Phi }_j^{\top }{W}_r^{\top }{W}_r\mathbf{y},$$

(A1)

where ${\Phi }_j$ is the Chebyshev design matrix for variable j defined previously. Iterating this update cyclically over all j and r constitutes a standard ALS algorithm for fitting the model.

Appendix B.2. Gradient-Based Optimization

The model is differentiable, and

$$\begin{array}{c}{\displaystyle \frac{\partial L}{\partial c_{i,j,r}}}=-2{\displaystyle \sum _{k=1}^M}\left(f^{\left(k\right)}-f_{N,R}^{\left(k\right)}\right)\hfill \\ \hfill \times {\lambda }_r{T}_i\left(x_j^{\left(k\right)}\right)\prod _{j^{\prime }\ne j}\left({\mathbf{c}}_{j^{\prime },r}^{\top }\mathbf{T}\left(x_{j^{\prime }}^{\left(k\right)}\right)\right),\end{array}$$

applying the chain rule. This allows optimization via Adam, RMSprop, or L-BFGS.

Appendix C.1. Orthogonal Basis and Series Expansion

The Chebyshev tensor basis functions, defined for a multi-index $\mathbf{i}=(i_1,\dots ,i_d)\in {\mathbb{N}}^d$ as

$$T_{\mathbf{i}}\left(\mathbf{x}\right)={\displaystyle \prod _{j=1}^d}{T}_{i_j}\left(x_j\right),$$

form a complete orthogonal system in the weighted space $L^2({[-1,1]}^d,w)$ with the product weight $w\left(\mathbf{x}\right)={\prod }_{j=1}^d{(1-x_j^2)}^{-1/2}$. Consequently, any function $f\in L^2({[-1,1]}^d,w)$ admits a unique expansion

$$f\left(\mathbf{x}\right)=\sum _{\mathbf{i}}{a}_{\mathbf{i}}{T}_{\mathbf{i}}\left(\mathbf{x}\right),$$

which converges in the $L^2$ norm.

Appendix C.2. Spectral Convergence of Polynomial Truncation

Truncating the series to maximum degree N in each variable yields the approximation

$$f_N\left(\mathbf{x}\right)={\displaystyle \sum _{i_1=0}^N}\dots {\displaystyle \sum _{i_d=0}^N}{a}_{i_1,\dots ,i_d}{T}_{\mathbf{i}}\left(\mathbf{x}\right).$$

This truncation is the orthogonal projection onto the subspace of polynomials of degree at most N, guaranteeing that $\parallel f-f_N{\parallel }_{L^2}\to 0$ as $N\to \infty$. The convergence rate depends on the smoothness of f:

• If f has s continuous derivatives, the error decays algebraically: $\parallel f-f_N{\parallel }_{L^2}=\mathcal{O}\left(N^{-s}\right)$.
• If f is analytic on a Bernstein polyellipse, the error decays exponentially: there exists $\rho >1$ such that $\parallel f-f_N{\parallel }_{L^2}=\mathcal{O}\left({\rho }^{-N}\right)$.

Appendix C.3. Low-Rank Tensor Approximation

Let ${\mathcal{A}}_N\in {\mathbb{R}}^{{(N+1)}^d}$ denote the coefficient tensor of $f_N$. To circumvent the exponential growth of terms, we employ a rank-R CANDECOMP/PARAFAC (CP) approximation ${\mathcal{A}}_{N,R}$. A fundamental property of the CP format is that any finite tensor can be approximated arbitrarily well with sufficient rank, i.e.,

$$\parallel {\mathcal{A}}_N-{\mathcal{A}}_{N,R}{\parallel }_F\to 0\mathrm{as}R\to \infty .$$

Since the $L^2$ norm of the corresponding function approximation is controlled by the Frobenius norm of the coefficient tensor, we have

$$\parallel f_N-f_{N,R}{\parallel }_{L^2}\le C{\parallel {\mathcal{A}}_N-{\mathcal{A}}_{N,R}\parallel }_F\to 0.$$

Appendix C.4. Total Error Decomposition and Main Result

The total approximation error decomposes via the triangle inequality:

$$\parallel f-f_{N,R}{\parallel }_{L^2}\le \underset{\mathrm{truncation}\mathrm{error}}{\underbrace{\parallel f-f_N{\parallel }_{L^2}}}+\underset{\mathrm{low}-\mathrm{rank}\mathrm{error}}{\underbrace{\parallel f_N-f_{N,R}{\parallel }_{L^2}}}.$$

The first term vanishes as the polynomial degree N increases, and the second term vanishes as the tensor rank R increases. This leads to our main convergence guarantee in Theorem 1.