Deep Learning Alternatives Of The Kolmogorov Superposition Theorem

## Abstract

Abstract This paper explores alternative formulations of the Kolmogorov Superposition Theorem (KST) as a foundation for neural network design. The original KST formulation, while mathematically elegant, presents practical challenges due to its limited insight into the structure of inner and outer functions and the large number of unknown variables it introduces. Kolmogorov-Arnold Networks (KANs) leverage KST for function approximation, but they have faced scrutiny due to mixed results compared to traditional multilayer perceptrons (MLPs) and practical limitations imposed by the original KST formulation. To address these issues, we introduce ActNet, a scalable deep learning model that builds on the KST and overcomes some of the drawbacks of Kolmogorov’s original formulation. We evaluate ActNet in the context of Physics-Informed Neural Networks (PINNs), a framework well-suited for leveraging KST’s strengths in low-dimensional function approximation, particularly for simulating partial differential equations (PDEs). In this challenging setting, where models must learn latent functions without direct measurements, ActNet consistently outperforms KANs across multiple benchmarks and is competitive against the current best MLP-based approaches. These results present ActNet as a promising new direction for KST-based deep learning applications, particularly in scientific computing and PDE simulation tasks.

## 1 Introduction

The Kolmogorov Superposition Theorem (KST) is a powerful and foundational result in mathematics, originally developed to solve Hilbert’s 13th problem (Kolmogorov, ). Over time, it has gained significant attention as a tool for function approximation, showing that any multivariate continuous function can be represented exactly using finite sums and compositions of univariate functions. This theoretical insight has inspired numerous applications in mathematics and computational sciences since its inception in 1957, particularly in problems requiring efficient representation of complex functions in low-dimensional domains (Sprecher, ; Kŭrková, ; Köppen, ).

Recently, Kolmogorov-Arnold Networks (KANs) (Liu et al., ) have renewed interest in the practical use of KST, especially in the context of physics-informed machine learning (PIML) (Karniadakis et al., ). Several papers have been published on this topic in the span of just a few months (Shukla et al., ; Wang et al., ; Howard et al., ; Shuai & Li, ; Koenig et al., ; Rigas et al., ; Patra et al., ), demonstrating the current relevance of the KST for PIML research.

The KST’s strength in approximating functions aligns well with problems in scientific computing and partial differential equation (PDE) simulation, where we often encounter functions representing physical fields (e.g., velocity, temperature, pressure) that live in low-dimensional domains. Physics-Informed Neural Networks (PINNs) (Raissi et al., ; Wang et al., ) have emerged as a popular deep learning framework for such problems, offering a unique approach to learning latent functions without direct measurements by minimizing PDE residuals. By exploiting the KST, KANs offer a mathematically grounded approach to function approximation, sparking new research and development in the area of neural networks for numerically solving PDEs.

Despite its mathematical elegance, the original formulation of the KST is not necessarily well-suited for many practical applications, including PINNs Shukla et al. (); Yu et al. (). The limitations involved in implementing the univariate functions required by the theorem can be cumbersome, especially when dealing with real-world problems. As a result, alternative formulations and modifications of KST may be better equipped to address the challenges faced in modern approximation tasks, motivating ongoing research into more adaptable and scalable methods.

Alternative formulations of KST, which we will collectively refer to as Superposition Theorems, offer greater flexibility and more robust guarantees, making them better suited for modern deep learning applications, particularly in scientific computing. These variants relax some of the original theorem’s constraints, enabling more efficient computation and training of neural networks. By allowing for different types of function decompositions or less restrictive conditions, Superposition Theorems open up new possibilities for designing deep learning architectures that can improve performance and scalability in challenging settings like PINNs.

The summary of contributions in this paper are as follows:

- •
We argue for using alternative versions of the KST for building neural network architectures, as shown in table 1. Despite being the most popular iteration, Kolmogorov’s original formula is the least suited for practical implementations.
- •
We propose ActNet, a novel neural network architecture leveraging the representation theorem from Laczkovich () instead of the original KST employed in KANs. We prove ActNet’s universal approximation properties using fixed depth and width, and propose a well-balanced initialization scheme, which assures that the activations of the network scale well with the size of the network, regardless of depth or width. Furthermore, ActNet does not have vanishing derivatives, making it ideal for physics-informed applications.
- •
We evaluate ActNet in the context of PINNs, demonstrating its potential in a setting that uniquely challenges current deep learning practices. Our physics-informed experiments showcase ActNet’s ability to learn latent functions by minimizing PDE residuals, a task that aligns well with KST’s strengths. Although we do not claim state-of-the-art performance on all benchmarks, ActNet outperforms KANs in every single experiment we examined, and appears to be competitive against the current best MLP-based approaches.

Taken together, these contributions not only demonstrate the potential of ActNet as a novel architecture for tackling challenging problems in scientific computing and PDE simulation, but also open up a broader question: which formulation of the KST is most suitable for deep learning applications? While ActNet, leveraging Laczkovich’s theorem (Laczkovich, ), shows promising results, it represents just a first step in exploring this rich space of possibilities. By addressing key limitations of existing KST-based approaches while maintaining their theoretical strengths, our work with ActNet aims to pave the way for future research into optimal KST formulations for neural networks.

## 2 Superposition Theorems For Representing Complex Functions

Proposed in 1957 by Kolmogorov (), the Kolmogorov Superposition Theorem (KST, also known as Kolmogorov-Arnold representation theorem) offers a powerful theoretical framework for representing complex multivariate functions using sums and compositions of simpler functions of a single variable. This decomposition potentially simplifies the challenge of working with multivariate functions, breaking them down into manageable components.

In its original formulation, the KST states that any continuous function $f:[0,1]^{d}\to\mathbb{R}$ can be represented exactly as

$$ $\displaystyle f(x_{1},\dots,x_{d})=\sum_{q=0}^{2d}\Phi_{q}\left(\sum_{p=1}^{d} \phi_{q,p}(x_{p})\right),$ (1) $$

where $\Phi_{q}:\mathbb{R}\to\mathbb{R}$ and $\phi_{q,p}:[0,1]\to\mathbb{R}$ are univariate continuous functions known as outer functions and inner functions, respectively.

The first version of the KST was created to address Hilbert’s 13th problem Hilbert () regarding the representability of solutions of 7th-degree equations in terms of continuous functions of at most two variables. Since its inception in the mid 20th century, the theorem has since inspired attempts to apply it in fields like approximation theory and neural networks as a potential foundation for representing and approximating multivariate functions.

Over time, a multitude of different versions of KST have emerged, each with slightly different formulations and conditions, reflecting the evolving attempts to adapt the theorem to different applications and mathematical frameworks. Some of the most important variants are outlined in table 1. For more information, we recommend the excellent book written by Sprecher (), who has worked on KST related problems for the past half century.

**Table 1: Superposition formulas for a continuous function $f(x_{1},\dots,x_{d})$. Kolmorogov’s original result lead to many different versions of superpositions based on univariate functions. Although KANs employ the original formulation, other versions suggest distinct architectures. Other than the underlying formula, different theorems also use different assumptions about the domain of $f$ and regularity conditions of the inner/outter functions.**
| Version | Formula | Inner<br>Functions | Outer<br>Functions | Other<br>Parameters |
| --- | --- | --- | --- | --- |
| Kolmogorov<br>(1957) | $\sum_{q=0}^{2d}\Phi_{q}\left(\sum_{p=1}^{d}\phi_{q,p}(x_{p})\right)$ | $2d^{2}$ | $2d$ | N/A |
| Lorentz<br>(1962) | $\sum_{q=0}^{2d}g\left(\sum_{p=1}^{d}\lambda_{p}\phi_{q}(x_{p})\right)$ | $2d$ | $1$ | $\lambda\in\mathbb{R}^{d}$ |
| Sprecher<br>(1965) | $\sum_{q=0}^{2d}g_{q}\left(\sum_{p=1}^{d}\lambda_{p}\phi(x_{p}+qa)\right)$ | $1$ | $2d$ | $a\in\mathbb{R},\ \lambda\in\mathbb{R}^{d}$ |
| Laczkovich<br>(2021) | $\sum_{q=1}^{m}g\left(\sum_{p=1}^{d}\lambda_{pq}\phi_{q}(x_{p})\right)$ | $m$ | $1$ | $\lambda\in\mathbb{R}^{d\times m}$ |

##### KST And Neural Networks.

There is a substantial body of literature connecting the Kolmogorov Superposition Theorem (KST) to neural networks (Lin & Unbehauen, ; Köppen, ), notably pioneered by Hecht-Nielsen in the late 1980s (Hecht-Nielsen, ). KST was seen as promising for providing exact representations of functions, unlike approximation-based methods such as those of Cybenko (). More recently, Kolmogorov-Arnold Networks (KANs) (Liu et al., ) were proposed, attracting attention as a potential alternative to MLPs. Despite the initial enthusiasm, KANs have faced scrutiny due to mixed empirical performance and questions about experimental fairness (Yu et al., ; Shukla et al., ).

As detailed in table 1, the KST exists in multiple versions, each with unique conditions and formulations, offering significant flexibility over how functions can be decomposed and represented. These KST variants provide a theoretical foundation for designing new neural network architectures, rooted in well-established function approximation theory. Despite serving as the basis for KANs, the original formulation of the KST in equation (1) offers the least information about the structure of the superposition formula. As can be seen in table 1, Kolmogorov’s formulation is in fact, the only one with $\mathcal{O}(d^{2})$ unknown functions that have to be inferred, whereas other formulations scale linearly with the dimension $d$. This discrepancy motivates the use of other formulas from table 1 in order to design neural network architectures based on superposition theorems.

##### Kolmogorov’s Theorem Can Be Useful, Despite Its Limitations.

Previous efforts to apply KST in function approximation have had mixed results (Sprecher, ; ; Köppen, ; Kŭrková, ) as KST’s non-constructive nature and irregularity of the inner functions (e.g., non-differentiability) can pose significant challenges. Vitushkin () even proved that there exists certain analytical functions $f:[0,1]^{d}\to\mathbb{R}$ that are $r$ times continuously differentiable which cannot be represented by $r$ times continuously differentiable functions with less than $d$ variables, suggesting that in certain cases, the irregularity of inner functions is unavoidable.

The limitations of KST-based representations seem to be at direct odds with the many attempts over the past 50 years to use it as the basis for computational algorithms. The primary argument against the use of KST for computation is the provable pathological behavior of the inner functions, as detailed in Vitushkin (). However, this impediment may not be directly relevant for most applications for a few different reasons.

- 1.
Approximation vs. Exact Representation: The pathological behavior of inner functions mainly applies to exact representations, not approximations, which are more relevant in machine learning (Kŭrková, ).
- 2.
Function-Specific Inner Functions: Instead of using the same inner functions for all functions as suggested by most superposition theorems, allowing them to vary for specific target functions might address the irregularity problem in some cases.
- 3.
Deeper Compositions: Deep learning allows for more layers of composition, which could help alleviate the non-smoothness in approximation tasks

While further research is needed to verify these claims, they offer practical insights that motivate continued exploration, including the work presented in this paper.

## 3 ActNet - A Kolmogorov Inspired Architecture

As an area of research that recently regained widespread interest, there are many promising directions towards effectively integrating KST and deep learning, and many unanswered questions regarding the applicability and limitations of the theorem.

In light of these questions, we propose ActNet, a deep learning architecture inspired by the Laczkovich version of the KST presented in Laczkovich (). This formulation implies several differences from KANs and other similar methods. The basic building block of an ActNet is an ActLayer, whose formulation is derived from the KST and can additionally be thought of as a multi-head MLP layer, where each head has a trainable activation function (hence the name ActNet).

### 3.1 Theoretical Motivation

Theorem 1.2 in Laczkovich () presents the following version of Kolmogorov’s superposition theorem, which is the basis of the ActNet architecture:

###### Theorem 3.1 .

Let $C(\mathbb{R}^{d})$ denote the set of continuous functions from $\mathbb{R}^{d}\rightarrow\mathbb{R}$ and $m>(2+\sqrt{2})(2d-1)$ be an integer. There exists positive constants $\lambda_{ij}>0$, $j=1,\dots,d$; $i=1,\dots,m$ and $m$ continuous increasing functions $\phi_{i}\in C(\mathbb{R})$, $i=1,\dots,m$ such that for every bounded function $f\in C(\mathbb{R}^{d})$, there is a continuous function $g\in C(\mathbb{R})$ such that

$$ $\displaystyle f(x_{1},\dots,x_{d})$ $\displaystyle=\sum_{i=1}^{m}g\left(\sum_{j=1}^{d}\lambda_{ij}\phi_{i}(x_{j}) \right).$ $$

Using vector notation, this can equivalently be written as

$$ $\displaystyle f({\bm{x}})$ $\displaystyle=\sum_{i=1}^{m}g\left({\bm{\lambda}}_{i}\cdot\phi_{i}({\bm{x}}) \right),$ (2) $$

where the $\phi_{i}$ are applied element-wise and ${\bm{\lambda}}_{i}=(\lambda_{i1},\dots,\lambda_{id})\in\mathbb{R}^{d}$.

As a corollary, we have that this expressive power is maintained if we allow the entries of $\lambda$ to be any number in $\mathbb{R}$, instead of strictly positive, and do not restrict the inner functions to be increasing(^1^11An earlier version of ActNet attempted to impose these constraints to $\lambda$ and $\phi$, but this approach did not seem to lead to any practical gains, and potentially limited the expressive power of the network.). Furthermore, if we define the matrices ${\bm{\Phi}}(x)_{ij}:=\phi_{i}(x_{j})$ and ${\bm{\Lambda}}_{ij}:=\lambda_{ij}$ and let $S:\mathbb{R}^{a\times b}\rightarrow\mathbb{R}^{a}$ be the function that returns row sums of a matrix, or, if the input is a vector, the sum of its entries, then equation 2 can be equivalently expressed as

$$ $\displaystyle f(x)$ $\displaystyle=\sum_{i=1}^{m}g\left({\bm{\lambda}}_{i}\cdot\phi_{i}({\bm{x}})\right)$ (3) $\displaystyle=\sum_{i=1}^{m}g\left(\left[S({\bm{\Lambda}}\odot{\bm{\Phi}}(x)) \right]_{i}\right)$ (4) $\displaystyle=S\left(g\left[S({\bm{\Lambda}}\odot{\bm{\Phi}}(x))\right]\right),$ (5) $$

where $\odot$ is the Haddamard product between matrices.

From the formulations described in table 1, we believe this one from Laczkovich () is the best-suited for practical applications in deep learning. Some of the reasons for this are

- •
This version of the theorem is valid across all of $\mathbb{R}^{d}$, as opposed to being only applicable to the unit cube $[0,1]^{d}$, as is the case for most superposition theorems.
- •
It presents a flexible width size $m$, which yields exact representation as long as $m>(2+\sqrt{2})(2d-1)$, which scales linearly with the dimension $d$.
- •
It requires $\mathcal{O}(d)$ inner functions, as opposed to $\mathcal{O}(d^{2})$ in the case of KANs.
- •
This formulation has other attractive interpretations in the light of recent deep learning advances, as detailed in section 3.4.

###### Theorem 3.1 .

### 3.2 ActNet Formulation

The main idea behind the ActNet architecture is implementing and generalizing the structure of the inner functions implied by theorem 3.1. That is, we create a computational unit called an ActLayer, which implements the forward pass $S({\bm{\Lambda}}\odot{\bm{\Phi}}({\bm{x}}))$, and then compose several of these units sequentially. This is done by parametrizing the inner functions $\phi_{1},\dots,\phi_{m}$ as linear compositions of a set of basis functions.

Given a set of univariate basis functions $b_{1},\dots,b_{N}:\mathbb{R}\rightarrow\mathbb{R}$, let ${\bm{B}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{N\times d}$ return the basis expansion matrix defined by ${\bm{B}}(x)_{ij}=b_{i}(x_{j})$. The forward pass of the ActLayer component given an input ${\bm{x}}\in\mathbb{R}^{d}$ is then defined as

$$ $\displaystyle\text{ActLayer}_{\beta,{\bm{\Lambda}}}({\bm{x}})=S\left({\bm{ \Lambda}}\odot\beta{\bm{B}}({\bm{x}})\right),$ (6) $$

where the output is a vector in $\mathbb{R}^{m}$ and the trainable parameters are $\beta\in\mathbb{R}^{m\times N}$ and ${\bm{\Lambda}}\in\mathbb{R}^{m\times d}$. The graphical representation of this formula can be seen in Figure 1.

Under this formulation, we say that the matrix ${\bm{\Phi}}({\bm{x}})=\beta{\bm{B}}({\bm{x}})\in\mathbb{R}^{m\times d}$ where ${\bm{\Phi}}(x)_{ij}=\phi_{i}(x_{j})$ represent the inner function expansion of the layer and write

$$ ${\bm{\Phi}}({\bm{x}})=\left(\phi_{1}({\bm{x}}),\dots,\phi_{m}({\bm{x}})\right) =\beta{\bm{B}}({\bm{x}}),$ (7) $$

where $\phi_{k}:\mathbb{R}\to\mathbb{R}$ are real-valued functions that are applied element-wise, taking the form

$$ $\phi_{k}(t)=\sum_{j=1}^{N}\beta_{kj}b_{j}(t).$ (8) $$

Thus, the $k^{\text{th}}$ entry of an ActLayer output can be written as

$$ $\displaystyle\left(\text{ActLayer}_{\beta,{\bm{\Lambda}}}({\bm{x}})\right)_{k}$ $\displaystyle=\sum_{i=1}^{d}\lambda_{ki}\sum_{j=1}^{N}\beta_{kj}b_{j}(x_{i})$ (9) $\displaystyle=\sum_{i=1}^{d}\lambda_{ki}\phi_{k}(x_{i})=\lambda_{k,:}\cdot\phi _{k}({\bm{x}}),$ (10) $$

where $\lambda_{k,:}=(\lambda_{k1},\dots,\lambda_{kd})\in\mathbb{R}^{d}$ and $\phi_{k}({\bm{x}})=(\phi_{k}(x_{1}),\dots,\phi_{k}(x_{d}))\in\mathbb{R}^{d}$.

Additionally, before passing the input to ActLayers and before outputting the final result, ActNet uses linear projections from the latent dimension used for the ActLayers. We have empirically found this improved the performance of ActNet in some tasks by mixing the information of the different entries of ${\bm{x}}$. It is also possible to include an additive bias parameter at the end of the forward pass of each ActLayer. A visualization of the ActNet architecture can be seen in Figure 1.

Thus, the hyperparameters of an ActNet are: the embedding dimension $m\in\mathbb{N}$, which is used as the input and output dimension of the ActLayers; the number $N\in\mathbb{N}$ of basis functions used at each ActLayer; and the amount $L\in\mathbb{N}$ of ActLayers. A network of this size has $\mathcal{O}(Lm(m+N))$ trainable parameters. For comparison, a KAN with $L$ layers of width $m$, spline order $K$ and grid definition $G$ has $\mathcal{O}(Lm^{2}(G+K))$ trainable parameters, which grows faster than ActNet. On the other hand, an MLP of $L$ layers of width $m$ has $\mathcal{O}(Lm^{2})$ trainable parameters. Our experiments indicate that $N$ can often be set as low as $N=4$ basis functions, in which case the number of parameters of an ActNet is roughly similar to that of an MLP of the same width and depth, whereas KANs typically use $G>100$, implying a much faster rate of growth. A comparison of the number of trainable parameters and flops between ActNet, KAN and MLPs can be seen in Appendix B, and table 5.

Figure: Figure 1: (left) Visual representation of an individual ActLayer. The ActLayer architecture can be seen as a MultiHead MLP layer with tunable activations. (right) Visual representation of the ActNet architecture. The input vector ${\bm{x}}$ is first projected to an embedding dimension, then passed into $L$ composed blocks of ActLayer, and finally linearly projected into the desired output dimension.
Refer to caption: extracted/5896997/actnet_figures/ActLayer_diagram.png

Another advantage of the ActLayer over a KAN layer is that of simplicity in implementation. For example, using the einsum function found in packages like NumPy, PyTorch and JAX the forward pass of the ActLayer can even be implemented using a single line of code as out=einsum(’ij,jk,ik->k’,$B({\bm{x}})$,$\beta$,$\lambda$).

### 3.3 Universality

###### Definition 3.2 .

Given a set $\Omega\subseteq\mathbb{R}^{d}$, a neural network architecture is called a universal approximator on $\Omega$ if for every $f\in C(\Omega)$ and every $\varepsilon>0$ there exists a choice of hyperparameters and trainable parameters $\theta$ such that the function $f_{\theta}$ computed by this neural network satisfies

$$ $\displaystyle\sup_{x\in\Omega}|f(x)-f_{\theta}(x)|<\varepsilon.$ $$

With this definition, we now state the universality of ActNet, with a proof presented in Appendix E.

###### Theorem 3.3 .

A composition of two ActLayers is a universal approximator for any compact set $X\subseteq\mathbb{R}^{d}$. In other words, for any continuous function with compact domain $f:X\to\mathbb{R}$ and any $\epsilon>0$, there exists a composition of two ActLayers, where $|f(x)-Act(x)|<\epsilon$ for all $x\in X$.

In the construction implied by this proof, the first ActLayer has output dimension $m>(2+\sqrt{2})(2d-1)$, while the second ActLayer has output dimension 1. This means that we obtain universality using a width size that scales linearly with the input dimension $d$, potentially at the cost of requiring a large number $N$ of basis funcitons.

###### Definition 3.2 .

###### Theorem 3.3 .

### 3.4 Other Interpretations of the ActLayer

Although the representation formula from theorem 3.1 is the primary inspiration for the ActLayer, this unit can also be seen through other mathematically equivalent lenses. For example, the strategy of using several small units in parallel, followed by concatenating their outputs is exactly the idea behind multi-head attention employed in the transformer architecture from Vaswani (). In the same sense, the ActLayer can be seen as a Multi-Head MLP Layer where each head has a different, tunable, activation function $\phi_{i}$. In this analogy, each head first applies its individual activation function $\phi_{i}$, then applies a linear layer ${\bm{\lambda}}_{i}={\bm{\Lambda}}_{i,:}$ with output dimension 1. Finally, the output of each head is concatenated together before being passed to the next layer.

Alternatively, the ActLayer can also be seen as a lighter version of a KAN layer, where each inner function $\phi_{pq}$ has the form $\phi_{pq}(x)=\lambda_{pq}\phi_{q}(x)$. Assuming all inner functions are parametrized by the same set of basis functions, this is then equivalent to employing a KAN layer using Low-Rank Adaptation (LoRA) from Hu et al. (). This comparison however, overlooks the fact that the basis functions used for KAN are cubic B-splines, whereas ActNet employs a sinusoidal basis by default, as described in section 3.5.

All three interpretations are mathematically equivalent, and showcase potential benefits of the ActLayer through different perspectives.

### 3.5 Choice of Basis functions and Initialization

Although ActNet can be implemented with any choice of basis functions, empirically we have observed robust performance using basis functions $b_{1},\dots,b_{N}:\mathbb{R}\to\mathbb{R}$ of the form

$$ $b_{i}(t)=\frac{\sin(\omega_{i}t+p_{i})-\mu(\omega_{i},p_{i})}{\sigma(\omega_{i },p_{i})},$ (11) $$

where frequencies $w_{i}\sim N(0,1)$ are initialized randomly at each layer from a standard normal distribution and phases $p_{i}$ are initialized at 0. We also add a small $\epsilon>0$ to the denominator in (11) for numerical stability when $\sigma(\omega_{i},p_{i})$ is small. In the equation above, the constants $\mu(\omega_{i},p_{i})$ and $\omega(\sigma_{i},p_{i})$ are defined as the mean and standard deviation, respectively, of the variable $Y=\sin(\omega_{i}X+p_{i})$, when $X\sim N(0,1)$. These values are expressed in closed form as

$$ $\displaystyle\mu(\omega_{i},p_{i})$ $\displaystyle=\mathbb{E}[Y]=e^{\frac{-\omega^{2}}{2}}\sin(p_{i}),$ (12) $\displaystyle\sigma(\omega_{i},p_{i})$ $\displaystyle=\sqrt{\mathbb{E}[Y^{2}]-\mathbb{E}[Y]^{2}}=\sqrt{\frac{1}{2}- \frac{1}{2}e^{-2\omega_{i}^{2}}\cos(2p_{i})-\mu(\omega_{i},p_{i})^{2}}.$ (13) $$

This ensures that, if the inputs to the ActNet are normally distributed with mean 0 and variance 1 (a common assumption in deep learning), then the value of these basis functions will also have mean 0 and variance 1. After properly initializing the $\beta$ and $\lambda$ parameters (detailed in Appendix F), the central limit theorem then tells us that the output of an ActLayer will roughly follow a standard normal at initialization, which results in stable scaling of the depth and width of the network.

This can be formally stated as the theorem below, which implies that the activations of an ActNet will remain stable as depth and width increase. The proof can be seen in Appendix F.

###### Theorem 3.4 .

At initialization, if the input ${\bm{x}}\in\mathbb{R}^{d}$ is distributed as $N(0,{\bm{I}}_{d})$, then each entry of the output $ActLayer({\bm{x}})$ has mean 0 and variance 1. In the limit as either the basis size $N\to\infty$ or the input dimension $d\to\infty$ we get that each output is distributed as a standard normal.

###### Theorem 3.4 .

## 4 Experiments

We believe that the framework of Physics Informed neural networks (PINNs) (Raissi et al., ) provides an ideal testbed for evaluating ActNet and other KST-inspired architectures. PINNs align well with the strengths of the Kolmogorov Superposition Theorem in several key aspects:

- 1.
Low-dimensional domains: PINNs often deal with functions representing physical fields (e.g., velocity, temperature, pressure) that typically live in low-dimensional domains, where KST’s function approximation capabilities have been better studied.
- 2.
Derivative information: PINNs leverage derivative information through PDE residuals, a unique aspect that aligns with KST’s ability to represent functions and their derivatives. In particular, as detailed in Appendix C, the derivative of an ActNet is another ActNet, which tells us the higher order derivatives of the network will not vanish, a frequent problem with ReLU MLPs and many of its variants.
- 3.
Complex functions: PDEs often involve highly nonlinear and complex solutions, challenging traditional neural networks but potentially benefiting from KST-based approaches.
- 4.
Lack of direct supervision: Unlike traditional supervised learning, PINNs often learn latent functions without direct measurements, instead relying on minimizing PDE residuals. This indirect learning scenario presents a unique challenge that has posed difficulties to existing architectures (Wang et al., ) and may benefit from KST-inspired architectures.

By focusing our experiments on PINNs, we not only aim to demonstrate ActNet’s capabilities in a challenging and practically relevant setting, but also explore how KST-based approaches can contribute to scientific computing and PDE simulation. This choice of experimental framework allows us to evaluate ActNet’s performance in scenarios that closely align with the theoretical strengths of the KST, potentially revealing insights that could extend to broader applications in the future.

Motivated by these observations, several papers (Shukla et al., ; Wang et al., ; Howard et al., ; Shuai & Li, ; Koenig et al., ; Rigas et al., ; Patra et al., ) have been published attempting to specifically exploit KANs for physics informed training, with varying degrees of success. In the experiments that follow, we demonstrate that ActNet consistently outperforms KANs, both when using conducting our own experiments and when comparing against the published literature. Additionally, we also show that ActNet can be competitive against strong MLP baselines such as Siren Sitzmann et al. (), as well as state-of-the-art architectures for PINNs such as the modified MLP from Wang et al. ().

### 4.1 Ablations Studies

We compare the effectiveness of ActNet, KANs and MLPs in minimizing PDE residuals. In order to focus on PDE residual minimization, we enforce all initial and boundary conditions exactly, following Sukumar & Srivastava (), which removed the necessity for loss weight balancing schemes (Wang et al., ), and thereby simplifies the comparisons. For all experiments, we train networks of varying parameter counts(^2^22$\{10k,20k,40k,80k\}$ for Poisson and Helmholtz PDEs, and $\{25k,50k,100k\}$ for Allen-Cahn.), and for each parameter size, with 12 different hyperparemeter configurations for each architecture. Additionally, for the sake of robustness, for each hyperparameter configuration we run experiments using 3 different seeds and take the median result. This means that for each PDE, 144 experiments were carried out per architecture type. We do this as our best attempt at providing a careful and thorough comparison between the methods.

Additionally, for the Poisson and Helmholtz PDEs, we consider 6 different forcing terms, with frequency parameters $w\in\{1,2,4,8,16,32\}$. This means that for the Poisson PDE alone, 2,592 networks were trained independently, and similarly for Helmholtz. As the frequency parameter $w$ increases, solving the PDE becomes harder, as the model needs to adapt to very oscillatory PDE solutions. This is intended to showcase how ActNet’s expressive power and training stability directly translate to performance in approximating highly oscillatory functions.

Table 2 summarizes the results of these ablations, where we find that ActNet outperforms KANs across all benchmarks, often by several orders of magnitude, and also outperforms strong MLP baselines like Siren (Sitzmann et al., ) on most cases.

**Table 2: Best relative L2 errors and residual losses obtained by ablation experiments, as described in section 4.1. More details can be seen in Appendix G.**
|  | Relative L2 Error ($\downarrow$) | Residual Loss ($\downarrow$) |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Benchmark | ActNet | KAN | Siren | ActNet | KAN | Siren |
| Poisson ($w=16$) | 6.3e-4 | 1.8e-1 | 1.5e-1 | 4.6e-3 | 2.6e0 | 8.9e-1 |
| Poisson ($w=32$) | 6.3e-2 | 1.6e0 | 1.8e-1 | 8.0e-1 | 1.4e+3 | 8.4e+1 |
| Helmholtz ($w=16$) | 1.3e-3 | 2.8e-1 | 1.2e-1 | 5.3e-3 | 2.6e0 | 9.2e-1 |
| Helmholtz ($w=32$) | 1.2e-1 | 1.6e0 | 1.3e-1 | 5.7e-1 | 1.5e+3 | 8.1e+1 |
| Allen-Cahn | 6.8e-5 | 5.3e-4 | 2.1e-5 | 1.2e-8 | 3.8e-8 | 1.9e-8 |

##### Poisson Equation.

We consider the 2D Poisson PDE on $[-1,1]^{2}$ with Dirichlet boundary and forcing terms which yield exact solution $u(x,y)=\sin(\pi wx)\sin(\pi wy)$ for $w$ in $\{1,2,4,8,16,32\}$, which are shown in Appendix G.3 as Figure 5, along with a full description of the PDE. ActNet consistently outperforms the other two types of architectures for all frequencies $w$. In particular, the difference in performance between methods becomes starker as the frequency $w$ increases. For the highest two frequencies of $w=16$ and $w=32$, ActNet can outperform KAN and Siren by as much as two orders of magnitude. The full results comparing ActNet, KANs and Siren can be seen on the table in figure 7 of the Appendix.

##### Inhomogeneous Helmholtz Equation.

We consider the inhomogeneous Helmholtz PDE on $[-1,1]^{2}$ with Dirichlet boundary condition, which is fully described in Appendix G.4. Similarly to what was done for the Poisson problem, we set forcing terms so this PDE has exact solution $u(x,y)=\sin(wx)\sin(wy)$, which can be seen in Figure 5 for different values of $w$ and $\kappa=1$. Once again, ActNet outperforms both KANs and MLPs, across all frequencies (full details can be seen in figure 9 in the Appendix), with a starker discrepancy for high values of $w$.

Figure: Figure 2: Example predictions for the Helmholtz equation using $w=16$. The relative L2 errors for the ActNet, Siren and KAN solutions above are 1.04e-03, 8.82e-2 and 2.64e-1, respectively.
Refer to caption: extracted/5896997/experiment_figures/Helmholtz16_all_predictions_horizontal.png

##### Allen-Cahn Equation

We consider the Allen-Cahn PDE, which is fully described in Appendix G.5. Across all network sizes, ActNet outperforms KANs by around one order of magnitude. ActNet and Siren perform comparably, with ActNet yielding better results for smaller network sizes, and Siren performing better for larger network sizes, as can be seen in Figure 12 in the Appendix. However, for the larger networks, the relative L2 error is close to what can be achieved with single precision computation for neural networks, so it is possible that the discrepancy arises due to floating point error accumulation during training. Despite yielding slightly larger error in some cases, ActNet achieves the lowest final PDE residual loss across all network sizes when compared to both Siren and KAN, as can be seen in the Appendix, figure 13.

### 4.2 Comparisons against current State-of-the-Art

The goal of this section is to compare the performance of ActNet against some of the state-of-the-art results from the JaxPi library (Wang et al., ). As can be seen in table 3, ActNet is capable of improving the best results available in the literature for two very challenging PINN problems.

**Table 3: Relative L2 errors of ActNet versus state-of-the-art results from Wang et al. ().**
| Benchmark | ActNet | JaxPi |
| --- | --- | --- |
| Advection ($c=80$) | 9.50e-5 | 6.88e-4 |
| Kuramoto–Sivashinsky (first time window) | 1.34e-5 | 1.42e-4 |
| Kuramoto–Sivashinsky (full solution) | 8.53e-2 | 1.61e-1 |

##### Advection Equation.

We consider the 1D advection equation with periodic boundary conditions, which has been extensively studied in Wang et al. (); Daw et al. (); Krishnapriyan et al. () and is further described in Appendix G.6. Following Wang et al. (), we use initial conditions $g(x)=\sin(x)$ and the high transport velocity constant $c=80$. This yields a challenging problem for PINNs, with a highly oscillatory solution function. We compare an ActNet against the results presented in Wang et al. () and obtain a relative L2 error of $9.50\times 10^{-5}$, compared to their $6.88\times 10^{-4}$ result, see figure 3.

Figure: Figure 3: ActNet predictions for the advection equation ($c=80$). The relative L2 error is 9.50e-5, whereas the best result found in the literature is 6.88e-4 (Wang et al., ).
Refer to caption: x1.png

##### Kuramoto–Sivashinsky Equation.

The Kuramoto–Sivashinsky equation (fully described in Appendix G.7) is a chaotic fourth-order non-linear PDE that models laminar flames. Due to its chaotic nature, slight inaccuracies in solutions quickly lead to diverging solutions. Training PINNs for this PDE usually requires splitting the time domain into smaller windows, and sequentially solving one window at a time. Considering only the first window $t\in[0,0.1]$, ActNet obtains a relative L2 error an order of magnitude better than what is reported in Wang et al. (). Considering the full solution, ActNet also improves on the current best results from Wang et al. (), see figure 4.

Figure: Figure 4: ActNet predictions for the chaotic Kuramoto–Sivashinsky PDE. The relative L2 error is 8.53e-2, whereas the best result found in the literature is 1.61e-1 (Wang et al., ).
Refer to caption: x2.png

## 5 Discussion

Summary.
In this paper, we explore using alternative formulations of the Kolmogorov Superposition Theorem (KST) to develop neural network architectures and introduce ActNet. This architecture is inspired by Laczkovich’s version of KST, and is faster and simpler than Kolmogorov-Arnold Networks (KANs). We prove ActNet’s universality for approximating multivariate functions and propose an initialization scheme to maintain stable activations. In physics-informed experiments, ActNet consistently outperforms KANs and is competitive with leading MLP architectures, even surpassing state-of-the-art models like the modified MLP from Wang et al. () in some benchmarks.

Limitations.
The primary limitation of ActNet and other KST-based approaches to deep learning appears to be the slower computational speed when compared to plain MLPs, although ActNet mitigates this to a considerable degree when compared to KANs. Additionally, it might still be possible to optimize some aspects of our implementation, including potentially creating hardware-level implementations that are better suited for ActNet’s forward pass, as has been done with Transformers in recent years with FlashAttention (Dao et al., ).

Future Work.
Potential future directions in this field are determining the performance of ActNet and other KST-based architectures on data-driven problems, as well as studying the effect of replacing MLPs with these architectures inside larger neural networks such as transformers and U-Nets. We believe further study into KST-based architectures and specifically ActNet have the potential to advance not only scientific machine learning applications, but deep learning as a whole.

#### Ethics Statement

Our contributions enable advances in deep learning architectures. This has the potential to significantly impact a wide range of downstream applications. While we do not anticipate specific negative impacts from this work, as with any powerful predictive tool, there is potential for misuse. We encourage the research community to consider the ethical implications and potential dual-use scenarios when applying these technologies in sensitive domains.

#### Reproducibility Statement

The code used to carry out the ablation experiments is provided as a supplementary material for this submission. In the future, we plan on making our entire code-base publicly available on GitHub, along with integrating our ActNet implementation into the open-source JaxPi framework.

#### Author Contributions

L.F.G. conceived the methodology and conducted the experiments. P.P. provided funding and supervised this study. All authors helped conceptualize experiments and reviewed the manuscript.

#### Acknowledgments

We would like to acknowledge support from the US Department of Energy under the Advanced Scientific Computing Research program (grant DE-SC0024563) and the US National Science Foundation (NSF) Soft AE Research Traineeship (NRT) Program (NSF grant 2152205). We also thank Maryl Harris for helpful feedback when reviewing the writing of the manuscript and the developers of software that enabled this research, including JAX (Bradbury et al., ), Flax (Heek et al., ) Matplotlib (Hunter, ) and NumPy (Harris et al., ).

## Appendix A Mathematical Notation

Table 4 summarizes the symbols and notation used in this work.

**Table 4: Summary of the symbols and notation used in this paper.**
| Symbol | Meaning |
| --- | --- |
| ${\bm{x}}$ | A vector ${\bm{x}}=(x_{1},\dots,x_{d})$ in $\mathbb{R}^{d}$ |
| $C(\Omega)$ | The set of continuous functions from a set $\Omega$ to $\mathbb{R}$ |
| $g$ | Outer function implied by a superposition theorem (see table 1) |
| $\phi_{i}$ | Inner functions implied by a superposition theorem (see table 1) |
| ${\bm{\Phi}}({\bm{x}})$ | Inner function expansion matrix with $ij$ entry equal to $\phi_{i}(x_{j})\in\mathbb{R}$ |
| $b_{i}$ | Basis function from $\mathbb{R}$ to $\mathbb{R}$ |
| ${\bm{B}}({\bm{x}})$ | Basis expansion matrix with $ij$ entry equal to $b_{i}(x_{j})\in\mathbb{R}$ |
| ${\bm{\Lambda}}$ | Matrix with $ij$ entry equal to $\lambda_{ij}\in\mathbb{R}$ |
| $S$ | Function that returns row sums of a matrix, or the sum of entries of a vector |
| $\odot$ | Hadamard (element-wise) product between vectors/matrices |
| $\theta$ | Parameters of a neural network model |
| ${\bm{I}}_{d}$ | The identity matrix of size $d$ |
| $\mathbb{E}$ | Expectation operand |

## Appendix B Comparison on Parameters & Flops

From equation (9), each of the $k=1,\dots m$ outputs of an ActLayer with input dimension $d$, output dimension $m$ and $N$ basis functions is defined as

$$ $\displaystyle\left(\text{ActLayer}_{\beta,{\bm{\Lambda}}}({\bm{x}})\right)_{k}$ $\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{N}\lambda_{ki}\beta_{kj}b_{j}(x_{i}).$ $$

This sum entails $\mathcal{O}(dN)$ operations per output. Since there are $m$ outputs total, the computational complexity of this ActLayer is $\mathcal{O}(mdN)$. This means that this algorithmic complexity is a little slower than that of the traditional dense layer of an MLP (which is $\mathcal{O}(dm)$), but around the same order of magnitude for low values of $N$. In our experiments, we have observed that the basis size $N$ usually does not need to be large, with $N=4$ usually yielding the best results. In particular, in all experiments where ActNet beats state-of-the-art results, we do so by setting $N=4$ basis function, as can be seen in section 4.2. This setting makes an ActNet around $\approx$2-3 times slower than an MLP of the same size. To see this phenomena in practical terms, we direct the reader to table 9, where computational times for the Allen-Cahn PDE are compared for setting $N$ in $\{4,8,16\}$.

In the case of a KAN layer, the algorithmic complexity is $\mathcal{O}(mdK(G+K))$ Yu et al. (), where $G$ is the grid size hyper-parameter, and $K$ is the spline order (usually taken to be $K=3$). Unlike ActNet, however, typical values of $G$ can range from around 5 to several hundred, with the original KAN paper (Liu et al., ) using experiments with as high as $G=1000$, which makes the architecture impractical for realistic values of widths and depths. Once again, to see this phenomena in practical terms, we direct the reader to table 9, where we set $G=N$ in $\{4,8,16\}$. This comparison showcases another advantage of ActNets over KANs.

A summary of the hyperparameters of ActNet, KAN and MLP, along with their implied parameter count and flops of a forward pass can be seen in table 5.

**Table 5: Table comparing parameter and algorithmic complexity of ActNet, KAN and MLP. Since the typical value of $N$ is low, ActNet performs at slower, but comparable speed to MLPs, whereas KANs become very slow as the value of $G$ increases.**
| Architecture | Hyperparameters | Parameter Count | FLOPs |
| --- | --- | --- | --- |
| ActNet | depth $L$<br>hidden dim. $m$<br>basis size $N$ | $\mathcal{O}(Lm(m+N))$ | $\mathcal{O}(Lm^{2}N)$ |
| KAN | depth $L$<br>hidden dim. $m$<br>grid size $G$<br>spline order $K$ | $\mathcal{O}(Lm^{2}(G+K))$ | $\mathcal{O}(Lm^{2}K(G+K))$ |
| MLP | depth $L$<br>hidden dim. $m$ | $\mathcal{O}(Lm^{2})$ | $\mathcal{O}(Lm^{2})$ |

## Appendix C The Derivative of an ActNet is Another ActNet

From equation (9), we get that the partial derivative of the $k$th output of an ActLayer with respect to $x_{l}$ is

$$ $\displaystyle\frac{\partial}{\partial x_{l}}\left(\text{ActLayer}_{\beta,{\bm{ \Lambda}}}({\bm{x}})\right)_{k}$ $\displaystyle=\frac{\partial}{\partial x_{l}}\left(\sum_{i=1}^{d}\lambda_{ki} \sum_{j=1}^{N}\beta_{kj}b_{j}(x_{i})\right)$ $\displaystyle=\lambda_{kl}\sum_{j=1}^{N}\beta_{kj}\frac{\partial}{\partial x_{ l}}b_{j}(x_{l})$ $\displaystyle=\lambda_{kl}\sum_{j=1}^{N}\beta_{kj}b_{j}^{\prime}(x_{l})$ $\displaystyle=\lambda_{kl}\phi_{k}^{\prime}(x_{l}),$ $$

where $\phi_{k}^{\prime}(t)=\sum_{j=1}^{N}\beta_{kj}b_{j}^{\prime}(x_{l})$ is the derivative of the $k$th inner function $\phi_{k}(t)=\sum_{j=1}^{N}\beta_{kj}b_{j}(x_{l})$.

This formula then tells us that the Jacobian $J_{\text{ActLayer}_{\beta,{\bm{\Lambda}}}}({\bm{x}})$ of an ActLayer is

$$ $\displaystyle J_{\text{ActLayer}_{\beta,{\bm{\Lambda}}}}({\bm{x}})={\bm{ \Lambda}}\odot{\bm{\Phi}}^{\prime}({\bm{x}}),$ (14) $$

where ${\bm{\Phi}}^{\prime}({\bm{x}})\in\mathbb{R}^{m\times d}$ is the matrix defined by ${\bm{\Phi}}^{\prime}({\bm{x}})_{ij}=\phi_{i}^{\prime}(x_{j})$. This formulation is precisely what we compute in the forward pass $\text{ActLayer}_{\beta,{\bm{\Lambda}}}({\bm{x}})=S({\bm{\Lambda}}\odot{\bm{
\Phi}}({\bm{x}}))$, with the caveat that this Jacobian is now a matrix instead of a single vector.

If the basis functions $b_{1},\dots,b_{N}$ are picked so that their derivatives do not deteriorate, as is the case with the sinusoidal basis proposed in section 3.5, then the universality result of theorem 3.3 tells us that a sufficiently large ActNet will be able to approximate derivatives up to any desired precision $\epsilon>0$. This statement once again makes the case for the potential of ActNet for Physics Informed Machine Learning, where a network needs to learn to approximate derivatives of a function, instead of point values.

## Appendix D ActNet Implementation Details

Our implementation of ActNet is carried out in Python using JAX (Bradbury et al., ) and Flax (Heek et al., ), although the architecture is also easily compatible with other deep learning frameworks such as PyTorch (Ansel et al., ) and TensorFlow Abadi et al. (). The code is available in the supplementary material provided with this manuscript.

### D.1 Frequency Parameter ω 0 subscript 𝜔 0 \omega_{0} italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT

In a manner analogous to what is done in Siren (Sitzmann et al., ) and random Fourier features (Tancik et al., ), we have observed that ActNet benefits from having a $\omega_{0}>0$ parameter that multiplies the original input to the network. Using such a parameter better allows the architecture to create solutions that attain to a specific frequency content, using higher values of $\omega_{0}$ for approximating highly oscillatory functions and smaller values for smoother ones. This parameter can be trainable, but we generally found it best to set it to a fixed value, which may be problem dependent, much like the case with Siren and Fourier Features.

### D.2 Training Frequencies and Phases

Additionally, it is possible to allow the frequencies $\omega_{i}$ and phases $p_{i}$ of that basis functions $b_{i}$ from equation (11) to be trainable parameters themselves, making the basis functions at each layer customizable to the problem at hand. However, the parameters $\omega_{i}$ and $p_{i}$ of the basis functions appear to be very sensitive during training, and optimization via stochastic gradient-based methods sometimes lead to unstable behavior. This issue, however, can be mitigated by using Adaptive Gradient Clipping (AGC) proposed in Brock et al. (), where each parameter unit has a different clip depending on the magnitude of the parameter vector. We empirically observe that a hyperparameter of $1e-2$ works well for training basis functions of an ActNet.

### D.3 Initialization Scheme

The parameters of an ActLayer are initialized according to table 6. As will be shown later in the proof of theorem 3.4, using these initializations and applying the central limit theorem then tells us that the output of an ActLayer, as defined in (6), will be normally distributed with mean 0 and variance 1, as long as either $N$ or $d$ are sufficiently large. This same argument can then be applied inductively on the composed layers of the ActNet to get a statement on the stability of the statistics of activations throughout the entire network.

The $\beta$ and ${\bm{\Lambda}}$ parameters can be initialized using either a Gaussian or uniform distribution. Having said that, for the sake of consistency we use uniform distributions for initializing these parameters for all experiments in section 4.

**Table 6: Initialization schemes for an ActLayer with input dimension $d$, output dimension $m$ and $N$ basis functions. Each entry of a parameter is sampled in an i.i.d. manner.**
| Trainable<br>Parameter | Dimension | Mean<br>of Entries | Std. Dev.<br>of Entries | Distribution |
| --- | --- | --- | --- | --- |
| $\beta$ | $\mathbb{R}^{m\times N}$ | 0 | $\frac{1}{\sqrt{N}}$ | Gaussian/uniform |
| ${\bm{\Lambda}}$ | $\mathbb{R}^{m\times d}$ | 0 | $\frac{1}{\sqrt{d}}$ | Gaussian/uniform |
| $\omega$ | $\mathbb{R}^{N}$ | 0 | 1 | Gaussian |
| $p$ | $\mathbb{R}^{N}$ | 0 | 0 | constant |
| bias (optional) | $\mathbb{R}^{m}$ | 0 | 0 | constant |

## Appendix E Proof of theorem 3.3

###### Proof of theorem 3.3 .

By KST (theorem 3.1), we know that for any $m>(2+\sqrt{2})(2d-1)$, the function $f:\mathbb{R}^{d}\to\mathbb{R}$ has a representation of the form

$$ $\displaystyle f(x)=\sum_{q=1}^{m}g\left(\sum_{p=1}^{d}\lambda_{pq}\phi_{q}(x_{ p})\right),$ $$

where $g:\mathbb{R}\to\mathbb{R}$ and $\phi_{1},\dots,\phi_{m}:\mathbb{R}\to\mathbb{R}$ are continuous univariate functions and $\lambda_{pq}>0$ are positive scalars.

Now, let $Y:=\left\{\sum_{p=1}^{d}\lambda_{pq}\phi_{q}(x_{p})|x\in X,q=1,\dots,m\right\}$ be the set of possible inputs to $g$. Since $X$ is compact and each of the $\phi_{q}$ is continuous, this means that $Y$ must be compact as well.

Now let $\epsilon^{\prime}:=\epsilon/m$. Since $g$ is continuous and its inputs are confined to a compact domain $Y$, it must be uniformly continuous. That is, there exists some $\delta>0$ such that for all $y,z\in Y$, if $|z-y|<\delta$ then $|g(z)-g(y)|<\epsilon^{\prime}/2$. Additionally, due to Weierstrass theorem(^3^33Weierstrass theorem states that polynomials are dense in $C(X)$ for any compact $X\subset\mathbb{R}$.) we know that there exists a polynomial $poly_{g}:\mathbb{R}\to\mathbb{R}$ of degree $N_{g}$ such that for any $y\in Y$ we have that $|g(y)-poly_{g}(y)|<\epsilon^{\prime}/2$. Therefore, we get that

$$ $\displaystyle|g(z)-pol(y)|$ $\displaystyle=|g(z)-g(y)+g(y)-poly_{g}(y)|$ $\displaystyle\leq|g(z)-g(y)|+|g(y)-poly_{g}(y)|$ $\displaystyle<\epsilon^{\prime}/2+\epsilon^{\prime}/2=\epsilon^{\prime}.$ $$

Which implies that if $z_{1},\dots,z_{m}$ and $y_{1},\dots,y_{m}$ are scalars such that $|z_{q}-y_{q}|<\delta$ for all $q=1,\dots,m$, then we have

$$ $\displaystyle\left|\sum_{q=1}^{m}g(z_{q})-\sum_{q=1}^{m}poly_{g}(y_{q})\right|$ $\displaystyle\leq\sum_{q=1}^{m}\left|g(z_{q})-poly_{g}(y_{q})\right|$ $\displaystyle<\sum_{q=1}^{m}\epsilon^{\prime}=m\epsilon^{\prime}=\epsilon.$ $$

Therefore, by making sure that $\sum_{p=1}^{d}\lambda_{pq}\phi_{q}(x_{p})$ is at most $\delta$ far from the approximation $\sum_{p=1}^{d}\lambda_{pq}\hat{\phi}_{q}(x_{p})$ for all $q=1,\dots,m$, we can finish the proof of universality for ActNet.

To prove this final statement, we first set $M:=\max\{\lambda_{pq}|p=1,\dots,d;\ q=1,\dots,m\}$ and then note that by Weierstrass’ Theorem, for each $\phi_{q}$ there exists a polynomial $poly_{q}$ of degree $N_{q}$ such that $|\phi_{q}(z)-poly_{q}(z)|<\frac{\delta}{dM}$ for all $z\in Y$. Thus, we have that

$$ $\displaystyle\left|\sum_{p=1}^{d}\lambda_{pq}\phi_{q}(x_{p})-\sum_{p=1}^{d} \lambda_{pq}poly_{q}(x_{p})\right|$ $\displaystyle\leq\sum_{p=1}^{d}\lambda_{pq}\left|\phi_{q}(x_{p})-poly_{q}(x_{p })\right|$ $\displaystyle<\sum_{p=1}^{d}\lambda_{pq}\frac{\delta}{dM}$ $\displaystyle\leq\sum_{p=1}^{d}\frac{\delta}{d}$ $\displaystyle=\delta,$ $$

which completes the proof.
∎

###### Proof of theorem 3.3 .

## Appendix F Proof of Statistical Stability of Activations Through Layers

Before proving theorem 3.4, we prove the following lemma, where we compute the formulas stated in equations (12-13)

###### Lemma F.1 .

Let $\omega,p\in\mathbb{R}$ be fixed real numbers and $X\sim N(0,1)$ be a random variable distributed as a standard normal. Then the random variable $Y=\sin(\omega X+p)$ has mean and variance as follows:

$$ $\displaystyle\mathbb{E}[Y]$ $\displaystyle=e^{\frac{-\omega^{2}}{2}}\sin(p),$ (15) $\displaystyle Var[Y]$ $\displaystyle=\frac{1}{2}-\frac{1}{2}e^{-2\omega^{2}}\cos(2p)-e^{-\omega^{2}} \sin(p)^{2}.$ (16) $$

###### Proof.

First, we note that for any $a>0$ and $b\in\mathbb{C}$ we have that

$$ $\displaystyle\int_{-\infty}^{\infty}e^{-ax^{2}+bx}=\frac{\pi}{a}e^{\frac{b^{2} }{4a}}.$ (17) $$

Now, using this fact and the complex representation $\sin(\theta)=\frac{e^{i\theta}-e^{-i\theta}}{2i}$ we have that

$$ $\displaystyle\mathbb{E}[Y]$ $\displaystyle=\int_{-\infty}^{\infty}\sin(\omega x+p)\left(\frac{1}{\sqrt{2\pi }}e^{\frac{-x^{2}}{2}}\right)dx$ $\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\frac{e^{i(\omega x+ p)}-e^{-i(\omega x+p)}}{2i}e^{\frac{-x^{2}}{2}}dx$ $\displaystyle=\frac{1}{2i\sqrt{2\pi}}\int_{-\infty}^{\infty}\left(e^{i(\omega x +p)-\frac{x^{2}}{2}}-e^{-i(\omega x+p)-\frac{x^{2}}{2}}\right)dx$ $\displaystyle=\frac{1}{2i\sqrt{2\pi}}\left(\int_{-\infty}^{\infty}e^{-\frac{1} {2}x^{2}+i\omega x+ip}dx-\int_{-\infty}^{\infty}e^{-\frac{1}{2}x^{2}-i\omega x -ip}dx\right)$ $\displaystyle=\frac{1}{2i\sqrt{2\pi}}\left(e^{ip}\int_{-\infty}^{\infty}e^{- \frac{1}{2}x^{2}+i\omega x}dx-e^{-ip}\int_{-\infty}^{\infty}e^{-\frac{1}{2}x^{ 2}-i\omega x}dx\right)$ $\displaystyle=\frac{1}{2i\sqrt{2\pi}}\left(e^{ip}e^{\frac{-w^{2}}{2}}\sqrt{2 \pi}-e^{-ip}e^{-\frac{w^{2}}{2}}\sqrt{2\pi}\right)$ $\displaystyle=\frac{1}{2i}\left(e^{ip}e^{\frac{-w^{2}}{2}}-e^{-ip}e^{-\frac{w^ {2}}{2}}\right)$ $\displaystyle=e^{-\frac{w^{2}}{2}}\frac{e^{ip}-e^{-ip}}{2i}$ $\displaystyle=e^{-\frac{w^{2}}{2}}\sin(p),$ $$

which proves the claim for the expectation of $Y$. As for the variance, we first compute

$$ $\displaystyle\mathbb{E}[Y^{2}]$ $\displaystyle=\int_{-\infty}^{\infty}\sin(\omega x+p)^{2}\left(\frac{1}{\sqrt{ 2\pi}}e^{\frac{-x^{2}}{2}}\right)dx$ $\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{\infty}^{\infty}\frac{e^{2i(\omega x+ p)}+e^{-2i(\omega x+p)}-2}{-4}e^{\frac{-x^{2}}{2}}dx$ $\displaystyle=-\frac{1}{4\sqrt{2\pi}}\int_{-\infty}^{\infty}\left(e^{2i(\omega x +p)}e^{\frac{-x^{2}}{2}}+e^{-2i(\omega x+p)}e^{\frac{-x^{2}}{2}}-2e^{\frac{-x^ {2}}{2}}\right)dx$ $\displaystyle=-\frac{1}{4\sqrt{2\pi}}\left(\int_{-\infty}^{\infty}e^{2i(\omega x +p)}e^{\frac{-x^{2}}{2}}dx+\int_{-\infty}^{\infty}e^{-2i(\omega x+p)}e^{\frac{ -x^{2}}{2}}dx-\int_{-\infty}^{\infty}2e^{\frac{-x^{2}}{2}}dx\right)$ $\displaystyle=-\frac{1}{4\sqrt{2\pi}}\left(\int_{-\infty}^{\infty}e^{-\frac{1} {2}x^{2}+2i\omega x+2ip}dx+\int_{-\infty}^{\infty}e^{-\frac{1}{2}x^{2}-2i \omega x-2ip}dx-2\sqrt{2\pi}\right)$ $\displaystyle=-\frac{1}{4\sqrt{2\pi}}\left(e^{2ip}\int_{-\infty}^{\infty}e^{- \frac{1}{2}x^{2}+2i\omega x}dx+e^{-2ip}\int_{-\infty}^{\infty}e^{-\frac{1}{2}x ^{2}-2i\omega x}dx\right)+\frac{1}{2}$ $\displaystyle=-\frac{1}{4\sqrt{2\pi}}\left(e^{2ip}e^{\frac{-4\omega^{2}}{2}} \sqrt{2\pi}+e^{-2ip}e^{\frac{-4\omega^{2}}{2}}\sqrt{2\pi}\right)+\frac{1}{2}$ $\displaystyle=-e^{-2\omega^{2}}\left(\frac{e^{2ip}+e^{-2ip}}{4}\right)+\frac{1 }{2}$ $\displaystyle=-e^{-2\omega^{2}}\frac{\cos(2p)}{2}+\frac{1}{2}$ $\displaystyle=\frac{1}{2}-\frac{1}{2}e^{-2\omega^{2}}\cos(2p),$ $$

which means that

$$ $\displaystyle Var[Y]$ $\displaystyle=\mathbb{E}[Y^{2}]-\mathbb{E}[Y]^{2}$ $\displaystyle=\frac{1}{2}-\frac{1}{2}e^{-2\omega^{2}}\cos(2p)-e^{-w^{2}}\sin(p )^{2},$ $$

thus completing the proof.
∎

We are now ready to prove the main result of this section.

###### Proof of theorem 3.4 .

The $k$th output of an ActLayer is computed as

$$ $\left(\text{ActLayer}_{\beta,{\bm{\Lambda}}}({\bm{x}})\right)_{k}=\sum_{i=1}^{ d}\sum_{j=1}^{N}\lambda_{ki}\beta_{kj}\sin(\omega_{j}x_{i}+p_{j}).$ $$

Therefore, by linearity of expectation and then lemma F.1 we have

$$ $\displaystyle\mathbb{E}_{{\bm{x}}\sim N(0,I)}\left[\left(\text{ActLayer}_{ \beta,{\bm{\Lambda}}}({\bm{x}})\right)_{k}\right]$ $\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{N}\lambda_{ki}\beta_{kj}\mathbb{E}_{{ \bm{x}}\sim N(0,I)}\left[\sin(\omega_{j}x_{i}+p_{j})\right]$ $\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{N}\lambda_{ki}\beta_{kj}.$ $$

Since ${\bm{\Lambda}}$ and $\beta$ are independent and all entries have mean 0, this means that the $k$th output has mean 0. As for the variance $\sigma^{2}_{k}$, we now compute

$$ $\displaystyle\sigma^{2}_{k}$ $\displaystyle=\mathbb{E}\left[\left(\text{ActLayer}_{\beta,{\bm{\Lambda}}}({ \bm{x}})\right)_{k}^{2}\right]$ $\displaystyle=\mathbb{E}\left[\left(\sum_{i=1}^{d}\sum_{j=1}^{N}\lambda_{ki} \beta_{kj}\sin(\omega_{j}x_{i}+p_{j})\right)^{2}\right]$ $\displaystyle=\mathbb{E}\left[\left(\sum_{i=1}^{d}\sum_{j=1}^{N}\sum_{a=1}^{d} \sum_{b=1}^{N}\lambda_{ki}\beta_{kj}\sin(\omega_{j}x_{i}+p_{j})\lambda_{ka} \beta_{kb}\sin(\omega_{b}x_{a}+p_{b})\right)\right]$ $\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{N}\sum_{a=1}^{d}\sum_{b=1}^{N}\mathbb{ E}\left[\lambda_{ki}\beta_{kj}\sin(\omega_{j}x_{i}+p_{j})\lambda_{ka}\beta_{kb }\sin(\omega_{b}x_{a}+p_{b})\right]$ $\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{N}\mathbb{E}\left[\lambda_{ki}^{2} \beta_{kj}^{2}\sin(\omega_{j}x_{i}+p_{j})^{2}\right]$ $\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{N}\mathbb{E}\left[\lambda_{ki}^{2} \right]\cdot\mathbb{E}\left[\beta_{kj}^{2}\right]\cdot\mathbb{E}_{{\bm{x}}\sim N (0,I)}\left[\sin(\omega_{j}x_{i}+p_{j})^{2}\right]$ $\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{N}\frac{1}{d}\cdot\frac{1}{N}\cdot 1$ $\displaystyle=1.$ $$

This proves the first part of the theorem (outputs have mean 0 and variance 1 at initialization). Now to prove the second part (outputs are distributed normally as either $N$ or $d$ go to infinity), we note that from the deduction above, for any fixed $i=1,\dots,d$, the random variable $X_{i}^{d}\in\mathbb{R}$ defined by

$$ $X_{i}^{d}:=\sum_{j=1}^{N}\lambda_{ki}\beta_{kj}\sin(\omega_{j}x_{i}+p_{j}),$ $$

has mean 0 0 and variance $\frac{1}{d}$.

Therefore, by the central limit theorem, we have that

$$ $\left(\text{ActLayer}_{\beta,{\bm{\Lambda}}}({\bm{x}})\right)_{k}=\sum_{i=1}^{ d}X_{i}^{d},$ $$

will be distributed as a standard normal $N(0,1)$ as $d\rightarrow\infty$.
An analogous argument can be made, swapping the role of $d$ and $N$, which finishes the proof for the theorem.
∎

###### Lemma F.1 .

###### Proof.

###### Proof of theorem 3.4 .

## Appendix G Experimental Details

### G.1 KAN Implementation Details

All experiments were conducted in Python using JAX (Bradbury et al., ). Since the original implementation of KANs is not GPU-enabled and therefore very slow, we instead reference the “Efficient KAN” code found in [https://github.com/Blealtan/efficient-kan.git](https://github.com/Blealtan/efficient-kan.git). Since Efficient KANs were implemented originally in PyTorch (Ansel et al., ), we translated the code to Flax (Heek et al., ). This implementation can be found in the supplementary material submitted with this manuscript.

### G.2 Hyperparameter Ablation For Poisson, Helmholtz and Allen-Cahn PDEs

For every architecture, we disregard the parameters used in the first and last layers, as the input dimension of 2 and output dimension of 1 make them negligible compared to the number of parameters used in the intermediate layers. For all architectures and problems, we consider $\{1,2,4,6\}$ intermediate layers (not considering first and last layers), as well as three different architecture-specific hyperparameters, as detailed below. All model widths were then inferred to satisfy the model size at hand. This results in each architecture having 12 possible hyperparameter configurations for each parameter size (4 for depth, times 3 for architecture-specific hyperparameter).

Poisson & Helmholtz: For ActNet, we use $\{8,16,32\}$ values of $N$ and for KAN grid resolutions of $\{3,10,30\}$. For Siren we consider $\omega_{0}$ in $\{\frac{\pi w}{3},\pi w,3\pi w\}$.

Allen-Cahn: For ActNet, we use $\{4,8,16\}$ values of $N$ and for KAN grid resolutions of $\{4,8,16\}$. For Siren we consider $\omega_{0}$ in $\{10,30,90\}$.

### G.3 Poisson

We consider the 2D Poisson PDE with zero Dirichlet boundary condition, defined by

$$ $\displaystyle\Delta u(x,y)=f(x,y),\qquad$ $\displaystyle(x,y)\in[-1,1]^{2},$ $\displaystyle u(x,y)=0,\qquad$ $\displaystyle(x,y)\in\delta[-1,1]^{2}.$ $$

We use forcing terms of the form $f(x,y)=2\pi^{2}w^{2}\sin(\pi wx)\sin(\pi wy)$ for values of $w$ in $\{1,2,4,8,16\}$. Each of these PDEs has exact solution $u(x,y)=\sin(\pi wx)\sin(\pi wy)$, which are shown in Figure 5.

Each model was trained using Adam (Kingma & Ba, ) for 30,000 iterations, then fine tuned using LBFGS Liu & Nocedal () for 100 iterations. We use a batch size of 5,000 points uniformly sampled at random on $[0,1]^{2}$ at each step. For training using Adam, we use learning rate warmup from $10^{-7}$ to $5\times 10^{-3}$ over 1,000 iterations, then exponential decay with rate 0.75 every 1,000 steps, and adaptive gradient clipping with parameter $0.01$ as described in Brock et al. (). Furthermore, to avoid unfairness from different scaling between the residual loss and the boundary loss, we enforce the boundary conditions exactly for all problems by multiplying the output of the neural network by the factor $(1-x^{2})(1-y^{2})$, in a strategy similar to what is outlined in Sukumar & Srivastava ().

The final relative L2 errors can be seen in figure 7, while the final residual losses can be seen in figure 8. An example solution is plotted in figure 6 and sample computational times are reported in table 7.

Figure: Figure 5: Visualization of the different target functions for the Poisson and Helmholtz PDEs. As the frequency $w$ increases, the magnitude of gradients become larger and the laplacian $\Delta u(x,y)=2\pi^{2}w^{2}\sin(wx)\sin(wy)=2\pi^{2}w^{2}u(x,y)$ scales as $\propto w^{2}$, thus creating a challenging problem for PINNs.
Refer to caption: extracted/5896997/actnet_figures/poisson_target.png

Figure: Figure 6: Example predictions for the Poisson equation using $w=32$. The relative L2 errors for the ActNet, and Siren solutions above are 6.42e-02, 1.91e-1, respectively. Predictions for KANs are not plotted here, as they did not converge for this example.
Refer to caption: extracted/5896997/experiment_figures/Poisson32_predictions_horizontal.png

Figure: Figure 7: Results for the 2D Poisson problem using PINNs. For each hyperparameter configuration, 3 different seeds were used for initialization, and the median result was used. For each square, the best hyperparameter configuration (according to the median) is reported. The best performing method for each frequency $w$ and each number of parameters is underlined.
Refer to caption: extracted/5896997/experiment_figures/PoissonL2ErrorTable.png

Figure: Figure 8: Final residual loss for the Poisson problem using PINNs. For each hyperparameter configuration, 3 different seeds were used for initialization, and the median result was used. For each square, the best hyperparameter configuration (according to the median) is reported. The best performing method for each number of parameters is underlined.
Refer to caption: extracted/5896997/experiment_figures/PoissonLossTable.png

**Table 7: Average computational time per Adam training iteration for the Poisson problem on a Nvidia RTX A6000 GPU. All times are reported in milliseconds. Each cell averages across three seeds and four different widths and depths hyperparameters as reported on G.2.**
| Model<br>Size | ActNet<br>($N=8$) | ActNet<br>($N=16$) | ActNet<br>($N=32$) | KAN<br>($G=3$) | KAN<br>($G=10$) | KAN<br>($G=30$) | Siren |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 10k | 5.10 | 7.17 | 9.59 | 6.87 | 10.7 | 13.0 | 1.32 |
| 20k | 7.37 | 10.8 | 15.3 | 10.3 | 16.3 | 18.5 | 1.66 |
| 40k | 10.7 | 16.7 | 26.2 | 17.1 | 25.1 | 25.6 | 2.49 |
| 80k | 16.2 | 26.5 | 43.4 | 23.6 | 40.0 | 37.4 | 3.60 |

### G.4 Helmholtz

The inhomogeneous Helmholtz PDE with zero Dirichlet boundary condition is defined by

$$ $\displaystyle\Delta u(x,y)+\kappa^{2}u(x,y)=f(x,y,)\qquad$ $\displaystyle(x,y)\in[-1,1]^{2},$ $\displaystyle u(x,y)=0,\qquad$ $\displaystyle(x,y)\in\delta[-1,1]^{2}.$ $$

If we set $\kappa=1$ and the forcing term to be $f(x,y)=(\kappa-2\pi^{2}w^{2})\sin(\pi\omega x)\sin(\pi\omega y)$, then this PDE has exact solution $u(x,y)=\sin(wx)\sin(wy)$, which is pictured in figure 5.

Each model was trained using Adam (Kingma & Ba, ) for 30,000 iterations, then fine tuned using LBFGS Liu & Nocedal () for 100 iterations. We use a batch size of 5,000 points uniformly sampled at random on $[0,1]^{2}$ at each step. For training using Adam, we use learning rate warmup from $10^{-7}$ to $5\times 10^{-3}$ over 1,000 iterations, then exponential decay with rate 0.75 every 1,000 steps, and adaptive gradient clipping with parameter $0.01$ as described in Brock et al. (). Furthermore, to avoid unfairness from different scaling between the residual loss and the boundary loss, we enforce the boundary conditions exactly for all problems by multiplying the output of the neural network by the factor $(1-x^{2})(1-y^{2})$, in a strategy similar to what is outlined in Sukumar & Srivastava ().

The final relative L2 errors can be seen in figure 9, while the final residual losses can be seen in figure 10. An example solution is plotted in figure 2 and sample computational times are reported in table 8.

Figure: Figure 9: Results for the 2D Helmholtz problem using PINNs. For each hyperparameter configuration, 3 different seeds were used for initialization, and the median result was used. For each square, the best hyperparameter configuration (according to the median) is reported. The best performing method for each frequency $w$ and each number of parameters is underlined.
Refer to caption: extracted/5896997/experiment_figures/HelmholtzL2ErrorTable.png

**Table 8: Average computational time per Adam training iteration for the Helmholtz problem on a Nvidia RTX A6000 GPU. All times are reported in milliseconds. Each cell averages across three seeds and four different widths and depths hyperparameters as reported on G.2.**
| Model<br>Size | ActNet<br>($N=8$) | ActNet<br>($N=16$) | ActNet<br>($N=32$) | KAN<br>($G=3$) | KAN<br>($G=10$) | KAN<br>($G=30$) | Siren |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 10k | 5.03 | 7.12 | 9.53 | 6.84 | 10.7 | 13.0 | 1.31 |
| 20k | 7.33 | 10.8 | 15.2 | 10.2 | 16.2 | 18.5 | 1.67 |
| 40k | 10.6 | 16.6 | 26.0 | 17.1 | 25.0 | 25.7 | 2.50 |
| 80k | 16.2 | 26.3 | 43.1 | 23.6 | 40.0 | 37.4 | 3.57 |

Figure: Figure 10: Final residual loss for the Helmholtz problem using PINNs. For each hyperparameter configuration, 3 different seeds were used for initialization, and the median result was used. For each square, the best hyperparameter configuration (according to the median) is reported. The best performing method for each number of parameters is underlined.
Refer to caption: extracted/5896997/experiment_figures/HelmholtzLossTable.png

### G.5 Allen-Cahn

We consider the Allen-Cahn equation with Dirichlet boundary and initial condition $u(x,0)=x^{2}\cos(\pi x)$, defined as

$$ $\displaystyle\frac{\partial u}{\partial t}(x,t)-D\frac{\partial^{2}u}{\partial x ^{2}}(x,t)+5(u^{3}(x,t)-u(x,t)),\qquad$ $\displaystyle(x,t)\in[-1,1]\times[0,1],$ $\displaystyle u(x,0)=x^{2}\cos(\pi x),\qquad$ $\displaystyle x\in[-1,1],$ $\displaystyle u(-1,t)=u(1,t)=-1,\qquad$ $\displaystyle t\in[0,1].$ $$

Each model was trained using Adam (Kingma & Ba, ) for 100,000 iterations and batch size of 10,000 points uniformly sampled at random each step on $[-1,1]\times[0,1]$. We use learning rate warmup from $10^{-7}$ to $5\times 10^{-3}$ over 1,000 iterations, then exponential decay with rate 0.9 every 1,000 steps, stopping the decay once the learning rate hits $5\times 10^{-6}$. We also use adaptive gradient clipping with parameter $0.01$ as described in Brock et al. (). Furthermore, to avoid unfairness from different scaling between the residual loss and the boundary loss, we enforce the initial and boundary conditions exactly by setting the output of the model to $(1-t)(x^{2}\cos(\pi x))+t[(1-x^{2})u_{\theta}(x,t)-1]$, where $u_{\theta}$ is the output of the neural network, in a strategy similar to what is outlined in Sukumar & Srivastava ().

Finally, we also employ the causal learning strategy from Wang et al. (), updating the causal parameter every 10,000 steps according to the schedule [$10^{-1}$, $10^{0}$, $10^{1}$, $10^{2}$, $10^{3}$, $10^{3}$, $10^{3}$, $10^{4}$, $10^{4}$, $10^{4}$].

The final relative L2 errors can be seen in figure 12, while the final residual losses can be seen in figure 13. An example solution is plotted in figure 11 and sample computational times are reported in table 9.

**Table 9: Average computational time per Adam training iteration for the Allen-Cahn problem on a Nvidia RTX A6000 GPU. All times are reported in milliseconds. Each cell averages across three seeds and four different widths and depths hyperparameters as reported on G.2.**
| Model<br>Size | ActNet<br>($N=4$) | ActNet<br>($N=8$) | ActNet<br>($N=16$) | KAN<br>($G=4$) | KAN<br>($G=8$) | KAN<br>($G=16$) | Siren |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 25k | 4.16 | 6.29 | 10.8 | 9.01 | 24.5 | 56.5 | 2.04 |
| 50k | 5.94 | 10.2 | 18.2 | 12.3 | 35.5 | 84.1 | 2.76 |
| 100k | 9.65 | 16.9 | 30.1 | 17.1 | 48.8 | 122.9 | 4.30 |

Figure: Figure 11: ActNet predictions for the Allen-Cahn equation. The relative L2 error is 4.51e-05.
Refer to caption: extracted/5896997/experiment_figures/ActNet_AllenCahn_prediction.png

Figure: Figure 12: Results for the Allen-Cahn problem using PINNs. For each hyperparameter configuration, 3 different seeds were used for initialization, and the median result was used. For each square, the best hyperparameter configuration (according to the median) is reported. The best performing method for each number of parameters is underlined.
Refer to caption: extracted/5896997/experiment_figures/AllenCahnL2ErrorTable.png

Figure: Figure 13: Final residual loss for the Allen-Cahn problem using PINNs. For each hyperparameter configuration, 3 different seeds were used for initialization, and the median result was used. For each square, the best hyperparameter configuration (according to the median) is reported. The best performing method for each number of parameters is underlined.
Refer to caption: extracted/5896997/experiment_figures/AllenCahnFinalLossTable.png

### G.6 Advection

The 1D advection equation is defined by

$$ $\displaystyle\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0,\qquad$ $\displaystyle t\in[0,1],\ x\in(0,2\pi),$ $\displaystyle u(0,x)=g(x),\qquad$ $\displaystyle x\in(0,2\pi),$ $$

with initial condition $g(x)=\sin(x)$ and transport velocity constant $c=80$.

For our implementation, we use the JaxPi library from Wang et al. (; ). We copy over the configuration of their sota.py config file, along with their learnable periodic embedding, changing the architecture to an ActNet with 5 layers, embedding dimension of 256 and $\omega_{0}$ value of 5. We train using Adam for 300,00 iterations with adaptive gradient clipping (Brock et al., ) with parameter $0.01$. We use a starting learning rate of $10^{-3}$ and exponential rate decay of 0.9 every 5,000 steps. Each step took on average 40.0ms on a Nvidia RTX A6000 GPU.

### G.7 Kuramoto–Sivashinsky

The Kuramoto–Sivashinsky equation we considered is defined by

$$ $\displaystyle\frac{\partial u}{\partial t}+\alpha u\frac{\partial u}{\partial x }+\beta\frac{\partial^{2}u}{\partial x^{2}}+\gamma\frac{\partial^{4}u}{ \partial x^{4}}=0$ $\displaystyle\quad t\in[0,1],x\in[0,2\pi]$ $\displaystyle u(0,x)=u_{0}(x)$ $\displaystyle\quad x\in[0,2\pi]$ $$

with periodic boundary conditions and constants $\alpha=100/16$, $\beta=100/16^{2}$ and $\gamma=100/16^{4}$, along with initial condition $u_{0}(x)=\cos(x)(1+\sin(x))$.

For our implementation, we use the JaxPi library from Wang et al. (; ). We copy over the configuration of their sota.py config file, changing the architecture to an ActNet with 5 layers, embedding dimension of 256 and $\omega_{0}$ value of 5. We train using Adam (this time without adaptive gradient clipping) using a starting learning rate of $10^{-3}$ and exponential rate decay of 0.95 every 2,000 steps.

We split the time domain $[0,1]$ into 10 windows of equal length. Since the precision of the solution matters most in the initial time steps (most of the error from the latter windows is due to propagated error in the initial conditions), we train the first window for 250k steps, the second and third for 200k steps, and all other for 150k steps. This results in 15% less steps than what was done in JaxPi, where each window is trained for 200k steps. Each step took on average 66.7ms on a Nvidia RTX A6000 GPU.

## Peer Reviews

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: I would like to see, even in the annex, the behavior against Neural Operators.

is there other choices of the b(t) functions that works well?

How do you choose the hyper-parameters (N,m)?
soundness: 3
strengths: The paper is clear and justify its contribution, and put in context with the current state of art

The paper justifies the change in the KAN architecture and its relationship. The connection with the multi-head transformer is a bit stretched tho. 

The paper provides some experiments to show the potential of changing the representation of the KST.
summary: The paper presents an alternative formulation for a trainable network based on the KST. 

The paper presents the theoretical property of the model and experimentally validates the performance as an approximation function of PINN, on 3 PDEs.
weaknesses: The exposition is very good, the experiments show the advantage with respect to other KST architecture.

There is only the evaluation in the PINN context. It would be nice to have more experiments, for example against some neural operators and training from data.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. This paper presents a novel approach for solving PDEs, similar to PINN. However, since the PDEs are known, why not use traditional numerical methods, such as FEM or FVM, on GPUs? A fair comparison is needed, focusing on accuracy, efficiency, and implementation complexity.
2. How does the proposed method address high-dimensional input problems?
3. Line 161: Why is the pathological behavior of the inner function relevant to the exact representation rather than the approximation?
4. Table 2: A relative L2 error should suffice for accuracy comparison, so why include the residual loss? Additionally, why does the Allen-Cahn example in Table 2 lack a consistently best-performing model?
soundness: 3
strengths: 1. The authors provide a solid theoretical foundation for ActNet, detailing its universal approximation capabilities and presenting a well-motivated formulation based on KST. This paper compares the complexity of different formulations well.
2. Three 1D examples and two 2D examples are included to demonstrate that ActNet can achieve better performance. 
3. This paper compares the model performance with other open-source JaxPi frameworks to enhance reproducibility and fairness.
summary: Based on a new formulation of the Kolmogorov Superposition Theorem (KST), this paper introduces a neural network ActNet to solve PDEs. ActNet is proposed as a global approximation in solving PDEs, which is similar to PINN. The authors demonstrate ActNet's theoretical properties, including universality, and validate it on a range of PDE benchmarks, showing improved performance over KAN and SIREN.
weaknesses: 1. The selected examples primarily use sinusoidal forcing terms, including equations like Poisson, Helmholtz, and Allen-Cahn. The proposed model uses sine functions as basis functions, which justifies its improved performance over the spline basis functions used in KAN. In Figure 12, we can see that SIREN can have the best performance. Therefore, comprehensive benchmarks should be included for challenging 2D and 3D problems, including the Navier–Stokes equations and turbulence cases.
2. As this paper emphasizes the comparison between MLP and KAN, it is essential to include a general MLP model for performance evaluation alongside SIREN.

---

confidence: 5
contribution: 4
initial_rating: 6
presentation: 4
questions: 1. As you mentioned about JAXPI, Why you didn't compare with their latest result which is Piratenet [1] that has been pubilshed Feb 2024 rather than Causal PINN that from 2022. Also I saw that you have cited this paper as well.
2. What happen if you apply ActNet into more chaotic system, for example Navier–Stokes equations. is ActNet still performs better?
3.

[1] PirateNets: Physics-informed Deep Learning with Residual Adaptive Networks
soundness: 4
strengths: 1. **Originality**: The paper brings a fresh perspective to neural network design by leveraging alternative formulations of the Kolmogorov Superposition Theorem (KST). Introducing ActNet, based on Laczkovich's theorem, reflects a novel approach to overcoming the limitations of Kolmogorov-Arnold Networks (KANs), making KST more applicable to practical deep learning tasks, particularly within Physics-Informed Neural Networks (PINNs).

2. **Quality**: The research is thorough, with ActNet being tested across multiple benchmarks against established models like KANs and MLPs, demonstrating its advantages in function approximation and handling PDE simulations. Theoretical foundations are strong, with proofs of ActNet’s universal approximation properties, and empirical results consistently show ActNet’s improved performance in accuracy and stability.

3. **Clarity**: The paper is well-structured, presenting a clear motivation for the need for alternative KST formulations, followed by detailed descriptions of ActNet’s architecture and its theoretical underpinnings. The explanations of mathematical concepts and the positioning of ActNet within current scientific computing challenges are accessible and well-supported with illustrative figures and tables.

4. **Significance**: ActNet addresses critical limitations in applying KST to neural networks, opening new possibilities for KST-based models in scientific computing and PDE simulation. The model's competitive performance against leading MLP-based approaches highlights its potential impact on advancing scientific machine learning, particularly in low-dimensional function approximation and complex simulations, where existing architectures struggle.
summary: This paper investigates alternative approaches to the Kolmogorov Superposition Theorem (KST) for neural network design, addressing practical limitations of the original KST formulation, such as its complexity and lack of structural insights. Kolmogorov-Arnold Networks (KANs) utilize KST for function approximation but have shown inconsistent performance compared to traditional multilayer perceptrons (MLPs). The authors introduce ActNet, a scalable model that builds on KST while mitigating some of its original limitations. Evaluated within the Physics-Informed Neural Networks (PINNs) framework for PDE simulation, ActNet consistently outperforms KANs and competes well with leading MLP-based methods, marking it as a promising direction for KST-based deep learning in scientific computing.
weaknesses: 1.  Although ActNet performs well on PINNs, its comparisons are limited to specific benchmarks and do not consistently compare against the latest models for PINNs. 
2. The paper lacks detailed ablation studies on critical design choices within ActNet, such as the basis functions used or the impact of ActLayer depth. Including these analyses would clarify the sensitivity of ActNet’s performance to these parameters and help guide future implementations or adaptations of the model.
3. While the paper offers rigorous theoretical grounding, some sections—particularly those detailing the inner workings of the ActLayer—may be dense for readers unfamiliar with KST.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: * The original KAN paper boasts interpretability of the learnt network. Does your method offer any interesting results in terms of interpretability?
* Do you have any experiments out-of-distribution (OOD) regime? It would be interesting to see how the model reacts to change in initial conditions or boundary conditions for example.
* Why didn't you test on the function approximation tasks as presented in the original KAN paper? Given your method is a direct competitor to the KAN architecture, such results would be interesting to have.
* How does your method perform given a different function basis?
* Are $p$ and $\omega$ trainable parameters?
* My understanding of the proof of Theorem 3.3 is that you use Weierstrass theorem to create polynomial bases which approximate both i) the sum of $g$s ii) the weighted sum of $\phi$s up to the desired precision $\epsilon$. However, the link between this and the architecture is not immediately clear to me.

    From what I gather, each polynomial approximation can be represented by an ActLayer with a polynomial basis of size $N_g$/$N_q$ and well chosen $\beta$ values. Is my reasoning correct here? I think it would improve clarity if you mentioned how your proof relates to the proposed architecture.
soundness: 4
strengths: 1. I find the paper to be well-written. The exposition of KST and of the ActNet architecture are clearly presented. The contributions made by this work is also clearly outlined.
2. As far as I am aware, the technical contributions in this paper are novel. The proposed architecture outperforms other similar models (KST variants such as KANs and MLPs) on the chosen tasks, has better parameter efficiency and seems to benefit from good theoretical guarantees (universal approximation properties for a two layer network).
3. I found the experimental section of the paper to be strong. Information about implementation and the considered PDEs is clearly given either in the main paper or in the appendix. Multiple ablations are conducted to validate the model. The choice of experimental framework is well motivated and the code is publicly available on Github.

All in all, I find this paper presents an interesting alternative to KANs.
summary: The paper proposes ActNet, a novel neural network architecture based on Kolmogorov Superposition Theorem (KST), as an alternative to Kolmogorov-Arnold Networks (KANs). The authors use the Laczkovich version of KST to develop their architecture instead of the original KST emplyoed in KANs. They prove universal approximation properties and propose an initialization scheme for their method which they assures stability of network activations. 

Empirically, they evaluate their model on PDE simulation tasks using a PINN framework. They compare their model to KANs and MLPs on a range of PDE equations such as Poisson, Helmholtz, Allen-Cahn and Kuramoto-Sivashinsky. Their method has comparable accuracy to the benchmarks chosen for comparison.
weaknesses: 1. Although interesting when compared to KANs and other MLPs, this method is not competitive when compared to SoTA ML informed solvers which use more sophisticated inductive biases such as [1,2,3,4] for example.
2. I find the theoretical contribution of the paper to be a bit weak. Although an interesting result, it seems trivial to me that a method based on Superposition Theorems for representing multivariate functions would benefit from universal approximation properties. Moreover, it is not clear to me why the property presented in Theorem 3.4 is a selling point of your method. Even if this property didn't hold, we could simply use some form of LayerNorm to enforce this. Also, in the context of modelling physical equations, it is not clear to me that inputs would be distributed following a $(0,1)$ Normal distribution in the first place.
3. There are a few typos/minor mistakes in the proofs. I will list them by point below.

Comments about the proof of Theorem 3.3:
* Line 868: typo, I believe you meant $poly_g$ not $pol$.
* "by making $\sum \lambda \phi (x)$ sure that is at most $δ$ far from the approximation $\sum \lambda \hat{\phi} (x)$". You do not introduce/define this approximation before mentioning it here, I think it would make this step clearer if you defined what this approximation is.

Comments about the proof of Theorem 3.4
* "After properly initializing the $β$ and $λ$ parameters (detailed in Appendix F)" this is not detailed in the appendix. Thus, when defining the expectation and second moment of these variables, it is not clear where the values you obtain come from (around lines 981 and 1006 )
* Do you consider the that $p=0$ and $\omega \sim \mathcal{N}(0,1)$ for this proof? If so, the eq. at line 984 is not clear to me. Same for the step between lines 1006 and 1009.

Also here are some minor comments you may want to address for the final version: 
* "Table 1" not "table 1" (line 67 for instance).
* $m$ in Table 1 is not described until section 3. It would be nice to know (at least in big O) what this value is in terms of $d$ when looking at this table.
* $\beta$ should be bolded since it is a matrix, not a scalar (line 237).
* Inconsistent use of $\epsilon$ and $\varepsilon$ (around line 295 theorem 3.3 def 3.2).
* The grid size considered for the experiments is not reported in main paper and generally hard to find.
* The "Discussion" section appears to be more of a conclusion to me than a discussion of results.

[1] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., & Anandkumar, A. (2020). Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895.

[2] Gin, C., Lusch, B., Brunton, S. L., & Kutz, J. N. (2021). Deep learning models for global coordinate transformations that linearise PDEs. European Journal of Applied Mathematics, 32(3), 515-539.

[3] Lusch, B., Kutz, J. N., & Brunton, S. L. (2018). Deep learning for universal linear embeddings of nonlinear dynamics. Nature communications, 9(1), 4950.

[4] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Stuart, A., Bhattacharya, K., & Anandkumar, A. (2020). Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, 33, 6755-6766.

Meta-review: Based on a variant of Kolmogorov superposition theorem given by Miklos Laczkovich in 2021, the authors propose to use it as a deep neural network design. It used d inner functions, and one outer function. The ActNet based on this has d adaptive activation functions, different for each neuron. The numerical experiments are quite sparse in this case: 5 problems vs 4 methods However, the experiments themselves seem to be done correctly, for different number of parameters/seeds. 
Overall, although: 
1) No theoretical results (i.e. are the 'true' functions smooth, so they can be approximated well with simple basis functions?)
2) Numerical experiments are not 100% convincing that this approach is better

all the reviewers are quite positive in the end.

Decision: Accept (Spotlight)