Neural Networks are Integrable
In this study, we explore the integration of Neural Networks, a powerful
class of functions known for their exceptional approximation capabilities. Our
primary emphasis is on the integration of multi-layer Neural Networks, a
challenging task within this domain. To tackle this challenge, we introduce a
novel numerical method that consist of a forward algorithm and a corrective
procedure. Our experimental results demonstrate the accuracy achieved through
our integration approach.
Yucong Liu
# Neural Networks are Integrable

###### Abstract

In this study, we explore the integration of neural networks, a potent class of functions known for their exceptional approximation capabilities. Our primary emphasis is on the integration of multi-layer neural networks, a challenging task within this domain. To tackle this challenge, we introduce a novel numerical method that combines the forward algorithm with a corrective procedure. Our experimental results demonstrate the accuracy achieved through our integration approach.

keywords: Neural Network, Integration, Numerical Algorithm +
Footnote †: journal: Elsevier

## 1 Introduction

Deep learning models have demonstrated incredible power in various fields such as image and speech recognition, natural language processing, and autonomous driving in recent years. This work primarily centers on the deep neural network, also known as multi-layer perceptron or feed-forward neural network.

We give the definition of neural networks as below. A k-layer neural network from \(\mathbb{R}^{n_{0}}\) to \(\mathbb{R}^{n_{k}}\) is a layer-wise structure. The \(i\)-the layer, for \(i\in\{1,\cdots,k\}\), is defined as :

\[y^{(i)} =W^{(i)}x^{(i-1)}+b^{(i)},\] \[x^{(i)} =\sigma(y^{(i)}).\]

In each layer, \(W^{(i)}\in\mathbb{R}^{n_{i}\times n_{i-1}}\) is the weight matrix and \(b^{(i)}\in\mathbb{R}^{n_{i}}\) is the bias vector. The input of \(i\)-th layer \(x^{(i)}\in\mathbb{R}^{n_{i}}\) is the output of \((i-1)\)-th layer. Activation function \(\sigma\) is a nonlinear point-wise function, i.e \((\sigma(y^{(i)}))_{j}=\sigma(y^{(i)}_{j})\)
One of the most well-known activation function is rectified linear unit (ReLU), which is defined as \(\text{ReLU}(x)=\max(x,0)\).

There has been numerous works on theory, algorithms and methods for neural networks. One of the key factors contributing to the success of neural networks is the Universal Approximation Theorem (UAT), which states that these networks have the ability to approximate arbitrary continuous function.

Throughout the past century, several researchers such as Cybenko [3], Hornik [8], Pinkus [14], et al, directed their attention toward assessing the approximation capabilities of single-layer neural networks. However, recent work by Kidger and Lyons [9] has extended these findings. They demonstrated that multi-layer networks, each layer having a fixed width, possess the capacity to approximate arbitrary continuous functions within a compact set.

This work focuses on the topic of integration for both shallow and deep neural networks. While it is known that a neural network with a continuous activation function is integrable, the primary research emphasis has been on exploring the gradients of neural networks. The gradient's pivotal role in optimizing diverse loss functions has spurred extensive research efforts. In contrast, the integrability of neural networks, despite its significance in comprehending the models' theoretical properties, has garnered comparatively less attention. In this work, we aim to bridge this gap by providing explicit forms of integration for one-layer neural networks with any integrable activation function and deriving a piece-wise structure of the integration for multi-layer neural networks with ReLU activation function, along with a proposed algorithm with a corrector.

We will introduce our basic motivation in Section 2. For the integration of ReLU neural networks, we separate it into two cases, one-layer case in Section 3.1 and multi-layer case in Section 3.2. Our algorithm also works for Convolutional Neural Networks [10] and Residual Neural Networks [6], which will be introduced in Section 3.3. We will present our numerical experiments in Section 4 and discuss about future work in Section 5.

## 2 Basic Motivation

Let us begin by focusing on a classical numerical question: how can we approximate the integral of a function given only some samples? Suppose we have a integrable function \(f:[a,b]\rightarrow\mathbb{R}\). Let \(F(x)=\int_{a}^{x}f(t)dt\) denotes the integral of function \(f\), and suppose we have \(N\) samples of \(f\), i.e., \(\{(x_{n},f(x_{n}))\}_{n=1}^{N}\). How can we recover \(F\) with these samples without knowing
the corresponding values \(F(x_{n})\) or the formula of \(F\)? Especially, can we get an accurate estimation for \(F(b)=\int_{a}^{b}f(t)dt\)?

Different from numerical integration algorithms, we here propose one alternative way to solve this problem. With samples, we could approximate \(f\) using some integrable estimation \(\hat{f}\), then approximate \(F\) by integral over \(\hat{f}\). UAT guarantees that, for any \(\varepsilon>0\), there exists a neural network \(\psi\) is a good estimation of \(f\), such that \(|\psi-f|\leq\varepsilon/(b-a)\). Then the integral over \(\phi\) satisfies that \(|\int_{a}^{x}\phi(t)dt-F(x)|\leq\varepsilon\), which gives a good estimation of function \(F\).

The problem is how to get the integral of a neural network which will be explained in next section.

## 3 Integral and Algorithms

In this section, our focus is on obtaining the closed-form solutions for integrating neural networks. We first provide explicit integration forms for one layer neural networks with any integrable activation function. Then, we derive a piece-wise structure for the integration of multi-layer neural networks that use ReLU activation function, and propose an forward integral algorithm with a corrector to enhance the accuracy of the integration.

### Integral over one-layer Neural Networks

We start with a trivial case: one layer neural networks \(\psi:\mathbb{R}\rightarrow\mathbb{R}\)

\[\begin{split} y&=W^{(1)}x+b^{(1)},\\ \psi(x)&=W^{(2)}\sigma(y)+b^{(2)}.\end{split} \tag{1}\]

Notice that in this 1-D case, \(W^{(2)}\) is a row vector with length \(n_{1}\) and \(W^{(1)}=(w_{1},\ldots,w_{n_{1}})^{\intercal}\) is a column vector with length \(n_{1}\). Then we state the following Lemma.

**Lemma 3.1**.: _For a one-layer neural network \(\psi\) defined on a closed interval \([a,b]\), the integral of \(\psi\) can be expressed as:_

\[\int_{a}^{x}\psi(t)dt=W^{(2)}z+b^{(2)}(x-a),\]

_where \(z=(z_{1},\ldots,z_{n_{1}})^{\intercal}\in\mathbb{R}^{n_{1}}\) and_

\[z_{i}=\int_{a}^{x}\sigma(w_{i}t)dt.\]
Proof.: By direct calculation, this lemma holds. 

For \(n\)-dimensional case, the result is also similar. We only have to repeat Lemma 3.1 for n time. Let \(\psi\) be an one-layer neural network: \(\mathbb{R}^{n}\rightarrow\mathbb{R}\). Then,

\[\begin{split}&\int_{[a_{1},b_{1}]\times\ldots[a_{n},b_{n}]}\psi(x )dx_{1}\ldots dx_{n}\\ &=\int_{[a_{2},b_{2}]\times\ldots[a_{n},b_{n}]}W^{(2)}[\int_{a_{1 }}^{b_{1}}\sigma(W^{(1)}x+b^{(1)})dx_{1}]+b^{(2)}(b_{1}-a_{1}).\end{split} \tag{2}\]

As long as we have the explicit form of \(\tilde{\sigma}=\int\sigma\), we can write down the integral as

\[\begin{split}&\int_{[a_{1},b_{1}]\times\ldots[a_{n},b_{n}]}\psi(x )dx_{1}\ldots dx_{n}\\ &=\int_{[a_{2},b_{2}]\times\ldots[a_{n},b_{n}]}W^{(2)}\tilde{ \sigma}(\tilde{W}^{(1)}\tilde{x}+\tilde{b}^{(1)})+\tilde{b}^{(2)},\end{split} \tag{3}\]

for some constant term \(\tilde{b}^{(1)}\) relative with \([a_{1},b_{1}]\), \(b^{(1)}\), \(W^{(1)}\) and \(\sigma\). Here \(\tilde{\sigma}\) is the anti-derivative of \(\sigma\), \(\tilde{x}=(x_{2},\ldots,x_{n})\) and \(\tilde{b}^{(2)}=b^{(2)}(b_{1}-a_{1})\). As a result, we can represent the formula with a new activation function \(\tilde{\sigma}\). Thus, repeating Lemma 3.1 leads to the final result for this one-layer case.

[1] focused on the integral of one layer neural network with sigmoid activation function. We demonstrate that our analysis here works for general integrable activation function.

### Numerical Integral over multi-layer Neural Networks

In this section, we will be delving into the integration of a multi-layer ReLU neural network, which is considerably more complex than the one-layer case. In a single-layer neural network, the network can be expressed as a simple weighted sum of several activation functions, making integration relatively straightforward. However, for a multi-layer neural network, the situation is quite different. Even with ReLU activation function, we do not possess explicit knowledge of the areas where the output of a neuron is zero or positive, particularly for neurons in higher layers. As a result, integration of multi-layer neural networks is a challenging task. To address this challenge, we propose a numerical method.

We start with an observation for neural networks. Every ReLU neural network is a piece-wise linear function. Arora et al. [2] also showed that every
piece-wise linear function is also an ReLU neural network. We cite their Theorem as reference.

**Theorem 3.2** ([2]).: _Every \(\mathbb{R}^{n}\rightarrow\mathbb{R}\) ReLU neural network represents a piecewise linear function, and every piece-wise linear function \(\mathbb{R}^{n}\rightarrow\mathbb{R}\) can be represented by a ReLU neural network with at most \(\lceil\log 2(n+1)\rceil+1\) depth_

Then to integral a ReLU neural network is equivalent with to integral a piece-wise linear function without explicitly knowing the break points. In each piece, the neural network can be represented in the form of \(\alpha x+\beta\), where \(\alpha\) and \(\beta\) are coefficients determined by the weights and biases in the network structure. Note that both \(\alpha\) and \(\beta\) could be vector-valued. Knowing these coefficients enables us to compute the integral of the neural network separately in each piece. Thus, we design a forward integral algorithm. We demonstrate that this algorithm works for batch of input and can be accelerated by GPUs.

However, by using a piece-wise integral, it is likely that we cannot obtain a continuous function since there are usually jumps at the breakpoints. Nonetheless, it is worth noting that these jumps are actually constant. To address this issue, we propose a numerical algorithm that corrects these jumps and ensures the continuity of the entire integral.

Figure 1 illustrates the entire neural network integration process, and we will delve into each component in detail in the subsequent sections.

#### 3.2.1 Forward Integral Algorithm

In this section, we will introduce our forward integral algorithm.

As mentioned before, given arbitrary input \(x\), the key to calculate the integral would be the coefficient \(\alpha\) and the constant term \(\beta\). By Theorem 3.2, the output of \(i\)-th layer could be represented as a linear combination of input \(x\) plus a constant term, i.e \(\alpha_{i}x+\beta_{i}\). Note that the output of \(i\)-th layer could be a

Figure 1: Integration Process
vector, which means \(\alpha_{i}\) could be matrix-valued and \(\beta_{i}\) could be vector-valued. Then the output of \((i+1)\)-th layer \(\alpha_{i+1}x+\beta_{i+1}=\sigma(W^{(i+1)}(\alpha_{i}x+\beta_{i})+b^{(i+1)})\). Recall the definition of ReLU, it's equivalent with identity when the input is positive. Otherwise, it always outputs \(0\), having no effect on \(\alpha\) and \(\beta\). Given \(\alpha_{i}\), \(\beta_{i}\), \(W_{(i+1)}\), \(b^{(i+1)}\) and the sign of the term inside \(\sigma\), it's easily to get \(\alpha_{i+1}\) and \(\beta_{i+1}\). Thus, following this procedure for each layer, we could summarize our Forward Integral Algorithm 1. Here \(\mathbb{1}\) denotes indicator function and \(\odot\) denotes Hadamard product.

```
Input:\(x\in\mathbb{R}^{n}\), \(k\)-layer neural network with weights \(W^{(i)}\), bias \(b^{(i)}\) Output: First order coefficient \(\alpha\) and constant term \(\beta\) \(\alpha\gets 0\) \(\beta\gets 0\) \(y\gets x\) for\(i=0\)to\(k\)do if\(i==0\)then \(y\leftarrow\sigma(W^{(i)}y+b^{(i)})\) \(z\leftarrow\mathbb{1}_{\{y>0\}}\) \(\alpha\gets z\odot W^{(0)}\) \(\beta\gets z\odot b^{(0)}\) elseif\(i==k\)then \(\alpha\gets W^{(k)}\times\alpha\) \(\beta\gets W^{(k)}\times\beta+b^{(k)}\) else \(y\leftarrow\sigma(W^{(i)}y+b^{(i)})\) \(z\leftarrow\mathbb{1}_{\{y>0\}}\) \(\alpha\leftarrow(z\odot W^{(i)})\times\alpha\) \(\beta\leftarrow(z\odot W^{(i)})\times\beta+(z\odot b^{(0)})\) endif endfor Return:\(\alpha\) and \(\beta\)
```

**Algorithm 1** Forward Integral Algorithm

We demonstrate here, this algorithm is just a modified version of the forward algorithm to get output by adding Hadamard and indicator operators. As a result, it can be accelerated by GPUs.

As long as we have \(\alpha\) and \(\beta\), we can directly get the integral at the piece where input \(x\) located in, i.e \(\int_{a_{i}}^{x}\psi(x)dx_{i}\) can be represented a Polynomial
term plus a constant term. Since the Polynomial term only depends on \(\alpha\) and \(\beta\), we denote it as Poly[\(\alpha\), \(\beta\)] for convenience.

#### 3.2.2 Numerical Corrector

By Forward Integral Algorithm, we know the integral at the each piece while ignoring a constant term. In this section, we introduce a numerical method to correct the jump at the break point between pieces.

Given an interval, we first select a partition \(\{z_{i}\}_{i=1}^{N}\) of such interval. For each \(z_{i}\). After getting corresponding \(\alpha_{i}\) and \(\beta_{i}\) for each \(z_{i}\) by our forward integral algorithm, we can connect the integral at \(z_{i}\) by adding a constant term \(\text{Poly}[\alpha_{i},\beta i](z_{i})-\text{Poly}[\alpha_{i-1},\beta i-1](z_{i})\). Our corrector algorithm is shown in Algorithm 2.

```
Input: Partition \(\{z_{i}\}_{i=1}^{N}\) on interval \([a,b]\), k-layer neural network \(\psi\) Output:\(\int_{a}^{b}\psi dx_{j}\) \(\alpha_{0}\gets 0\) \(\beta_{0}\gets 0\) \(C_{0}\gets 0\) for\(i=1\)to\(N\)do \(\alpha_{i},\beta_{i}\leftarrow\text{Forward Integral Algorithm}(z_{i})\) \(C_{i}\leftarrow\text{Poly}[\alpha_{i},\beta i](z_{i})-\text{Poly}[\alpha_{i-1},\beta i-1](z_{i})+C_{i-1}\) endfor Return:\(\text{Poly}[\alpha_{N},\beta_{N}](z_{N})-\text{Poly}[\alpha_{0},\beta_{0}](z_{0})+C_{N}\)
```

**Algorithm 2** Numerical Corrector

### Extension to Convolutional Neural Networks and Residual Neural Networks

Convolutional Neural Networks [10] and Residual Neural Networks [6] are more powerful in practice. Our algorithms also work for them, when the activation function is ReLU.

The convolutional layer can be expressed as \(\sigma(\kappa*x+b)\), where \(*\) denotes the convolution operation, \(\kappa\) is the kernel and \(b\) is the bias vector. A convolutional neural network is one consist of convolutional layers and fully connected layers.

A Residual Block is usually consists of convolutional layers, fully-connected layers and a residual connection, i.e \(\Psi(x)=x+\psi(x)\), where \(\psi(x)\) can be expressed as a convolutional neural network. A Residual Neural Network is one
taking structures of Residual Blocks, convolutional layers and fully-connected layers.

We first observe that the convolutional layer and the residual block would not affect the piece-wise structure for ReLU neural networks. As a result, our Forward Integral Algorithm works for both ReLU Convolutional Neural Network and ReLU Residual Neural Network. We demonstrate here the convolutional layer works as same as the fully connected layer in our Forward Integral Algorithm because the convolution operator by kernel \(\kappa\) acts as a weighted linear combination of input. For Residual Neural Networks, we only have to add an all one vector on \(\alpha\) at the end of Residual block when we running forward algorithm.

## 4 Experiments

We constructed extensive experiments to demonstrate that our algorithms work well on numerical integral, compared with traditional numerical algorithms.

We set the domain to be \([0,5]\), \(f(x)=\cos(x)-x^{2}+4-1/(x+1)\). Then the integration \(F(x)=\int_{0}^{x}f(t)dt=\sin(x)-1/3*x^{3}+4x-\log(x+1)\). We trained a 2-layer neural network with width 100 to approximate \(f\) with 51 data points \(\{x_{i},f(x_{i})\}_{i=0}^{50}\), where \(x_{i}\) evenly spaced over \([0,5]\). We set epochs to be 200, learning rate 0.001 with batch size 20, using Stochastic Gradient Descent and mean squared error loss.

For comparisons, we also implemented Euler-Forward method [5], Explicit Runge-Kutta method of order 5(4) (RK45) [4] and Explicit Runge-Kutta method of order 8 (DOP853)[16].

We present our experiment results in Figure 2. In the left figure, we display our neural network's approximation capabilities. The 2-layer neural network exhibits outstanding approximations of the integrand \(f\), and its integral simultaneously provides an excellent approximation of \(F\). In the right figure, we track the approximation error, denoted as \(\|F(x)-\hat{F}(x)\|\), where \(\hat{F}(x)\) represents the estimated integral obtained through numerical algorithms. Notably, the integral computed by the neural network displays the smallest approximation error, highlighting its superior performance in our experiment.
## 5 Discussion

To compute the integral of neural networks, we have developed novel algorithms, including the forward algorithm and a corrector. In conclusion, we outline several promising avenues for future research.

One immediate area of interest is the path integrals of neural networks, offering potential solutions to challenges in fields such as Molecular Dynamics [13] and Quantum Mechanics [7].

Furthermore, Physics-Informed-Neural-Networks (PINNs) [15] have garnered significant attention. While current approaches focus on using derivatives of neural networks to satisfy partial differential equations (PDEs), exploring the integration of neural networks may open up new avenues for addressing these problems.

Neural networks have also demonstrated their efficacy in tackling challenges related to density estimation and Bayesian inference, as demonstrated in prior research [12; 1; 11]. The integrability of neural networks has the potential to significantly enhance these applications, particularly in tasks like expectation and variance estimation.
