Paper deep dive
Spectral Higher-Order Neural Networks
Gianluca Peri, Timoteo Carletti, Duccio Fanelli, Diego Febbe
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 94%
Last extracted: 3/31/2026, 2:44:49 AM
Summary
The paper introduces Spectral Higher-Order Neural Networks (SHONNs), a novel architecture that incorporates higher-order (triadic) interactions into feedforward neural networks. By utilizing a spectral reparameterization of weight matrices, SHONNs reduce the parameter scaling complexity from O(N^3) to O(N^2), matching standard Multi-Layer Perceptrons (MLPs) while maintaining superior expressivity and stability. The authors prove universal approximation capabilities and demonstrate empirical performance gains on MNIST, Fashion-MNIST, and CIFAR-10 datasets, including integration into MLP-Mixer architectures.
Entities (5)
Relation Signals (3)
SHONNs → reducescomplexityto → O(N^2)
confidence 95% · The computational cost (as measured by the parameters scaling) is reduced to O(N^2) due to an effective weight reparametrization
SHONNs → generalizes → MLP
confidence 90% · The objective of our work is to address the intrinsic limitations that have hindered the broader adoption and advancement of higher-order neural networks... through Spectral Higher-Order Neural Networks (shonns), where triadic interactions are integrated into the standard mlp forward pass.
SHONNs → outperforms → MLP
confidence 90% · The results show a clear advantage for the triadic architecture.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Neural networks are fundamental tools of modern machine learning. The standard paradigm assumes binary interactions (across feedforward linear passes) between inter-tangled units, organized in sequential layers. Generalized architectures have been also designed that move beyond pairwise interactions, so as to account for higher-order couplings among computing neurons. Higher-order networks are however usually deployed as augmented graph neural networks (GNNs), and, as such, prove solely advantageous in contexts where the input exhibits an explicit hypergraph structure. Here, we present Spectral Higher-Order Neural Networks (SHONNs), a new algorithmic strategy to incorporate higher-order interactions in general-purpose, feedforward, network structures. SHONNs leverages a reformulation of the model in terms of spectral attributes. This allows to mitigate the common stability and parameter scaling problems that come along weighted, higher-order, forward propagations.
Tags
Links
- Source: https://arxiv.org/abs/2603.28420v1
- Canonical: https://arxiv.org/abs/2603.28420v1
Trouble viewing inline? Open PDF directly →
Full Text
66,216 characters extracted from source content.
Expand or collapse full text
Spectral Higher-Order Neural Networks Gianluca Peri1, Timoteo Carletti2, Duccio Fanelli1, Diego Febbe1 1University of Florence, Italy 2University of Namur, Belgium gianluca.peri, duccio.fanelli, diego.febbe@unifi.it timoteo.carletti@unamur.be Abstract Neural networks are fundamental tools of modern machine learning. The standard paradigm assumes binary interactions (across feedforward linear passes) between inter-tangled units, organized in sequential layers. Generalized architectures have been also designed that move beyond pairwise interactions, so as to account for higher-order couplings among computing neurons. Higher-order networks are however usually deployed as augmented graph neural networks (gnns), and, as such, prove solely advantageous in contexts where the input exhibits an explicit hypergraph structure. Here, we present Spectral Higher-Order Neural Networks (shonns), a new algorithmic strategy to incorporate higher-order interactions in general-purpose, feedforward, network structures. shonns leverages a reformulation of the model in terms of spectral attributes. This allows to mitigate the common stability and parameter scaling problems that come along weighted, higher-order, forward propagations. 1 Introduction Higher-Order Neural Networks (honns) can claim a long history in the field of machine learning, despite having being overlooked for many years. To the best of our knowledge, the first conceptual attempt to account for higher order schemes dates back to 1986, with the introduction of the so called sigma-pi units [22]. These are structures that first group the inputs into clusters, then compute the product within each cluster, and finally sum the obtained products. Sigma-pi units have been conceived, from the very beginning, as neural networks. The acronym honn became popular only one year after their actual formalization [6, 13]. The field of modern deep learning has however steered to progressively favor dense neural network models: pairwise interactions among neurons are combined into weighted sums before being fed, as an input, to apposite non linear functions. This orientation has persisted despite the fact that architectures featuring higher-order interactions potentially offer greater expressivity than those restricted to binary couplings. This is likely attributable to the ensuing parameters scaling and associated computational costs. For a typical feedforward neural network model, also termed Multi Layer Perceptron (mlp), one can expect a parameter scaling of the order O(N2)O(N^2), for any given pair of nested layers of size N. Conversely, for a higher-order network with triadic interactions (include two body correlations from the departing layer) the total parameter count scales as O(N3)O(N^3). The latter approach is impractical due to the prohibitively high training times required for real-world applications. Furthermore, the necessity of navigating an extensive parameter space often results in suboptimal solutions. Despite various efforts to mitigate the parameter scaling problem, honns have generally failed to achieve the same level of mainstream adoption as multi-layer perceptron mlp models. Two notable examples include the aforementioned pi-sigma units [22], and Π -nets [5, 4]. In recent years, several machine learning subfields have emerged from the integration of various higher-order processing techniques. Prominent among these is the attention mechanism employed in modern transformer architectures; other notable examples include Gated Linear Units (glu) and factorization machines for recommendation systems [23, 21, 19]. The field’s most direct adoption of higher-order neural networks arguably occurs within Hypergraph Neural Networks (hgnns) and Simplicial Neural Networks (snns) [11, 7]. Both architectures generalize Graph Neural Networks (gnns) to inputs with higher-order structures by adapting the message-passing algorithms to hypergraphs. It should be noted, however, that these models not only suffer from the curse of dimensionality (O(N3)O(N^3) parameter scaling, a lower bound which applies to the most conservative setting where just triadic interactions are incorporated) but are also inherently constrained to datasets possessing an explicit hypergraph structure. The objective of our work is to address the intrinsic limitations that have hindered the broader adoption and advancement of higher-order neural networks. This is achieved through Spectral Higher-Order Neural Networks (shonns), where triadic interactions are integrated into the standard mlp forward pass. The computational cost (as measured by the parameters scaling) is reduced to O(N2)O(N^2) due to an effective weight reparametrization that exploits the spectral attributes of the involved transfer matrices [12]. Moreover, we demonstrate that shonns possess full expressivity, as they are capable of universal approximation for any continuous function. In a manner similar to hgnns and snns, our proposed scheme operates on complete hypergraphs of neurons; consequently, higher-order interactions are distributed throughout the architecture rather than being confined to isolated functional blocks. In the following, we present the mathematical foundation of honns, focusing on the specific case of triadic interactions. These serve as the fundamental building blocks of the higher-order neural structure. We will then turn to demonstrate their universal approximation properties, as anticipated above. Dedicated benchmark tests have been performed, both in classification and regression modalities. Our results suggest that higher-order networks, integrated with the spectral parametrization, achieve superior performance over traditional low-order counterparts, without incurring additional computational overhead. 2 Results 2.1 Higher-Order Layers x1x_1x2x_2x3x_3x4x_4x5x_5y1y_1y2y_2y3y_3yk=∑i=0Nwkixi y_k= _i=0^Nw_kix_i (a) Standard connectivity of a neural network. x1x_1x2x_2x3x_3x4x_4x5x_5y1y_1y2y_2y3y_3yk=∑i=0Nwkixi+∑0≤i≤jN−1w~kijxixj aligned &y_k= _i=0^Nw_kix_i+\\ & _0≤ i≤ j^N-1 w_kijx_ix_j aligned (b) Standard higher-order (triadic) forward propagation structure. x1x_1x2x_2x3x_3x4x_4x5x_5y1y_1y2y_2y3y_3yk=∑i=0Nwkixi+∑i≤j(λ~ij(in)−λ~k(out))ϕ~kijxixj aligned &y_k= _i=0^Nw_kix_i+\\ & _i≤ j( λ_ij^(in)- λ_k^(out)) φ_kijx_ix_j aligned (c) Forward propagation of a spectral higher-order neural network. Figure 1: Cartoon representing different neural architectures. The standard neural networks (panel 1(a)) with input neurons, xix_i, generating output signals, yky_k, via weighted averages, i.e., linear combinations (colored arrows) followed by the application of a local nonlinearity (not shown). The standard triadic higher-order network (panel 1(b)) is obtained by adding to the previous architecture, weighted sums of hyperlinks, mimicking the two body interaction xi1xi2x_i_1x_i_2 (symbolized by the curved arrows connecting a couple (xi1,xi2)(x_i_1,x_i_2) to an output yky_k, each pair with its own specific color). The spectral higher-order networks (panel 1(c)), reduces the number of used parameters by exploiting the spectral decomposition. This yields an effective parameter sharing among hyperlinks (the curved arrows connecting a couple (xi1,xi2)(x_i_1,x_i_2) to several outputs yky_k, share the same color). Consider a standard mlp and focus on two adjacent layers. Assume xix_i to denote the signal associated to the departing layer, made of N individual units. The signal yky_k referred to the k-th neuron at the arriving layer can be computed as follows the usual forward propagation, namely: yk=∑i=0Nwkixi,y_k= _i=0^Nw_kix_i, (1) where wkiw_ki is the k,ik,i-th element of the rectangular weight matrix WK×N+1(1)W_K× N+1^(1) and K stands for the size of the second layer (see panel (a) Fig. 1). Notice that the sum in (1) runs on N+1N+1 elements, because we have chosen to deploy the scalar bias via a dedicated bias neuron, xNx_N. An element-wise non-linearity σ:ℝK→ℝKσ:R^K ^K is then applied to the obtained activation vector =(y1,…,yK)⊤ y=(y_1,…,y_K) , to expand the model’s hypothesis space. We can extend the above scheme beyond 11-simplexes by accounting for higher-order interactions between layers. We begin by adding a quadratic term and thus recast (1) in the form: yk=∑i=0Nwkixi+∑0≤i≤jN−1w~kijxixj,y_k= _i=0^Nw_kix_i+ _0≤ i≤ j^N-1 w_kijx_ix_j\,, (2) where w~kij w_kij is now the element of a tridimensional tensor K×N×NK× N× N. See middle panel of Fig. 1 for a graphical representation of the ensuing pattern of interaction. Without loss of generality, we have chosen not to include a bias term in the higher-order interaction defined in (2). Henceforth, we shall assume this distinction without further mention. Notably, the non-linearity is integrated directly into the network topology, thereby obviating the requirement for an explicit non linear post-processing activation step - a sensible departure from standard architectures. Model (2) can be seen as a second order approximation of the full honn expansion [13] yk=wk(0)+∑iwki(1)xi+∑i,jwkij(2)xixj+∑i,j,ℓwkijℓ(3)xixjxℓ+….y_k=w_k^(0)+ _iw_ki^(1)x_i+ _i,jw_kij^(2)x_ix_j+ _i,j, w_kij ^(3)x_ix_jx_ +…\,. (3) However, it should be noted that in practical implementations, the aforementioned sum must be truncated; under these conditions, the expressive power of (3) remains to be proven. Additionally, the rapid growth in the parameters’ count renders any practical implementation computationally prohibitive. From a more general perspective, model (3) can be interpreted as the Taylor expansion 111By setting n=2n=2 one can approximate f(xi,xj;)f(x_i,x_j; θ) as follows f(xi,xj;)=f(;)+∂f∂xixi+∂f∂xjxj+12[∂2f∂xi2xi2+2∂2f∂xi∂xjxixj+∂2f∂xj2xj2]+….f(x_i,x_j; θ)=f( 0; θ)+ ∂ f∂ x_ix_i+ ∂ f∂ x_jx_j+ 12 [ ∂^2f∂ x_i^2x_i^2+2 ∂^2f∂ x_i∂ x_jx_ix_j+ ∂^2f∂ x_j^2x_j^2 ]+...\,. (4) The 0-th and 11-th orders can be identified with the bias and the linear term, i.e., a standard mlp. The truly novel term, is the second order one, xixjx_ix_j, i,j∈1,2i,j∈\1,2\. of the following non linear update rule yk=∑i=0Nwkixi+∑0≤i≤j≤…≤nN−1f(xi,xj,…,xn;(k,i,j,…,n)),y_k= _i=0^Nw_kix_i+ _0≤ i≤ j≤...≤ n^N-1f (x_i,x_j,...,x_n; θ(k,i,j,...,n) )\,, (5) where f is a generic coupling function and (k,i,j,…,n) θ(k,i,j,...,n) sets the pattern of active higher order correlations. As we shall prove, model choice (2) possesses universal approximation capabilities, provided the transformation is iteratively applied through a deep sequence of hidden layers. Stated differently, one can approximate any continuous function by coupling sufficiently many quadratic interactions of the type accommodated for in (2) (see Appendix A). This implies a drastic decrease in complexity compared to (3), thus enabling more scalable implementations. A primary factor hindering the widespread adoption of the latter models within the deep learning community has to do with a phenomenon closely related to the curse of dimensionality. The number of parameters needed for an usual mlp network of fixed depth grows as N2N^2, where N refers to size of the layers. The triadic hyperedge model suffers from an O(N3)O(N^3) parameter explosion, which substantially increases the computational burden in practical applications. The main result of this work is a solution to this latter problem, via an efficient reparameterization scheme that induces weight sharing across the architecture. The recipe represents a generalization of the spectral parametrization which has been recently proposed as a viable alternative to the usual optimization frameworks [3] (see section 4.1). Applied to the classical mlp (1), the spectral formulation yields yk=∑i(λi(in)−λk(out))ϕkixi,y_k= _i( _i^(in)- _k^(out)) _kix_i\,, (6) where elements ϕki _ki define the eigenvectors of a square transfer operator, (in)∈ℝN λ^(in) ^N and (out)∈ℝK λ^(out) ^K the associated eigenvalues, referred to departure and arrival nodes, respectively. By optimizing only the eigenvalues and preserving the original eigenvectors the cost is O(N)O(N), we achieve thus a significant efficiency gain over the original setting, which requires direct weight optimization. By extending the spectral paradigm to accommodate higher-order forward propagation through triadic interactions, as specified in (2), one eventually obtains: yk=∑i(λi(in)−λk(out))ϕkixi+∑i≤j(λ~ij(in)−λ~k(out))ϕ~kijxixj.y_k= _i( _i^(in)- _k^(out)) _kix_i+ _i≤ j( λ_ij^(in)- λ_k^(out)) φ_kijx_ix_j\,. (7) The derivation of the above result is given in the Methods section. A pictorial representation of the spectral higher order algorithm is provided in the right panel of Fig. 1. Hold constant the scalar quantities ϕ~kij φ_kij defined within the triadic summation, and just train the generalized eigenvalues λ~ij(in) λ_ij^(in), λ~k(out) λ_k^(out) together with the entire parameter suite associated with the linear forward pass (both the eigenvectors and the eigenvalues). Consequently, the parametric complexity required to optimize the generalized model—including two-body interactions—is reduced to O(N2)O(N^2), in contrast to the O(N3)O(N^3) scaling inherent in a naive implementation. We refer to Appendix B for a comprehensive analysis of the involved parameter scalings. In summary, by adopting the spectral ansatz, the triadic higher-order model can be trained with the same parametric complexity as a standard mlp 222Recall that, under the family of gd optimizers routinely used for deep learning application, every omitted parameter represents a corresponding reduction in the total count of partial derivatives required during backpropagation.. In the following section, we will begin to explore the practical side of dealing with a spectral triadic architecture. As a first pedagogical example we shall focus on a perceptron, a simple two layers module, run against benchmark datasets for classification tasks. 2.2 Triadic Perceptrons (a) Results on mnist (b) Results on fashion-mnist Figure 2: Results of the perceptrons’ training on mnist and fashion-mnist (via Adam optimizer). We performed learning rate warm-up, followed by a reduce-learning-rate-on-plateau protocol, to guard against a possible dependence of the results on an unlucky hyperparameter choice. The direct space triadic model is plagued by instability and confidence saturation problems, while the standard perceptron lacks in expressivity. At variance, the spectral triadic model (i) achieves the same performances of the standard triadic network, with a substantially more efficient parameter scaling and (i) it is also way more stable. To assess the fundamental performance of the shonn framework, we first conduct tests on a single layer perceptron, operated in classification mode against the mnist dataset [20, 15]. The results are reported in Figure 2(a), and show clear performance and stability gaps, favoring the spectral triadic model. Note that for these, and subsequent, experiments we trained 33 separate models for each architecture, to acquire an uncertainty measure. To ensure the robustness of our findings we evaluated the perceptron models on fashion-mnist [24]. The results, displayed in Figure 2(b), corroborate our previous conclusions. Even from these preliminary tests, it is clear that the spectral paradigm not only substantially improves parameter scaling but also mitigates the numerical instability and gradient saturation issues inherent in standard higher-order neural networks trained in direct space. In the next section we will turn to considering the multi-layered generalization of the two layered models. 2.3 Triadic mlps Figure 3: Standard mlp vs. spectral triadic mlp on CIFAR-10. The models were trained with the Adam optimizer, following a halving-lr-on-plateau scheduler. The results show a clear advantage for the triadic architecture. Figure 4: Standard mlp-mixer vs a spectral triadic version of it, on CIFAR-10. For this experiment the Adam optimizer was used with fixed learning rate. From the results, it seems that the spectral forward propagation not only shows a sufficient degree of expressivity, but also acts as an implicit regularizer, preventing overfitting. Perceptrons can be staked recursively to create Multi Layer Perceptrons with increased expressivity power. The same applies to triadic perceptrons of the type analyzed above. A standard 33-layer neural network is the simplest mlp architecture yk=∑iHwki(2)σ(∑jwij(1)xj),y_k= _i^Hw_ki^(2)σ ( _jw_ij^(1)x_j )\,, (8) where wij(1),wki(2)w_ij^(1),w_ki^(2) are the elements of the weight matrices of the first and second layer respectively, and σ is an appropriate non-linear filter (usually a ReLU function) to potentiate the degree of expressivity beyond the trivial linear setting. It can be proven that (8) approximate any continuous function with arbitrary precision as long as the hidden dimension H is large enough, given the complexity of the task [16]. This high degree of expressiveness, however, entails an interpretability trade-off; it is notoriously difficult to elucidate the inner workings of this class of models. From (2), it is evident that an elementary triadic perceptron can represent quadratic mappings, offering a clear advantage in expressivity over standard perceptron models. Furthermore, by hierarchically stacking triadic perceptrons a multi-layer configuration reveals a notable progression in the model’s hypothesis space, e.g., a 22-layer triadic mlp is able to represent polynomial mappings up to the fourth degree. In this sense, it proves equivalent to a higher-order perceptron with pentadic interactions (5-body interactions, a further step up in the higher-order hierarchy), but with a significantly lowered parameter count. It can be shown that a N-layer triadic mlp is able to represent any polynomial with degree 2N−12^N-1 (see Methods section for more details). Building upon this result, we demonstrate that any continuous function can be represented by a sufficiently wide and deep triadic mlp. This holds for both the direct-space parametrization in Eq. (2) (see Appendix A) and the spectral parametrization, as discussed in Section 4.3. Recall that this latter setting comes with a concomitant compression of the relevant parameter space. We also remark that the incremental expressivity, as gained by progressively adding more triadic layers, could enhance the interpretability of the proposed framework. To assess empirically the capabilities of spectral higher-order mlps, we tested their performances on CIFAR-10 [14]. Specifically, we compared the run on a standard, 44 layers, mlp architecture, equipped with a non-linear activation function (ReLU), against an equivalently shaped spectral triadic mlp. The results are shown in Figure 3, and remarkably indicate better performances for our novel architecture, despite the absence of a point-wise non topological activation function σ. Having established the reliability of the triadic mlp, it is important to highlight its potential for integration into existing architectures. One may in fact consider replacing conventional mlp modules with triadic layers to effectively boost model capabilities. To substantiate this claim, we considered a state-of-the-art architecture for vision, the mlp-mixer [8], which we modified by replacing mlp modules with triadic spectral perceptrons. Even in such vastly different context, spectral higher-order networks exhibit a significant performance advantage, as clearly depicted in Figure 4. We refer the reader to Appendix C for extended regression analysis, which confirms the above conclusions. 3 Discussion We here introduced the spectral parameterization for higher-order, triadic, forward propagation. This new framework opens up the perspective to incorporate hyperlinks into feedforward neural networks, without sacrificing the customary O(N2)O(N^2) scaling. Indeed, to the best of our knowledge, the proposed spectral parametrization exhibits the most favorable parameter scaling among all existing higher-order neural network architectures. Our experiments demonstrate that the newly proposed hyper-linked neural architectures exhibit high efficacy even in the absence of explicit, pointwise activation functions, while simultaneously mitigating the stability issues typically associated with triadic configurations trained in direct space. The introduced class of networks also possess a highly configurable hypothesis space: indeed, the maximum complexity of the expressed mapping can be tailored by adjusting the network depth. This tunability can be extremely useful in context where the polynomial order of the targeted mapping is approximately known. For this reason, we identify the application of these networks to physical contexts (e.g. via pinns [18]) as a key future research direction. Furthermore, the results referred to deep architectures (mlps, mlp-mixers) already indicate that higher-order spectral networks outperform their standard, direct space, counterparts. At last, we have derived universal approximation theorems for this class of triadic architectures, rigorously granting their expressive power. The collective bulk of theoretical and empirical results here presented clearly shows that elements of hypergraph theory can be fruitfully integrated into deep learning frameworks to enhance model performance. Future research is required to comprehensively evaluate the full potential of this interdisciplinary framework. 4 Methods 4.1 The Spectral Parametrization Assume a weight matrix WK×NW_K× N of the links between neural layers of size N and K respectively. The spectral parametrization approach [12] focuses on the implied adjacency matrix between the N+KN+K neurons of the underlying bipartite graph: A=[N×N×KWK×NK×K]A= bmatrixO_N× N&O_N× K\\ W_K× N&O_K× K bmatrix (9) The A matrix operates the forward pass on the global activation vector of the bipartite graph, usually called a. Following the action of A on a, the signal that interests the neurons of the last layer is eventually stocked in the final K entries of a. The forward pass of information mediated by A is independent of its diagonal elements. This allows us to define an augmented version of A with non-zero diagonal entries, which serves as the fundamental operator for implementing the spectral parametrization. Specifically we can express the modified adjacency matrix in the form A=ΦΛΦ−1A= ^-1, with Λ∈ℝ(N+K)×(N+K) ^(N+K)×(N+K) diagonal eigenvalues matrix, and with Φ : Φ=[N×N×KϕK×NK×K] = bmatrixI_N× N&O_N× K\\ _K× N&I_K× K bmatrix (10) where I denotes the identity matrix. At last, given the properties of Φ−1 ^-1, is possible to show [12] that the elements wkiw_ki which populate the sub-diagonal block of A can be rewritten as: wki=(λi(in)−λk(out))ϕkiw_ki=( _i^(in)- _k^(out)) _ki (11) Where λi(in) _i^(in) is the i-th element of the first portion of the eigenvalue diagonal, and λk(out) _k^(out) is the k-th of the last one. This novel parametrization for the network’s weights offers the versatility to achieve a multitude of different goals, spanning pruning techniques, input feature relevance detection, architecture search strategies, and, as we shall see, parameters’ reduction paradigms [3, 1, 2, 17]. 4.2 Spectral Parametrization of Triadic Networks In equation (2) the w~kij w_kij tensor element can be thought as a matrix element, w~ki′ w_ki : here i′i represents an appropriate index which enumerates the set of unordered pairs (i,j), 0≤i≤j<N(i,j),\ 0≤ i≤ j<N in lexicographical order: (0,0)→0,(0,1)→1,…,(0,N−1)→N−1,(1,1)→N,…,(N−1,N−1)→N(N+1)2−1(0,0)→ 0,\ (0,1)→ 1,\ ...,\ (0,N-1)→ N-1,\ (1,1)→ N,\ ...,\ (N-1,N-1)→ N(N+1)2-1 In close form i′i is given by: i′=j−i+∑k=0i−1(N−k).i =j-i+ _k=0^i-1(N-k). (12) The lexicographical bijection (12) allows the deployment of the spectral parametrization for triadic interactions: we first use (12) to map the tridimensional weight tensor [w~kij][ w_kij] into the weight matrix [w~ki′][ w_ki ]; we can then simply express [w~ki′][ w_ki ] via the spectral parametrization. From (2) we get: yk=∑i(λi(in)−λk(out))ϕkixi+∑i′(λ~i′(in)−λ~k(out))ϕ~ki′xi′y_k= _i( _i^(in)- _k^(out)) _kix_i+ _i ( λ_i ^(in)- λ_k^(out)) φ_ki x_i (13) with xi′=xixjx_i =x_ix_j. Finally we can expand i′i back to 22 dimensions: yk=∑i(λi(in)−λk(out))ϕkixi+∑0≤i≤j(λ~ij(in)−λ~k(out))ϕ~kijxixj,y_k= _i( _i^(in)- _k^(out)) _kix_i+ _0≤ i≤ j( λ_ij^(in)- λ_k^(out)) φ_kijx_ix_j, (14) to eventually obtain equation (7). 4.3 Expressiveness of a Triadic Multilayer Perceptron In this section, we examine the expressivity of the triadic mlp and establish a universal approximation result for the proposed architecture. Let us consider a neural network layer with propagation rule from input x to output y given by Eq (2). For sake of clarity and without loss of generality, the input bias term is now (and just in this section) incorporated in the triadic interaction. The following expression thus holds: yk=∑iwkixi+∑i≤jw~kijxixjy_k= _iw_kix_i+ _i≤ j w_kijx_ix_j (15) Further bias terms of subsequent triadic layers will not be considered, just as exposed in Sec. 2.1. In short, we will denote the above transformation with the compact notation =HONN() y=HONN( x). According to Eq. (7), impose the spectral parametrization of the weights wki\w_ki\ and w~kij\ w_kij\ amounts to write: wki w_ki =(λi(in)−λk(out))ϕki = ( _i^(in)- _k^(out) ) _ki (16) w~kij w_kij =(λ~ij(in)−λ~k(out))ϕ~kij, = ( λ_ij^(in)- λ_k^(out) ) φ_kij, where λi(in),λk(out),ϕki,λ~ij(in),λ~k(out) are tunable, _i^(in), _k^(out), _ki, λ_ij^(in), λ_k^(out) are tunable, (17) ϕ~kij are fixed to the initialization values. φ_kij\, are fixed to the initialization values. As already mentioned, this choice enables one to effectively reduce the number of tunable parameters, which now scales as O(N2)O(N^2). In essence, training the spectral version of an higher order neural network of the type defined by Eq. (15) comes with a computational cost (as quantified by parameters’ scaling) comparable to that associated with standard mlp. Starting from this setting we will prove the following Theorem. Theorem 4.1. Given ∈X⊂ℝn x∈ X ^n, the space ℋH of functions h(1,)h(1, x), defined by iteratively composing honn layers of the form of Eq. (15) with the parametrization defined by Eqs. (16), (17), is dense in C(X,ℝn)C(X,R^n), namely the space of real-valued continuous functions on the compact set X. Leveraging the celebrated Stone-Weierstrass Theorem, the proof of Theorem 4.1 can be articulated by showing that: m(x1,⋯,xn)⊆ℋ,P_m(x_1,·s,x_n) \,, (18) where m(x1,⋯,xn)P_m(x_1,·s,x_n) stands for a generic m-degree polynomial in ℝnR^n. In Appendix A we will present also a version of this result in direct space, namely without the constraints imposed by Eqs (16)-(17). In this latter setting, the architecture aligns with those outlined in [5, 4]. While this result may be intuitive, to the best of our knowledge, no formal proof has yet been documented in the literature. Before demostrating Theorem 4.1 in its full generality, we first outline the proof strategy in a simplified case where the input is restricted to a scalar, one-dimensional variable, namely =(x)∈ℝx=(x) , supplemented with a trivial bias term, 11. The proof proceeds by iteratively constructing polynomials of increasing degree through the composition of successive layers in the honn. The first step (see Lemma 4.1) is devoted to the creation of polynomials of second degree. Let us remember that we start with 11, x (and their combinations), hence we pass from a first degree polynomial to polynomials with degrees one step higher, so to increase expressivity. It should be noted that we require these second-order polynomials to be capable of representing any arbitrary second-order polynomial. Our goal is to demonstrate the existence of a parameter configuration for a honn with three output nodes which is capable to generate three linearly independent second-order polynomials, denoted as y0(x)y_0(x), y1(x)y_1(x) and y2(x)y_2(x). The sought claim will be hence indirectly proven. Indeed we will show that a second set of weights, formally a new layer added to the honn, exists such that the standard basis elements of 2(x)P_2(x), i.e., 11, x and x2x^2, can be obtained from linear combinations of yj(x)y_j(x) that make use of the above mentioned weights (see Fig. 5 for a pictorial representation of the claim). Once this step is achieved, we can iterate the construction by building polynomials of degree three upon adding extra layers, and so forth up to any desired degree. Let us conclude this introduction with two remarks. First, the construction here presented can be extended to build polynomials of arbitrary degree in n variables (see Fig. 6 for a pictorial representation and Lemma 4.2 and 4.3 for a formal proof). Second, the mechanism by which we prove expressivity does not necessarily reflect the configuration produced by an optimization process acting directly on the honn weights. This, however, poses no limitation, as our primary objective is to establish the existence of such expressivity rather than the recovery of the "optimal" parameter setting. In the case of input (1,x)(1,x), the higher-order nonlinear interaction generates a term that scales like x2x^2, in the output function =HONN(1,x) y=HONN(1,x). In the following we will show, that the spectral parameters governing the update rule can be assigned so to map y into a basis for 2(x)=Span1,x,x2P_2(x)=Span\1,x,x^2\. 11xxy0y_0y1y_1y2y_211x2x^2W−1W^-1 Figure 5: Graphical sketch of Lemma 4.1, illustrating the inter-layer connections (1,x)→(1,x)→ y described in Eq. (21), under the triangular ansatz and in the simple setting where the only nonlinear interactions arise from the last column of matrix W. Subsequently, the elements of y can be linearly remapped to the canonical basis of 2(x)P_2(x). The black dashed arrows represent the linear passage corresponding to the first sum term of Eq. (15), while the orange solid lines represent the higher-order interactions, namely second sum in Eq. (15), here just restricted to x and itself. Lemma 4.1. Consider the output of a spectral layer of the form stipulated by Eq. (15), with the parametrization given by Eqs. (16)-(17) applied to a one-dimensional input x (together with the bias term), namely: =HONN(1,x), y=HONN(1,x), (19) where =(y0,y1,y2) y=(y_0,y_1,y_2). Hence, a solution exists (in terms of the trained spectral parameters) such that 2(x)=Span1,x,x2⊆Spany0,y1,y2,P_2(x)=Span\1,x,x^2\ \y_0,y_1,y_2\, i.e., the elements yi∈2(x)y_i _2(x) are linearly independent. Proof. We begin by writing the explicit expressions for y0,y1,y2y_0,y_1,y_2: y0 y_0 =(w00+w~000) 1+(w01+w~001)x+w~011x2, =(w_00+ w_000)1+(w_01+ w_001)\,x+ w_011\,x^2, (20) y1 y_1 =(w10+w~100) 1+(w11+w~101)x+w~111x2, =(w_10+ w_100)1+(w_11+ w_101)\,x+ w_111\,x^2, y2 y_2 =(w20+w~200) 1+(w21+w~201)x+w~211x2. =(w_20+ w_200)1+(w_21+ w_201)\,x+ w_211\,x^2. From Eq. (20) we can define the associated matrix W acting from the space Span1,x,x2Span\1,x,x^2\ to Spany0,y1,y2Span\y_0,y_1,y_2\: W=(w00+w~000w01+w~001w~011w10+w~100w11+w~101w~111w20+w~200w21+w~201w~211).W= pmatrixw_00+ w_000&w_01+ w_001& w_011\\ w_10+ w_100&w_11+ w_101& w_111\\ w_20+ w_200&w_21+ w_201& w_211 pmatrix. (21) Imagine that one can find a (at least one) specific combination of the weights, under the imposed spectral parametrization, such that det(W)≠0 (W)≠ 0. Then, we can prove that y0,y1,y2\y_0,y_1,y_2\ generates the space 2(x)P_2(x), as claimed, and consequently invert the above matrix to remap the output into the canonical basis through a linear layer whose weights incorporate W−1W^-1 (see Fig. 5 for a visual representation of the proposed procedure). To reach this goal, we show that W can be forced into an upper triangular form, with non zero elements across the main diagonal, which in turn amounts to set out det(W)≠0 (W)≠ 0, as wished. More concretely, we will demonstrate that the spectral parametrization is sufficiently flexible to turn to zero all the entries under the diagonal of matrix (21), while keeping the diagonal terms nonzero. In formulae: w10=−w~100⟹(λ0(in)−λ1(out))ϕ10=−w~100⟹ϕ10=−w~100λ0(in)−λ1(out),λ0(in)≠λ1(out), w_10=- w_100 ( _0^(in)- _1^(out) ) _10=- w_100 _10=- w_100 _0^(in)- _1^(out),\> _0^(in)≠ _1^(out), (22) w20=−w~200⟹(λ0(in)−λ2(out))ϕ20=−w~200⟹ϕ20=−w~200λ0(in)−λ2(out),λ0(in)≠λ2(out), w_20=- w_200 ( _0^(in)- _2^(out) ) _20=- w_200 _20=- w_200 _0^(in)- _2^(out),\> _0^(in)≠ _2^(out), w21=−w~201⟹(λ1(in)−λ2(out))ϕ21=−w~201⟹ϕ21=−w~201λ1(in)−λ2(out),λ1(in)≠λ2(out), w_21=- w_201 ( _1^(in)- _2^(out) ) _21=- w_201 _21=- w_201 _1^(in)- _2^(out),\> _1^(in)≠ _2^(out), and w00≠−w~000⟹(λ0(in)−λ0(out))ϕ00≠−w~000⟹ϕ00≠−w~000λ0(in)−λ0(out), w_00≠- w_000 ( _0^(in)- _0^(out) ) _00≠- w_000 _00≠- w_000 _0^(in)- _0^(out), (23) w11≠−w~101⟹(λ1(in)−λ1(out))ϕ11≠−w~101⟹ϕ11≠−w~101λ1(in)−λ1(out), w_11≠- w_101 ( _1^(in)- _1^(out) ) _11≠- w_101 _11≠- w_101 _1^(in)- _1^(out), w~211≠0⟹(λ~11(in)−λ~2(out))ϕ~211≠0⟹λ~2(out)≠λ~11(in)andϕ~211≠0. w_211≠ 0 ( λ_11^(in)- λ_2^(out) ) φ_211≠ 0 λ_2^(out)≠ λ_11^(in)\>\>and\>\> φ_211≠ 0. The conditions stemming from Eqs. (22) and (23) can be straightforwardly fulfilled by (i) setting ϕ~211≠0 φ_211≠ 0. Notice that ϕ~211 φ_211 is a non tunable parameter and thus it gets permanently frozen to the value assigned following initialization. (i) exploiting the fact that ϕij _ij can be tuned at will to set them to the values prescribed in Eq. (17). Under the above mentioned conditions, matrix W is reduced to a upper triangular form with non zero diagonal entries. As such, it can be inverted which immediately yields that Spany0,y1,y2=2(x)Span\y_0,y_1,y_2\=P_2(x), namely the proof of the Lemma. ∎ The proof can be adapted easily to deal with the subsequent step that aims at generating 3(x)P_3(x) from 2(x)P_2(x), with the inclusion of one additional cubic term. By iterating forward the process, with the inclusion of just one non-linear monomial per step, it readily yields m(x)P_m(x), for any generic m degree. Instead of detailing the derivation of the above results with reference to the scalar setting, we shall turn to discuss the general case of polynomial m()P_m( x), where ∈ℝ x as addressed in Theorem 4.1. First, we start by computing the number D of monomials which are associated to a complete polynomial of degree m in ℝnR^n. This is given by: D=∑k=0m(n+k−1k)=∑k=0m(n+k−1n−1)=∑k=n−1n+m−1(kn−1)=(m+n),D= _k=0^m n+k-1k= _k=0^m n+k-1n-1= _k=n-1^n+m-1 kn-1= m+nn, (24) where in the last chain of identities of Eq. (24) we used the hockey stick identity. Then, to prove the Theorem we proceed with a ladder strategy that we outline in the following. Given e¯0,⋯,e¯q\ e_0,·s, e_q\, q+1q+1 elements of the canonical basis of m(x1,⋯,xn)P_m(x_1,·s,x_n), we provide an iterative recipe to yield a novel and independent element e¯q+1 e_q+1 until reaching the sought dimension D. The first step of the ladder strategy (or the initial level of an induction proof) refers to the input features (potentiated with the bias term): 1,x1,⋯,xn1,x_1,·s,x_n, which also belong to the canonical basis of m(x1,⋯,xn)P_m(x_1,·s,x_n). The following Lemma holds: 11x1x_1x2x_2y0y_0y1y_1y2y_2y3y_311x1x_1x2x_2x1x2x_1x_2W−1W^-1y0y_0y1y_1y2y_2y3y_3y4y_4⋯·s Figure 6: Graphical sketch of the ladder demonstration scheme proposed in Lemmas 4.2, 4.3. Here the input (1,x1,x2)(1,x_1,x_2) is provided to the honn with output (y0,y1,y2,y3)(y_0,y_1,y_2,y_3), which can be linearly remapped, via W−1W^-1, to the new canonical basis augmented with a new nonlinear term, for instance x1x2x_1x_2. The procedure then continues by creating new terms: x12x_1^2, x22x_2^2, and so on. As in Fig. 5, dashed black lines indicate linear transfer, while orange ones indicate nonlinear interactions. Lemma 4.2. The first step of the induction proof of Theorem 4.1 Given (1,)=(1,x1,⋯,xn)(1, x)=(1,x_1,·s,x_n), then there exists a set of weights wkiw_ki and w~kij w_kij parametrized according to Eqs. (16), (17) such that Span1,x1,⋯,xn,xh¯xk¯⊆Spany0,⋯,yn+1Span\1,x_1,·s,x_n,x_ hx_ k\ \y_0,·s,y_n+1\ where =HONN(1,) y=HONN(1, x). Proof. Consider =(y0,y1,⋯,yn+1)=HONN(1,) y=(y_0,y_1,·s,y_n+1)=HONN(1, x), where HONN stands for a neural network layer defined by Eq. (15) with n+1n+1 nodes, as the input and n+2n+2 nodes for the output. Restrict to the case where the just one non-linear term is produced (at each iteration or, equivalently, added layer) by the interaction among xh¯x_ h and xk¯x_ k, for any generic h¯ h and k¯ k (raging from 11 to n). Then, the general output y can be written as: y0 y_0 =(w00+w~000)1+⋯+(w0j+w~00j)xj+⋯+w~0h¯k¯xh¯xk¯, =(w_00+ w_000)1+·s+(w_0j+ w_00j)x_j+·s+ w_0 h kx_ hx_ k, (25) y1 y_1 =(w10+w~100)1+⋯+(w1j+w~10j)xj+⋯+w~1h¯k¯xh¯xk¯, =(w_10+ w_100)1+·s+(w_1j+ w_10j)x_j+·s+ w_1 h kx_ hx_ k, ⋮ yi y_i =(wi0+w~i00)1+⋯+(wij+w~i0j)xj+⋯+w~ih¯k¯xh¯xk¯, =(w_i0+ w_i00)1+·s+(w_ij+ w_i0j)x_j+·s+ w_i h kx_ hx_ k, ⋮ yn+1 y_n+1 =(wn+1,0+w~n+1,00)1+⋯+(wn+1,j+w~n+1,0j)xj+⋯+w~n+1,h¯k¯xh¯xk¯, =(w_n+1,0+ w_n+1,00)1+·s+(w_n+1,j+ w_n+1,0j)x_j+·s+ w_n+1, h kx_ hx_ k, with an associated matrix, connecting y to the new basis 1,x1,⋯,xn,xh¯xk¯\1,x_1,·s,x_n,x_ hx_ k\, given by: W=(w00+w~000⋯w0j+w~00j⋯w~0h¯k¯⋮⋯wi0+w~i00⋯wij+w~i0j⋯w~ih¯k¯⋮⋯wn+1,0+w~n+1,00⋯wn+1,j+w~n+1,0j⋯w~n+1,h¯k¯).W= pmatrixw_00+ w_000&·s&w_0j+ w_00j&·s& w_0 h k\\ &&·s&&·s\\ w_i0+ w_i00&·s&w_ij+ w_i0j&·s& w_i h k\\ &&·s&&·s\\ w_n+1,0+ w_n+1,00&·s&w_n+1,j+ w_n+1,0j&·s& w_n+1, h k pmatrix. (26) As shown with reference to the simple scalar setting examined above, it is sufficient to show that a set of spectral weights exists such that det(W)≠0 (W)≠ 0. This ensures that y0,⋯,yn+1\y_0,·s,y_n+1\ generates Span1,x1,⋯,xn,xh¯xk¯Span\1,x_1,·s,x_n,x_ hx_ k\, namely that the produced output can be remapped onto the expanded canonical basis via a linear transformation associated to W−1W^-1 (see Fig. 6). Following the reasoning outlined above, we shall set to prove that a specific choice of the parameters exists that turns W into an upper triangular matrix, with non zero diagonal elements. Specifically, we seek to set to zero all the terms wij,i>jw_ij,\,i>j, while imposing wij≠0,i=j.w_ij≠ 0,\,i=j. To this end write: wij+w~i0j=0,ifi>j, w_ij+ w_i0j=0,\>if\>\>i>j, (27) wij+w~i0j≠0,ifi=j, w_ij+ w_i0j≠ 0,\>if\>\>i=j, w~n+1,h¯k¯≠0. w_n+1, h k≠ 0. The first of the above conditions, expressed in terms of the spectral parametrization Eq. (16), yields: (λj(in)−λi(out))ϕij=−w~i0j⟹ϕij=−w~i0jλj(in)−λi(out),λj(in)≠λi(out),∀i>j, ( _j^(in)- _i^(out) ) _ij=- w_i0j _ij=- w_i0j _j^(in)- _i^(out),\> _j^(in)≠ _i^(out),\>∀\,i>j, (28) The second set of conditions, prescribed in Eq. (27), results in (λj(in)−λi(out))ϕij≠−w~i0j⟹ϕij≠−w~i0jλj(in)−λi(out),∀i=j. ( _j^(in)- _i^(out) ) _ij≠- w_i0j _ij≠- w_i0j _j^(in)- _i^(out),\>∀\,i=j. (29) Both requirements (28), (29) can be easily fulfilled since the spectral eigenvector entries ϕij _ij (together with the associated eigenvalues) are fully trainable. The last condition of Eq. (27) reads, w~n+1,h¯k¯≠0⟹(λ~h¯k¯(in)−λ~n+1(out))ϕ~n+1,h¯k¯≠0⟹λ~n+1(out)≠λ~h¯k¯(in). w_n+1, h k≠ 0 ( λ_ h k^(in)- λ_n+1^(out) ) φ_n+1, h k≠ 0 λ_n+1^(out)≠ λ_ h k^(in). (30) This is the only condition imposed on the trainable higher-order spectral parameters (together with initializing ϕ~n+1,h¯k¯≠0 φ_n+1, h k≠ 0) and, as such, it can be accommodated for without clashes in the process of nested parameters handling. This ends the proof of the Lemma. ∎ As we can see from Eq. (28) the parameters wi00w_i00 are free to be set to zero, recovering the numerical setup where the biases are not included in the higher-order sum of Eq. (15). We can proceed forward to prove the following Lemma, which constitutes the main inductive step of the proof. In turn, the following proof combined with the preceding Lemma, will return a complete proof of the sought Theorem 4.1. Lemma 4.3. The general step of the induction proof of Theorem 4.1 Given e¯0,⋯,e¯q\ e_0,·s, e_q\, with q<D−1q<D-1 and e¯0≡1,e¯1≡x1,⋯,e¯n≡xn, e_0≡ 1,\, e_1≡ x_1,\,·s,\, e_n≡ x_n,. Then there exists a set of weights wkiw_ki and w~kij w_kij parametrized according to Eqs. (16), (17) such that Span1,e¯1,⋯,e¯q,e¯h¯e¯k¯⊆Spany0,⋯,yq+1Span\1, e_1,·s, e_q, e_ h e_ k\ \y_0,·s,y_q+1\ where =HONN(1,¯) y=HONN(1, e). Proof. The scheme applied to the proof of Lemma 4.2 can be straightforwardly generalized to deal with a set of monomials e¯0,⋯,e¯q\ e_0,·s, e_q\, with q<D−1q<D-1. In particular, e¯0≡1,e¯1≡x1,⋯,e¯n≡xn, e_0≡ 1,\, e_1≡ x_1,\,·s,\, e_n≡ x_n, and define e¯q+1≡eh¯ek¯ e_q+1≡ e_ he_ k, for any integers h¯ h and k¯ k (smaller than q). It is immediate to show that this setting maps exactly onto the one analyzed above in Lemma 4.2. In other words, the produced output can be turned into a canonical basis via a linear transformation that is granted by the inverse of an upper triangular matrix with non zero determinant. One can thus generate, any monomial q¯ q in the range q<q¯<Dq< q<D, as an iterative sequence of the algorithm described above. ∎ We have thus proved by induction that a basis for m(x1,⋯,xn)P_m(x_1,·s,x_n) can be iteratively constructed by an ad hoc choice of the spectral parameters, thus implying that, sufficiently large, spectral higher order neural networks of the type here considered, possess the expressivity required to approximate arbitrary polynomials of any degrees m and dimensions n. An illustrative scheme to depict the strategy employed in the proof is displayed in Figure 6. The function h()h( x), as obtained by stacking an arbitrary number of layers (with the needed number of nodes to store the basis elements), can thus mimic a general polynomial of arbitrary degree. Therefore, thanks to the Stone-Weierstrass theorem, it can approximate a target continuous function defined on a compact set to the desired level of approximation. Appendix A Proof of the expressiveness in Direct Space In the following we shall denote by k(ℝn)P_k(R^n) the vector space of all polynomials of n variables, up to k-th degree; we will also use the notation m(h)m^(h) to refer to an element of degree h, of the standard basis, i.e., composed by the standard monomials of k(ℝn)P_k(R^n). Before delving into the general theorem on the expressiveness of (non-spectral) deep triadic mlps, it is instructive to consider the basic case of a triadic network with just one input neuron. Given a scalar input x as a first key observation we note that the forward pass, see Eq. (2) in the main text, allows the model to output the standard basis of the space 2(ℝ)P_2(R), in the second layer (provided enough neurons are supplied). Indeed, if one equips the second layer with 33 neurons, it is manifestly possible to set the weights wki,w~kijw_ki, w_kij so as to obtain output (proportional to) 11 on the first neuron of the second layer, (proportional to) x on the second, and (proportional to) x2x^2 on the third. We have thus recovered the standard basis of 2(ℝ)P_2(R). The second key observation consists in noticing that this process can be iterated across subsequent layers: given a neural layer containing the standard basis of the space k(ℝ)P_k(R), the forward pass that follows Eq. (2) can generate the standard basis of the space 2k(ℝ)P_2k(R) in the next layer (once again, if enough neurons are present). The third, and final, observation is that, given the standard basis of k(ℝ)P_k(R), the forward pass dictated by Eq. (2) is clearly able to output any polynomial function at least up to degree k in the successive layer, no matter the size of the target output vector. These three observations, once put together, clearly imply that a single-input triadic mlp, if large enough, can regress any polynomial function of x. By leveraging on the Stone–Weierstrass theorem, we can thus conclude that a sufficiently large, single-input, triadic mlp can approximate any continuous function of x, even in the absence of a non-linear activation function. This result amounts to an universal approximation theorem for triadic neural networks with one-dimensional inputs, in direct space. A simple graphical sketch of the provided arguments is displayed in Fig. 7. x11x2x^211x2x^2x3x^3x4x^4p(4+)p^(4+)p(4+)p^(4+)w0,1,w~0,1,1=0w_0,1, w_0,1,1=0w1,0,w~1,1,1=0w_1,0, w_1,1,1=0w2,0,w2,1=0w_2,0,w_2,1=0 Figure 7: Graphical sketch of the universality proof for triadic multilayer perceptrons. The structure of the forward pass (2) allows our network to build standard basis of the polynomial vector spaces. Specifically, in this toy example it is easy to see that a sufficiently large k-th layer can yield the standard basis of the vector space 2k−1(ℝ)P_2^k-1(R). This in turn implies that the subsequent layer will be able to generate polynomial functions f:ℝ→ℝnf:R ^n at least up to the 2k−12^k-1 order. We keep track of the inherited level of expressivity with the p(∘+)p^( +) node tag (in the example in figure k=3k=3). The above qualitative description can be put on solid mathematical ground thanks to Theorem A.1, which holds true in the general setting of n variables. As a preliminary result, let us state and prove the following Lemma: Lemma A.1. A triadic perceptron, with a single output neuron, processing as an input the standard base of the vector space k(ℝn)P_k(R^n) can output any monomial m up to order 2k2k. Proof. The proof is established by contradiction. Let us focus on the single neuron of the next layer: suppose it exists a monomial m~ m of order α≤2kα≤ 2k that cannot be generated by the model forward pass, i.e. cannot be written in the form ∑iw0imi(k)+∑i,jw~0ijmi(k)mj(k), _iw_0im_i^(k)+ _i,j w_0ijm_i^(k)m_j^(k), (31) and so it does not belong to the span: m~∉span(mi(k)∪mi(k)mj(k)), m (\m^(k)_i\∪\m^(k)_im^(k)_j\ ), (32) where mi(k)m_i^(k) and mj(k)m_j^(k) are two monomials of degree k (note that the indicies i,ji,j span all the order k standard basis monomials). However, note that one can always write the selected monomial as: m~=m′⋅m′, m=m · m , (33) where m′m and m′m are two monomials of order α′,α′≤kα ,α ≤ k. This is clearly an element of the span of mi(k)mj(k)\m^(k)_im^(k)_j\, thus leading to a contradiction. ∎ Theorem A.1. A sufficiently large triadic multilayer perceptron is able to approximate any continuous function f:ℝn→ℝνf:R^n ^ν with an arbitrary degree of precision. Proof. Let us start by proving that a large enough triadic multilayer perceptron, with input ∈ℝn x ^n, can generate a standard basis for any polynomial vector space k(ℝn)P_k(R^n), regardless of how large k is. This proof is by induction: the base case consists in observing that a triadic perceptron, with input layer of size v, by definition holds in its first layer a standard basis of the vector space 1(ℝn)P_1(R^n) (i.e. x,y,z,…x,y,z,...; and as we argued in previous sections the absence of the constant bias term as neural activation is irrelevant). For the inductive step we want to show that, given a collection of monomials mi(k)\m_i^(k)\ that forms the standard basis for the vector space k(ℝn)P_k(R^n), a triadic perceptron gathering such a collection as an input is able to output, at the next layer, the standard basis of the vector space 2k(ℝn)P_2k(R^n) (once again, assuming the next layer has enough neurons). Thanks to Lemma A.1, we know that the proposed forward pass is able to generate any monomial up to order 2k2k as output. Since the lemma applies to each next layer neuron (with new weights each time), it is sufficient to prove that the triadic perceptron can to generate any collection of monomials up to order 2k2k as an output, i.e. it is able to generate the standard basis of the polynomial vector space 2k(ℝn)P_2k(R^n), thus concluding the proof of the inductive step. A sufficiently large triadic multilayer perceptron is able to generate any standard basis mi(s)\m^(s)_i\ of the polynomial vector space, regardless of the value of s (maximum order) or n (number of variables, i.e input size). This in turn implies that, by adding a new neural layer on top of the existing ones, we will be able to represent any polynomial function p:ℝn→ℝνp:R^n ^ν up to order s, with ν size of the lastly added layer. Given that there is no bound to the order s (it can be made as large as wished by adding new layers), this means that a sufficiently large triadic multilayer perceptron is able to represent any polynomial function, of any order. This at last allows us to invoke the Stone–Weierstrass theorem to finally prove that the model, if given enough neurons in enough layers, is capable to approximating, with an arbitrary degree of precision, any continuous function. ∎ Appendix B Parameter scaling comparison with state-of-the-art models The novel hypergraph architecture can be manifestly re-interpreted as a deep polynomial architecture. In the literature similar models have been discussed [5, 4], but remarkably the proposed spectral variant yields a smaller parameter scaling when compared with other state-of-the-art polynomial architectures, as we shall see in this section. For a perceptron with N input and output neurons, a spectral triadic network has the following parametrization yk=∑iwkixi+∑i≤j(λ~ij(in)−λ~k(out))ϕ~kijxixj,y_k= _iw_kix_i+ _i≤ j( λ_ij^(in)- λ_k^(out)) φ_kijx_ix_j, (34) note that the output is a polynomial of order 22. Since the generalized eigenvectors ϕ~kij φ_kij are not trained this implies the handling of two parameter matrices (W,Λ~(in)W, ^(in)), and a single parameter vector (~(out) λ^(out)). The number of involved parameters is thus: Pspectral(N)=N(N+1)⏟W+N(N+1)2⏟Λ~(in)+N⏟~(out)=32N2+52N.P_spectral(N)= N(N+1)_W+ N(N+1)2_ ^(in)+ N_ λ^(out)= 32N^2+ 52N. (35) If instead we choose to parametrize spectrally also the linear transfer we get: yk=∑i(λi(in)−λk(out))ϕkixi+∑i,j(λ~ij(in)−λ~k(out))ϕ~kijxixj.y_k= _i( _i^(in)- _k^(out)) _kix_i+ _i,j( λ_ij^(in)- λ_k^(out)) φ_kijx_ix_j. (36) Notice that the triadic part is unaltered, but instead of parameterizing the linear transfer with a single weight matrix W we are now using one matrix (ϕφ) and two vectors ((in),(out) λ^(in), λ^(out)). This, of course, implies the following parameter scaling: Pfullspectral(N)=3N+1⏟(in),(out),~(out)+N(N+1)⏟ϕ+N(N+1)2⏟Λ~(in)=32N2+92N+1P_fullspectral(N)= 3N+1_ λ^(in), λ^(out), λ^(out)+ N(N+1)_φ+ N(N+1)2_ ^(in)= 32N^2+ 92N+1 (37) Notably, this resultant scaling remains superior to the most parsimonious deep polynomial architectures currently documented in the literature. This is Pccp(N)=3N2+NP_ccp(N)=3N^2+N for the ccp decomposition version of the ProdPoly network [5], as elaborated in the next sub-section. B.1 CCP Decomposition In the CCP decomposition a polynomial of order 22 is instead obtained via the following: =+CU[1]T+C(U[2]⊙U[1])T(⊙) y= β+CU^T_[1] x+C(U_[2] U_[1])^T( x x) (38) where ⊙ indicates the Khatri-Rao product. Given, as before, an input dimension and an output dimension equal to N, this new parametrization introduces a degree of ambiguity since a new hidden dimension K is present. In fact, C is a N×KN× K matrix, and U[i]U_[i] are all N×KN× K matrices as well. This new variable could imply an increment or a reduction of the number of trainable parameters. For a fair comparison, we just take K to be equal to N, and this yields the following parameter scaling for the ccp decomposition: Pccp(N)=3N2+NP_ccp(N)=3N^2+N (39) The result is indeed expected since we deal with three N×N× N matrices (C,U[1]T,U[2]TC,U^T_[1],U^T_[2]), and a vector β of length N. Appendix C Regression Figure 8: Heatmap comparison between the function f(x,y)=x4+y4+2x2y2−2x2f(x,y)=x^4+y^4+2x^2y^2-2x^2 (on the left) and the reconstructed data by a triadic spectral mlp with hidden dimension hdim=20h_dim=20 and nℓ=3n_ =3 number of layers (by including input and output) and trained with the Adam optimizer and fixed learning rate. (a) (b) Figure 9: Regression of experimental data, opportunely rescaled, representing a spiking neural-like behaviour of a relaxation nonlinear circuit [10, 9]. On the left, we can see the charge of a capacitor indicated as x(t)x(t), fully reconstructed while on the right the regression of a spiking current denoted by y(t)y(t). This fit has been carried out by implementing a triadic mlp with hidden dimension hdim=20h_dim=20 and nℓ=20n_ =20 number of layers (including input and output) and trained with the Adam optimizer, following a halving-lr-on-plateau scheduler. Here, self-coupling pairs are not included into the implementation. Here we report on a gallery of regression tests based on the proposed architectures, both on mock and real data. These results serve primarily as a visual demonstration of the expressive power of the proposed network. A general theorem concerning universality is presented in Sec. 4.3. We first regressed the double-well profile f(x,y)=x4+y4+2x2y2−2x2f(x,y)=x^4+y^4+2x^2y^2-2x^2 with a triadic spectral mlp (see Fig. 8). As shown in the figure, the target function (depicted in the left panel) is accurately reconstructed (right panel), with a test mae of 9.5×10−69.5× 10^-6. We also present a fit of experimental data of a two-dimensional dynamical system representing a nonlinear circuit exhibiting neural-like spiking behavior (see [10, 9]). As it can be appreciated by visual inspection, a triadic mlp implementation is able to closely track the abrupt behavioral changes as displayed by the recorded data and reproduce the spiking dynamics, with a mean average error mae 8×10−38× 10^-3. Code Availability The code used for the tests reported in this paper is available at the following link: https://github.com/gianluca-peri/hyperspectral References [1] L. Buffoni, E. Civitelli, L. Giambagli, L. Chicchi, and D. Fanelli (2022) Spectral pruning of fully connected layers. Scientific Reports 12 (1), p. 11201. Cited by: §4.1. [2] L. Chicchi, L. Buffoni, D. Febbe, L. Giambagli, R. Marino, and D. Fanelli (2024) Estimating global input relevance and enforcing sparse representations with a scalable spectral neural network approach. arXiv preprint arXiv:2406.01183. Cited by: §4.1. [3] L. Chicchi, L. Giambagli, L. Buffoni, T. Carletti, M. Ciavarella, and D. Fanelli (2021) Training of sparse and dense deep neural networks: fewer parameters, same performance. Physical Review E 104 (5), p. 054312. Cited by: §2.1, §4.1. [4] G. G. Chrysos, M. Georgopoulos, J. Deng, J. Kossaifi, Y. Panagakis, and A. Anandkumar (2022) Augmenting deep classifiers with polynomial neural networks. In European Conference on Computer Vision, p. 692–716. Cited by: Appendix B, §1, §4.3. [5] G. G. Chrysos, S. Moschoglou, G. Bouritsas, J. Deng, Y. Panagakis, and S. Zafeiriou (2021) Deep polynomial neural networks. IEEE transactions on pattern analysis and machine intelligence 44 (8), p. 4021–4034. Cited by: Appendix B, Appendix B, §1, §4.3. [6] G. E. H. e. al. David E. Rumelhart (1986) Parallel distributed processing: explorations in the microstructure of cognition, vol. 1: foundations. The MIT Press. Cited by: §1. [7] S. Ebli, M. Defferrard, and G. Spreemann (2020) Simplicial neural networks. ArXiv. Cited by: §1. [8] I. O. T. et al. (2021) MLP-mixer: an all-mlp architecture for vision. In Neural Information Processing Systems, Cited by: §2.3. [9] D. Febbe, A. Di Garbo, R. Mannella, R. Meucci, and D. Fanelli (2024) Chaos and synchronization in the ujt relaxation oscillator. In 2024 IEEE Workshop on Complexity in Engineering (COMPENG), p. 1–7. Cited by: Figure 9, Figure 9, Appendix C. [10] D. Febbe, R. Mannella, R. Meucci, and A. Di Garbo (2024) Dynamical behaviour of a new model for the ujt relaxation oscillator. Chaos, Solitons & Fractals 183, p. 114906. Cited by: Figure 9, Figure 9, Appendix C. [11] Y. Feng, H. You, Z. Zhang, R. Ji, and Y. Gao (2018) Hypergraph neural networks. In AAAI Conference on Artificial Intelligence, Cited by: §1. [12] L. Giambagli, L. Buffoni, T. Carletti, W. Nocentini, and D. Fanelli (2021) Machine learning in spectral domain. Nature Communications 12 (1), p. 1330. Cited by: §1, §4.1, §4.1. [13] C. Giles and T. Maxwell (1987) Learning, invariance, and generalization in high-order neural networks.. Applied optics 26 23, p. 4972–8. External Links: Link Cited by: §1, §2.1. [14] A. Krizhevsky (2009) Learning multiple layers of features from tiny images. In Technical Report, University of Toronto, Cited by: §2.3. [15] Y. LeCun, C. Cortes, and C. Burges (1998) The mnist database of handwritten digits. Note: http://yann.lecun.org/exdb/mnist/ Cited by: §2.2. [16] M. Leshno, V. Ya. Lin, A. Pinkus, and S. Schocken (1993) Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. Neural Networks 6 (6), p. 861–867. Cited by: §2.3. [17] G. Peri, L. Chicchi, D. Fanelli, and L. Giambagli (2025-12-01) SPectral architecture search for neural network models. npj Artificial Intelligence 1 (1), p. 43. External Links: ISSN 3005-1460, Document, Link Cited by: §4.1. [18] M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 378, p. 686–707. Cited by: §3. [19] S. Rendle (2010) Factorization machines. 2010 IEEE International Conference on Data Mining, p. 995–1000. Cited by: §1. [20] F. Rosenblatt (1958) The perceptron: a probabilistic model for information storage and organization in the brain.. Psychological review 65 6, p. 386–408. Cited by: §2.2. [21] N. Shazeer (2020) GLU variants improve transformer. Cited by: §1. [22] Y. Shin and J. Ghosh (1991) The pi-sigma network: an efficient higher-order neural network for pattern classification and function approximation. IJCNN-91-Seattle International Joint Conference on Neural Networks i, p. 13–18 vol.1. Cited by: §1. [23] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin (2017) Attention is all you need. In Neural Information Processing Systems, Cited by: §1. [24] H. Xiao, K. Rasul, and R. Vollgraf (2017) Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms. ArXiv abs/1708.07747. External Links: Link Cited by: §2.2.