HOPE: High-order Polynomial Expansion of Black-box Neural Networks
Despite their remarkable performance, deep neural networks remain mostly
``black boxes'', suggesting inexplicability and hindering their wide
applications in fields requiring making rational decisions. Here we introduce
HOPE (High-order Polynomial Expansion), a method for expanding a network into a
high-order Taylor polynomial on a reference input. Specifically, we derive the
high-order derivative rule for composite functions and extend the rule to
neural networks to obtain their high-order derivatives quickly and accurately.
From these derivatives, we can then derive the Taylor polynomial of the neural
network, which provides an explicit expression of the network's local
interpretations. Numerical analysis confirms the high accuracy, low
computational complexity, and good convergence of the proposed method.
Moreover, we demonstrate HOPE's wide applications built on deep learning,
including function discovery, fast inference, and feature selection. The code
is available at https://github.com/HarryPotterXTX/HOPE.git.
Tingxiong Xiao, Weihang Zhang, Yuxiao Cheng, Jinli Suo
# HOPE: High-order Polynomial Expansion of Black-box Neural Networks

###### Abstract

Despite their remarkable performance, deep neural networks remain mostly "black boxes", suggesting inexpilability and hindering their wide applications in fields requiring making rational decisions. Here we introduce HOPE (High-order Polynomial Expansion), a method for expanding a network into a high-order Taylor polynomial on a reference input. Specifically, we derive the high-order derivative rule for composite functions and extend the rule to neural networks to obtain their high-order derivatives quickly and accurately. From these derivatives, we can then derive the Taylor polynomial of the neural network, which provides an explicit expression of the network's local interpretations. Numerical analysis confirms the high accuracy, low computational complexity, and good convergence of the proposed method. Moreover, we demonstrate HOPE's wide applications built on deep learning, including function discovery, fast inference, and feature selection. The code is available at [https://github.com/HarryPotterXY/HOPE.git](https://github.com/HarryPotterXY/HOPE.git).

 explainable artificial intelligence (XAI), high-order derivative, Taylor expansion, neural network, deep learning.

## 1 Introduction

Deep neural networks have gained widespread adoption due to their ability for universal approximation, as proved by numerous studies [1, 2, 3, 4]. However, deep learning networks are largely considered black boxes that hinder their practical applications. Therefore, understanding the rationale behind predictions is crucial when making relatively logical decisions based on the network output or deciding whether to deploy a new model. This requirement for understanding is particularly important in areas such as clinical decision-making [5][6], drug discovery [7], and physical law identification [8], so people often prioritize models that align with their intuition over accuracy [9, 10, 11]. Therefore, there is a growing need for explainable AI (XAI) approaches to make deep learning more transparent and convincing [12, 13, 14].

XAI can be divided into five categories [7]: feature attribution, instance-based, graph-convolution-based, self-explaining, and uncertainty estimation, among which the feature attribution family calculates the relevance of every input feature for the final prediction and has been the most widely used XAI in recent years. The implementations identifying the feature attribution can be further grouped into the following three categories.

_Gradient-based feature attribution_ approaches measure the impact of a change within an input neighborhood on the change in the output. These methods are mainly inspired by back-propagation [15], and some well-known methods include saliency map [16], Integrated Gradients [17], SmoothGrad [18], local explanation vectors [19], Grad-CAM [20], guided backpropagation [21], LRP [22], and deep Taylor decomposition [23]. The gradient-based methods calculate the first-order derivative of the model as the feature attribution, but omit the higher-order terms, which are important for nonlinear functions.

_Surrogate-model feature attribution_ aims at developing a surrogate explanatory model to the original function, to mirror the computational logic of the original model. The representative surrogate models include LIME [24], DeepLIFT [25], Shapley value [26], LRP [22], SHAP [27], and BETA [28]. However, such surrogate models are mostly linear and suffer from insufficient approximation accuracy.

_Perturbation-based feature attribution_ modifies or removes part of the input to measure the corresponding change in the model output, which reflects the feature importance of the neural network input. Methods like feature masking [29], perturbation analysis [30], response randomization [31], and conditional multivariate models [32] fall into this category. Although intuitive, perturbation-based methods are computationally slow when the number of input features increases [32], and the final result tends to be strongly influenced by the number of perturbed features [33].

In summary, an approach capable of approximating a general deep neural network with high accuracy, low cost and good interpretability is demanded. Here we present HOPE (High-order Polynomial Expansion), an approach to expand a deep neural network into a high-order Taylor polynomial on a reference input. The Taylor expansion is built on calculating the derivatives of the target neural network, which is intrinsically a nonlinear function. We first derive the high-order derivative rule for composite functions and extend this to neural networks to obtain the high-order derivatives quickly and accurately. Our method serves as a gradient-based method and a surrogate model, as it integrates the high-order derivatives and uses nonlinear polynomials to locally approximate the neural network. Our expansion is of high approximation accuracy, and our computational cost is far lower than perturbation-based methods because we can get all derivatives with only one back-propagation.
Actually, in recent years, some researchers have paid attention to the Taylor expansion of neural networks. LRP [22] and deep Taylor decomposition [23] provide a first-order approximation, but they neglect the high-order terms. Morala P et al. [34] explored a mathematical framework for expanding a single-hidden-layer neural network, but it does not apply to multi-layer networks. NN-Poly [35], on the other hand, can infer the Taylor expansion of single-layer neural networks, and obtain the Taylor polynomial of the multi-layer neural network through forward composition, but at extremely high computational complexity. SHoP [36] proposed a Taylor expansion framework to solve high-order partial differential equations, but it only considered fully connected neural networks and adopted a layer-wise approach for the derivation process, rather than a module-wise approach. Differently, our method is similar to back-propagation and computes the derivative of the final output with respect to the intermediate output, propagating from the output layer back to the input layer. Compared to forward composition, our method, HOPE, has significant advantages in terms of accuracy, speed, computational complexity, and memory consumption.

To summarize, the main technical contributions are listed as follows:

* We infer the high-order derivative rule and extend it for calculating the derivatives of deep neural networks with higher accuracy, higher speed, and less memory consumption than the conventional computational-graph-based counterpart.
* We propose a framework to expand a deep neural network into a Taylor series, providing an explicit explanation of the inner workings of the "black-box" model.
* We prove the equivalence between a neural network and its Taylor series, and analyze its convergence condition.
* We demonstrate the wide applications of Taylor expanded deep neural networks, such as function discovery, fast inference, and feature selection.

The paper is structured as follows: In Sec. 2, we present the high-order derivative rule for composite functions. In Sec. 3, we extend this rule to deep neural networks, which allows for efficient and accurate computation of their high-order derivatives. Based on the calculated derivatives, we propose the Taylor expansion framework--HOPE, and analyze its bounds, convergence, and computational complexity in Sec. 4. In Sec. 5, we conducted a series of experiments to show HOPE's significant advantages over the computational-graph-based method, verified the convergence condition, and demonstrated its applications in function discovery, fast inference, and feature selection. Finally, in Sec.6, we analyze the advantages and disadvantages of the proposed approach and discuss some possible future work.

## 2 High-order Derivative Rule for Composite Functions

Considering a composite function

\[\mathbf{y}=h(g(\mathbf{x})), \tag{1}\]

which is constructed by two functions \(\mathbf{z}=g(\mathbf{x})\) and \(\mathbf{y}=h(\mathbf{z})\) to map the input \(\mathbf{x}\in\mathbb{R}^{p}\) to the intermediate state variable \(\mathbf{z}\in\mathbb{R}^{s}\), and then to the final output \(\mathbf{y}\in\mathbb{R}^{o}\) sequentially. Assuming the two functions \(g(\cdot)\) and \(f(\cdot)\) are \(n\)-order differentiable at \(\mathbf{x}0\) and \(\mathbf{z}0=g(\mathbf{x}0)\) respectively, this section will derive the high-order derivative rule for the composite function \(\mathbf{y}=h(g(\mathbf{x}))\) in three systems: Single-Input Single-State Single-Output (SISSSO), Multiple-Input Multiple-State Single-Output (MIMSSO), and Multiple-Input Multiple-State Multiple-Output (MIMSMO).

### _Sissso_

For an SISSSO system, \(\mathbf{x},\mathbf{z},\mathbf{y}\in\mathbb{R}\). From the chain rule, the first three derivatives of \(\mathbf{y}=h(g(\mathbf{x}))\) can be calculated as

\[\left\{\begin{array}{l}\frac{\partial\mathbf{y}}{\partial \mathbf{x}}=\frac{\partial\mathbf{z}}{\partial\mathbf{x}}\frac{\partial \mathbf{y}}{\partial\mathbf{z}},\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}^{2}}=\frac{\partial^{2} \mathbf{y}}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{y}}{\partial \mathbf{z}}+\frac{(\partial\mathbf{z})}{(\partial\mathbf{x})}^{2}\frac{ \partial^{2}\mathbf{y}}{\partial\mathbf{z}^{2}},\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}^{2}}=\frac{\partial^{2} \mathbf{y}}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{y}}{\partial \mathbf{z}}+\frac{3}{\partial\mathbf{x}}\frac{\partial^{2}\mathbf{y}}{ \partial\mathbf{x}^{2}}\frac{\partial^{2}\mathbf{y}}{\partial\mathbf{z}^{2} }+(\frac{\partial\mathbf{z}}{\partial\mathbf{x}})^{3}\frac{\partial^{2} \mathbf{y}}{\partial\mathbf{z}^{3}}.\end{array}\right. \tag{2}\]

For more terms, we can convert \(\frac{\partial^{2}\mathbf{y}}{\partial\mathbf{z}^{2}-1}\frac{\partial \mathbf{x}}{\partial\mathbf{x}}\) to \(\frac{\partial\mathbf{z}}{\partial\mathbf{x}}\frac{\partial^{2}\mathbf{y}}{ \partial\mathbf{z}^{2}}\), and calculate \(\frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}^{2}}\) from \(\{\frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}},k=1,\dots,n\}\) and \(\{\frac{\partial^{2}\mathbf{x}}{\partial\mathbf{x}},k=1,\dots,n\}\). Hence, Eq. (2) turns into following matrix form

\[\left[\begin{array}{l}\frac{\partial\mathbf{y}}{\partial\mathbf{x}}\\ \vdots\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}^{2}}\end{array}\right]=\left[ \begin{array}{l}\frac{\partial\mathbf{z}}{\partial\mathbf{x}}&0&0\\ \frac{\partial^{2}\mathbf{z}}{\partial\mathbf{x}^{2}}&(\frac{\partial \mathbf{z}}{\partial\mathbf{x}})^{2}&0&0\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}^{2}}&3&\frac{\partial^{2} \mathbf{y}}{\partial\mathbf{x}^{2}}&(\frac{\partial\mathbf{z}}{\partial \mathbf{x}})^{3}&0\\ \vdots&\vdots&\vdots&\ddots\end{array}\right]\left[\begin{array}{l}\frac{ \partial\mathbf{y}}{\partial\mathbf{z}}\\ \vdots\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{z}^{2}}\end{array}\right], \tag{3}\]

which can be further abbreviated as

\[\mathbf{v}^{y,x}=\mathbf{M}^{z,x}\mathbf{v}^{y,z}. \tag{4}\]

In this equation \(\mathbf{v}^{y,x},\mathbf{v}^{y,z}\in\mathbb{R}^{n}\) are respectively the vectors composed of derivatives \(\{\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{x}}\}\) and \(\{\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{x}}\}\); \(\mathbf{M}^{z,x}\in\mathbb{R}^{n\times n}\) is the transformation matrix composed of \(\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{x}^{k}}\) and takes a lower triangular form. So far, the calculation of \(h(g(\mathbf{x}))\)'s \(n\)-order derivatives turns into the computation of \(\mathbf{M}^{z,x}\).

From Eq. (3) the \(i\)th and \(i+1\)th terms (\(i<n\)) are respectively

\[\frac{\partial^{i}\mathbf{y}}{\partial\mathbf{x}^{i}}=\sum_{j=1}^{n}\mathbf{M }_{i,j}^{z,x}\frac{\partial^{j}\mathbf{y}}{\partial\mathbf{z}^{j}} \tag{5}\]

and

\[\frac{\partial^{i+1}\mathbf{y}}{\partial\mathbf{x}^{i+1}}=\sum_{j=1}^{n}\mathbf{ M}_{i+1,j}^{z,x}\frac{\partial^{j}\mathbf{y}}{\partial\mathbf{z}^{j}}. \tag{6}\]

Taking derivatives over both sides of Eq. (5) we arrive at

\[\begin{split}\frac{\partial^{i+1}\mathbf{y}}{\partial \mathbf{x}^{i+1}}&=\sum_{j=1}^{n}\frac{\partial}{\partial\mathbf{x}}( \mathbf{M}_{i,j}^{z,x}\frac{\partial^{j}\mathbf{y}}{\partial\mathbf{z}^{j}})\\ &=\sum_{j=1}^{n}\frac{\partial}{\partial\mathbf{x}}\mathbf{M}_{i,j}^{z,x} \frac{\partial^{j}\mathbf{y}}{\partial\mathbf{z}^{j}}+\sum_{j=1}^{n}\frac{ \partial\mathbf{z}}{\partial\mathbf{x}}\mathbf{M}_{i,j}^{z,x}\frac{\partial^{j+1 }\mathbf{y}}{\partial\mathbf{z}^{j+1}}\\ &=\sum_{j=1}^{n}\left(\frac{\partial}{\partial\mathbf{x}}\mathbf{M}_{i,j}^ {z,x}+\frac{\partial\mathbf{z}}{\partial\mathbf{x}}\mathbf{M}_{i,j-1}^{z,x} \right)\frac{\partial^{j}\mathbf{y}}{\partial\mathbf{z}^{j}}\\ &\quad-\frac{\partial\mathbf{z}}{\partial\mathbf{x}}\mathbf{M}_{i,0}^ {z,x}\frac{\partial\mathbf{y}}{\partial\mathbf{z}}+\frac{\partial\mathbf{z}}{ \partial\mathbf{x}}\mathbf{M}_{i,n}^{z,x}\frac{\partial^{n+1}\mathbf{y}}{ \partial\mathbf{z}^{n+1}}.\end{split} \tag{7}\]

Because \(\mathbf{M}_{i,0}^{z,x}=0\) and \(\mathbf{M}_{i,n}^{z,x}=0\)\((i<n)\), Eq. (7) can be simplified into

\[\frac{\partial^{i+1}\mathbf{y}}{\partial\mathbf{x}^{i+1}}=\sum_{j=1}^{n} \left(\frac{\partial}{\partial\mathbf{x}}\mathbf{M}_{i,j}^{z,x}+\frac{ \partial\mathbf{z}}{\partial\mathbf{x}}\mathbf{M}_{i,j-1}^{z,x}\right) \frac{\partial^{j}\mathbf{y}}{\partial\mathbf{z}^{j}}. \tag{8}\]
Compare Eq. (6) and Eq. (8), we can get the recurrence formula of \(\mathbf{M}^{z,x}\) as

\[\left\{\begin{array}{ll}\mathbf{M}_{1,1}^{z,x}=\frac{\partial_{x}}{\partial x}, &i<j\\ \mathbf{M}_{1,j}^{z,x}=0,&i\geq j\\ \mathbf{M}_{i,j}^{z,x}=\frac{\partial\mathbf{M}_{i-1,j}^{z,x}}{\partial x}+ \frac{\partial_{z}}{\partial x}\mathbf{M}_{i-1,j-1}^{z,x},&i\geq j\end{array}\right. \tag{9}\]

which explicitly composes the \(n\)-order transformation matrix \(\mathbf{M}^{z,x}\) in Eq. (4).

### _Mimso_

Unmixed partial derivatives.For a MIMSSO system with \(p\)-dimensional input and \(s\) intermediate states, i.e., \(\mathbf{x}\in\mathbb{R}^{p}\), \(\mathbf{z}\in\mathbb{R}^{s}\), \(\mathbf{y}\in\mathbb{R}\), we can get the first three unmixed partial derivatives of \(\mathbf{y}=h(g(\mathbf{x}))\) as

\[\left\{\begin{array}{ll}\frac{\partial\mathbf{y}}{\partial\mathbf{x}_{i}}= \sum_{j=1}^{s}\left(\frac{\partial\mathbf{z}_{j}}{\partial\mathbf{x}_{i}}, \frac{\partial\mathbf{y}}{\partial\mathbf{z}_{j}}\right),\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}_{i}^{2}}=\sum_{j=1}^{s} \left(\frac{\partial^{2}\mathbf{z}_{j}}{\partial\mathbf{x}_{i}^{2}},\frac{ \partial^{2}\mathbf{y}}{\partial\mathbf{z}_{j}}+(\frac{\partial\mathbf{z}_{j }}{\partial\mathbf{x}_{i}})^{2}\frac{\partial^{2}\mathbf{y}}{\partial\mathbf{z }_{j}^{2}}\right),\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}_{i}^{2}}=\sum_{j=1}^{s} \left(\frac{\partial^{2}\mathbf{z}_{j}}{\partial\mathbf{x}_{i}^{2}},\frac{ \partial\mathbf{y}}{\partial\mathbf{z}_{j}}+3\frac{\partial\mathbf{z}_{j}}{ \partial\mathbf{x}_{i}}\frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}_{i}^{ 2}},\frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}_{j}^{2}}+(\frac{\partial \mathbf{z}_{j}}{\partial\mathbf{x}_{i}})^{3}\frac{\partial^{3}\mathbf{y}}{ \partial\mathbf{z}_{j}^{2}}\right).\end{array}\right. \tag{10}\]

To facilitate derivation, we define an operator \(\beta\) to save the information of their \(k\)-order unmixed partial derivatives

\[\frac{\beta^{k}}{\beta\mathbf{x}^{k}}=\left[\begin{array}{ll}\frac{\partial ^{k}}{\partial\mathbf{x}_{i}^{2}}\\ \vdots\\ \frac{\partial^{k}}{\partial\mathbf{x}_{p}^{k}}\end{array}\right], \tag{11}\]

and the following equations hold

\[\frac{\beta^{k}\mathbf{y}}{\beta\mathbf{x}^{k}}=\left[\begin{array}{ll}\frac {\partial^{k}\mathbf{z}_{j}}{\partial\mathbf{x}_{i}^{2}}&\cdots&\frac{\partial ^{k}\mathbf{z}_{j}}{\partial\mathbf{x}_{i}^{2}}\\ \vdots&\ddots&\vdots\\ \frac{\partial^{k}\mathbf{z}_{j}}{\partial\mathbf{x}_{p}^{k}}&\cdots&\frac{ \partial^{k}\mathbf{z}_{j}}{\partial\mathbf{x}_{p}^{k}}\\ \end{array}\right], \tag{12}\]

Based on the above definitions, Eq. (10) can be rewritten as

\[\left\{\begin{array}{ll}\frac{\partial\mathbf{y}}{\beta\mathbf{x}}=\frac{ \beta\mathbf{x}^{k}}{\beta\mathbf{x}^{k}}\frac{\partial\mathbf{y}}{\beta \mathbf{x}},&\\ \frac{\beta^{2}\mathbf{y}^{2}}{\beta\mathbf{x}^{2}}=\frac{\beta^{2}\mathbf{x}^{ 2}}{\beta\mathbf{x}^{2}}\frac{\partial\mathbf{y}}{\beta\mathbf{x}}+(\frac{ \beta\mathbf{z}^{2}}{\beta\mathbf{x}})^{\circ 2}\frac{\beta^{2}\mathbf{y}}{\beta\mathbf{z}^{2}},&\\ \frac{\beta^{3}\mathbf{y}}{\beta\mathbf{x}^{3}}=\frac{\beta^{3}\mathbf{x}^{ 2}}{\beta\mathbf{x}^{3}}\frac{\partial\mathbf{y}}{\beta\mathbf{z}}+(3\frac{ \beta\mathbf{z}^{2}}{\beta\mathbf{x}})^{\circ 2}\frac{\beta^{2}\mathbf{y}}{\beta\mathbf{x}^{2}} \frac{\partial^{2}\mathbf{y}}{\beta\mathbf{z}^{2}}+(\frac{\beta\mathbf{z}^{2} \mathbf{x}}{\beta\mathbf{x}})^{\circ 3}\frac{\beta^{3}\mathbf{y}}{\beta\mathbf{z}^{3}},&\end{array}\right. \tag{13}\]

where \(\circ k\) is Hadamard power, \((\mathbf{A}^{\circ k})_{i,j}=(\mathbf{A}_{i,j})^{k}\), and \(\odot\) is Hadamard product, \((\mathbf{A}\odot\mathbf{B})_{i,j}=\mathbf{A}_{i,j}\), \(\mathbf{B}_{i,j}\).

Similar to Eq. (3), we turn Eq. (13) into a matrix form

\[\left[\begin{array}{c}\frac{\partial\mathbf{y}}{\beta\mathbf{x}}\\ \vdots\\ \frac{\partial^{2}\mathbf{y}}{\beta\mathbf{x}^{k}}\end{array}\right]=\left[ \begin{array}{ll}\frac{\partial^{2}\mathbf{x}^{2}}{\beta\mathbf{x}^{2}}&0&0& 0\\ \frac{\partial^{2}\mathbf{x}^{2}}{\beta\mathbf{x}^{2}}&(\frac{\partial^{2} \mathbf{x}^{2}}{\beta\mathbf{x}^{2}})^{\circ 2}&0&0\\ \frac{\partial^{2}\mathbf{y}}{\beta\mathbf{x}^{2}}&3\frac{\partial\mathbf{x}^{2}}{ \beta\mathbf{x}^{2}}&(\frac{\partial^{2}\mathbf{x}^{2}}{\beta\mathbf{x}^{2}})^{ \circ 3}&0\\ \vdots&\vdots&\ddots\end{array}\right]\left[\begin{array}{c}\frac{\partial \mathbf{y}}{\beta\mathbf{x}}\\ \vdots\\ \frac{\partial^{3}\mathbf{y}}{\beta\mathbf{z}^{2}}\end{array}\right], \tag{14}\]

which is of a consistent form with Eq. (3), only with scalar elements replaced by matrices, and the power and multiplication turn into Hadamard power \(\circ k\) and Hadamard product \(\odot\) respectively. We further abbreviated the above equation as

\[\mathbf{v}^{y,x}=\mathbf{M}^{z,x}\mathbf{v}^{y,z}. \tag{15}\]

Mixed partial derivatives.The first module of the neural network is mostly a linear layer, such as a fully connected layer or convolutional layer that satisfies

\[\frac{\partial^{k}\mathbf{z}_{j}}{\partial\mathbf{x}_{i_{1}}\ldots\partial \mathbf{x}_{i_{k}}}=0\;(k>1). \tag{16}\]

The first three mixed derivatives of \(\mathbf{y}=h(g(\mathbf{x}))\) are calculated as

\[\left\{\begin{array}{ll}\frac{\partial\mathbf{y}}{\partial\mathbf{x}_{i_{1}} }=\sum_{j=1}^{s}\frac{\partial\mathbf{z}_{j}}{\partial\mathbf{x}_{i_{1}}}\frac{ \partial\mathbf{y}}{\partial\mathbf{z}_{i}},&\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}_{i_{1}}\partial\mathbf{x}_{i_{2}}}= \sum_{j=1}^{s}\frac{\partial\mathbf{z}_{j}}{\partial\mathbf{x}_{i_{1}}}\frac{ \partial\mathbf{z}_{j}}{\partial\mathbf{x}_{i_{2}}}\frac{\partial^{2}\mathbf{y}}{ \partial\mathbf{z}_{j}^{2}},&\\ \frac{\partial^{2}\mathbf{y}}{\partial\mathbf{x}_{i_{1}}\partial\mathbf{x}_{i_{2}}} \partial\mathbf{x}_{i_{3}}=\sum_{j=1}^{s}\frac{\partial\mathbf{z}_{j}}{ \partial\mathbf{x}_{i_{1}}}\frac{\partial\mathbf{z}_{j}}{\partial\mathbf{x}_{i_{2}}} \frac{\partial\mathbf{z}_{j}}{\partial\mathbf{x}_{i_{3}}}\frac{\partial^{3} \mathbf{y}}{\partial\mathbf{z}_{j}^{2}}.&\\ \end{array}\right. \tag{17}\]

Based on the definition in Eqns. (11)(12), the above equations turns into

\[\left\{\begin{array}{ll}\frac{\partial\mathbf{y}}{\beta\mathbf{x}}=\mathbf{Q} _{1}\frac{\partial\mathbf{y}}{\beta\mathbf{x}},&\\ \frac{\partial\mathbf{y}}{\beta\mathbf{x}}\otimes\frac{\partial\mathbf{y}}{\beta \mathbf{x}}=\mathbf{Q}_{2}\frac{\partial^{2}\mathbf{y}}{\beta\mathbf{x}^{2}},&\\ \frac{\partial\mathbf{y}}{\beta\mathbf{x}}\otimes\frac{\partial\mathbf{y}}{\beta \mathbf{x}}\otimes\frac{\partial\mathbf{y}}{\beta\mathbf{x}}=\mathbf{Q}_{3}\frac{ \partial^{2}\mathbf{y}}{\beta\mathbf{x}^{3}},&\\ \mathbf{Q}_{1}=\frac{\beta\mathbf{x}^{2}}{\beta\
## 3 High-order Derivative Rule for Neural Networks

This section mainly introduces the high-order derivative rule of a deep neural network. Since a multiple-output network can be regarded as multiple single-output networks, we consider only the single-output cases. Without loss of generality, we derive the back-propagation of high-order derivatives of the most common modules, with the network structure illustrated in Fig. 1.

Before proceeding with the detailed derivations, we define the following notations. The input is denoted as \(\mathbf{x}\), the output of \(m\)th module as \(\mathbf{y}^{(m)}\), the length of \(\mathbf{y}^{(m)}\) as \(o_{m}\), and the final output as \(\mathbf{y}\). To simplify the expression, we omit the superscripts of \(\mathbf{v}^{\mathbf{y},\mathbf{y}^{(m)}}\) and \(\mathbf{M}^{\mathbf{y}^{(m)},\mathbf{y}^{(m-1)}}\) respectively as \(\mathbf{v}_{m}\) and \(\mathbf{M}_{m}\).

### _Output Unit_

As for the final output \(\mathbf{y}=\mathbf{y}^{(7)}\in\mathbb{R}\), according to the definition in Eq. (12), its \(k\)-order derivatives can be calculated as

(23)

Further from Eq. (15), we can obtain the initial derivative in vector form

(24)

### _Fully Connected Layer_

For a fully connected layer, its input-output relationship is defined as

\[\mathbf{y}^{(m+1)}=\mathbf{W}^{(m+1)}\mathbf{y}^{(m)}+\mathbf{b}^{(m+1)}, \tag{25}\]

where \(\mathbf{W}^{(m+1)}\in\mathbb{R}^{o_{m+1}\times o_{m}}\) is the weight matrix, and \(\mathbf{b}^{(m+1)}\in\mathbb{R}^{o_{m+1}}\) is the bias vector. The \(k\)-order derivative of the \(i\)-th node of \(\mathbf{y}^{(m+1)}\) with respect to the \(j\)-th node of \(\mathbf{y}^{(m)}\) is

\[\frac{\partial^{k}\mathbf{y}_{i}^{(m+1)}}{\partial\mathbf{y}^{(m)}_{j}}^{k}= \left\{\begin{array}{cc}\mathbf{W}_{i,j}^{(m+1)},&k=1\\ 0,&k>1.\end{array}\right. \tag{26}\]

Combining with the definition in Eq. (12), we have

\[\frac{\partial^{k}\mathbf{y}^{(m+1)}}{\partial\mathbf{y}^{(m)}_{i}}^{k}= \left\{\begin{array}{cc}\mathbf{W}^{T},&k=1\\ \mathbf{0},&k>1.\end{array}\right. \tag{27}\]

On the one hand, we can get all the unmixed partial derivatives \(\{\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{y}_{i}^{(m)}}^{k},i=1,\dots, o_{m}\}\) by calculating \(\mathbf{v}_{m}\)

\[\mathbf{v}_{m}=\mathbf{M}_{m+1}\mathbf{v}_{m+1}, \tag{28}\]

with the transformation matrix \(\mathbf{M}_{m+1}\) in Eq. (14) rewritten as a block diagonal matrix

\[\mathbf{M}_{m+1}=diag(\mathbf{W}^{T},\mathbf{W}^{T^{\mathcal{O}2}},\dots, \mathbf{W}^{T^{\mathcal{O}n}}). \tag{29}\]

On the other hand, we can obtain all the mixed partial derivatives \(\{\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{y}_{i_{1}}^{(m)},\dots, \partial\mathbf{y}_{i_{k}}^{(m)}},i_{1},\dots,i_{k}=1,\dots,o_{m}\}\) by calculating \(\mathbf{v}_{m}^{\star}\) from

\[\mathbf{v}_{m}^{\star}=\mathbf{M}^{\star}_{m+1}\mathbf{v}_{m+1}. \tag{30}\]

In this equation, the transformation matrix \(\mathbf{M}^{\star}_{m+1}\) defined in Eq. (19) can be also rewritten as a block diagonal matrix

\[\mathbf{M}^{\star}_{m+1}=diag(\mathbf{Q}_{1},\mathbf{Q}_{2},\dots,\mathbf{Q} _{n}), \tag{31}\]

with \(\mathbf{Q}_{1}=\mathbf{W}^{T}\), \(\mathbf{Q}_{k}=\left(\mathbf{W}^{T}\otimes\mathbf{1}^{(o_{m}{}^{k-1})}\right) \odot(\mathbf{1}^{o_{m}}\otimes\mathbf{Q}_{k-1})\).

### _Convolutional Layer_

The input-output relationship of a convolutional layer can be described as

\[\mathbf{y}^{(m+1)}=\mathbf{y}^{(m)}\ast\mathbf{F}^{(m+1)}, \tag{32}\]

where \(\mathbf{F}^{(m+1)}\) is a convolutional kernel, and \(\ast\) denotes convolutional operation. Although a convolutional layer can be regarded as a fully connected layer with a sparse weight matrix and zero bias, taking derivative is both time consuming and memory demanding if transformed into a fully connected layer. Therefore, we derive the high-order derivative rule on the convolution representation.

The \(u\)-th output sums over the product of some elements in \(\mathbf{y}^{(m)}\) and all elements in \(\mathbf{F}^{(m+1)}\),

\[\mathbf{y}_{u}^{(m+1)}=\sum_{v}\mathbf{y}_{v}^{(m)}\mathbf{F}_{u,v}^{(m+1)}, \tag{33}\]

Fig. 1: **The framework for calculating the high-order derivatives of HOPE.** The network has seven modules: fully connected layer (fc), activation function (act), unflatten module (unflatten), convolutional layer (conv), pooling layer (pooling), flatten module (flatten), and single-output fully connected layer. The output of \(m\)th module is \(\mathbf{y}^{(m)}\) and the final output \(\mathbf{y}=\mathbf{y}^{(7)}\). HOPE is similar to back-propagation in that we calculate the derivatives of \(\mathbf{y}\) with respect to the intermediate output \(\mathbf{y}^{(m)}\), starting from the output layer and moving backward towards the input layer. In the intermediate layer, we only need to calculate \(\mathbf{v}_{m}\), which contains all the unmixed partial derivatives. In the input layer, we can calculate \(\mathbf{v}_{m}^{\star}\), which contains all the mixed partial derivatives.
where \(\mathbf{F}_{u,v}^{(m+1)}\) is the weight between \(\mathbf{y}_{u}^{(m+1)}\) and \(\mathbf{y}_{v}^{(m)}\). We can calculate the derivatives as

\[\frac{\partial\mathbf{y}}{\partial\mathbf{y}_{v}^{(m)}}=\sum_{u}\frac{\partial \mathbf{y}}{\partial\mathbf{y}_{u}^{(m+1)}}\frac{\partial\mathbf{y}_{u}^{(m+1) }}{\partial\mathbf{y}_{v}^{(m)}}=\sum_{u}\frac{\partial\mathbf{y}}{\partial \mathbf{y}_{u}^{(m+1)}}\mathbf{F}_{u,v}^{(m+1)}. \tag{34}\]

Because \(\frac{\partial\mathbf{y}_{v}^{(m+1)^{k}}}{\partial\mathbf{y}_{v}^{(m)}}=0\) (\(k>1\)), the high-order derivatives are

\[\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{y}_{v}^{(m)}}=\sum_{u}\frac{ \partial^{k}\mathbf{y}}{\partial\mathbf{y}_{u}^{(m+1)^{k}}}\left(\frac{ \partial\mathbf{y}_{u}^{(m+1)}}{\partial\mathbf{y}_{v}^{(m)}}\right)^{k}= \sum_{u}\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{y}_{u}^{(m+1)^{k}}} \mathbf{F}_{u,v}^{(m+1)^{k}}. \tag{35}\]

Taking the matrix form, Eq. (34) turn into

\[\frac{\partial\mathbf{y}}{\partial\mathbf{y}^{(m)}}=\frac{\beta\mathbf{y}}{ \partial\mathbf{y}^{(m+1)}}*rot180(\mathbf{F}). \tag{36}\]

Comparing Eq. (34) and Eq. (35), we can also get the matrix form of the high-order derivatives as

\[\frac{\partial^{k}\mathbf{y}}{\partial\mathbf{y}^{(m)}}^{k}=\frac{\beta^{k} \mathbf{y}}{\beta\mathbf{y}^{(m+1)^{k}}}*rot180(\mathbf{F}^{\circ k}). \tag{37}\]

Since the input of a convolution layer is usually image-like data, we only calculate unmized partial derivatives as

\[\mathbf{v}_{m}=\left[\begin{array}{cc}\left(\frac{\partial\mathbf{y}}{ \partial\mathbf{y}^{(m)}}\right)^{T}&\dots&\left(\frac{\beta^{n}\mathbf{y}}{ \partial\mathbf{y}^{(m)}}\right)^{T}\end{array}\right]^{T} \tag{38}\]

One can also calculate the mixed partial derivatives by converting the convolutional layer into a fully connected counterpart and then applying Eq. (30) directly.

### _Nonlinear Activation Layer_

Consider a nonlinear activation layer \(\mathbf{y}^{(m+1)}=\sigma(\mathbf{y}^{(m)})\), where \(\mathbf{y}_{i}^{(m+1)}=\sigma(\mathbf{y}_{i}^{(m)})\), we have

\[\frac{\partial^{k}\mathbf{y}_{i}^{(m+1)}}{\partial\mathbf{y}_{j}^{(m)}}^{k}= \left\{\begin{array}{cc}\sigma^{(k)}(\mathbf{y}_{j}^{(m)}),&i=j\\ 0,&i\neq j\end{array}\right. \tag{39}\]

where \(\sigma^{(k)}(\cdot)\) is the \(k\)-order derivative of this activation function. According to the definition in Eq. (12), we can get

\[\frac{\beta^{k}\mathbf{y}^{(m+1)^{T}}}{\beta\mathbf{y}^{(m)}}^{k}= \left[\begin{array}{ccc}\frac{\partial^{k}\mathbf{y}_{i}^{(m+1)}}{\partial \mathbf{y}_{i}^{(m)k}}&\dots&\frac{\partial^{k}\mathbf{y}_{i}^{(m+1)}}{ \partial\mathbf{y}_{i}^{(m)k}}\\ \vdots&\ddots&\vdots\\ \frac{\partial^{k}\mathbf{y}_{i}^{(m+1)}}{\partial\mathbf{y}_{i}^{(m)k}}& \dots&\frac{\partial^{k}\mathbf{y}_{i}^{(m+1)}}{\partial\mathbf{y}_{i}^{(m)k} }\end{array}\right] \tag{40}\] \[=diag\left(\sigma^{(k)}(\mathbf{y}_{1}^{(m)}),\dots,\sigma^{(k)}( \mathbf{y}_{o}^{(m)})\right).\]

After calculating \(\frac{\beta^{k}\mathbf{y}^{(m+1)^{T}}}{\beta\mathbf{y}^{(m)k}}\), we can further obtain the transformation matrix \(\mathbf{M}_{m+1}\) from Eq. (14) and the unmixed partial derivative vector \(\mathbf{v}_{m}\) from Eq.(28). The only left question is how to calculate the value of \(\sigma^{(k)}(x)\) with \(x\in\mathbb{R}\). Here we give the calculation of several widely used activation functions.

\(\bullet\)**Sine \(\boldsymbol{\sigma(x)=sin(x)}\).** The derivatives are

\[\sigma^{(k)}(x)=\left\{\begin{array}{ll}cos(x),&k\ mod\ 4=1\\ -sin(x),&k\ mod\ 4=2\\ -cos(x),&k\ mod\ 4=3\\ sin(x).&k\ mod\ 4=0\end{array}\right. \tag{41}\]

\(\bullet\)**ReLU \(\boldsymbol{\sigma(x)=max(0,x)}\).** The derivatives are

\[\frac{\partial^{k}\sigma(x)}{\partial x^{k}}=\left\{\begin{array}{ll}1,&if \ k=1\ and\ x>0\\ 0.&else\end{array}\right. \tag{42}\]

\(\bullet\)**Sigmoid \(\boldsymbol{\sigma(x)=\frac{e^{x}}{1+e^{x}}}\).** The first derivative is

\[\sigma^{(1)}(x)=\frac{e^{x}}{(1+e^{x})^{2}}=\sigma(x)-\sigma(x)^{2}. \tag{43}\]

Further, we can express \(\sigma^{(k)}(x)\) as the form containing only \(\sigma(x)\), i.e.,

\[\sigma^{(2)}(x) =\sigma^{(1)}(x)-2\sigma(x)\sigma^{(1)}(x)\] \[=[\sigma(x)-\sigma(x)^{2}]-2\sigma(x)[\sigma(x)-\sigma(x)^{2}] \tag{44}\] \[=\sigma(x)-3\sigma(x)^{2}+2\sigma(x)^{3}.\]

After calculating the other derivatives and organizing them into matrix form, we have

\[\left[\begin{array}{c}\sigma(x)\\ \sigma^{(1)}(x)\\ \sigma^{(2)}(x)\\ \vdots\\ \sigma^{(n)}(x)\end{array}\right]=\left[\begin{array}{cccc}1&0&0&0\\ 1&-1&0&0\\ 1&-3&2&0\\ \vdots&\vdots&\vdots&\ddots\end{array}\right]\left[\begin{array}{c}\sigma(x) \\ \sigma(x)^{2}\\ \vdots\\ \sigma(x)^{n+1}\end{array}\right]. \tag{45}\]

This square matrix takes a lower triangular form, and we abbreviate it as \(B\in\mathbb{R}^{n+1\times n+1}\). Similar to the derivation of the transformation matrix in Eq. (3), from above equation, the \(k\)-th and \((k+1)\)-th derivatives are

\[\sigma^{(k)}(x)=\sum_{i=1}^{k+1}B_{k+1,i}\sigma(x)^{i}, \tag{46}\]

\[\sigma^{(k+1)}(x)=\sum_{i=1}^{k+2}B_{k+2,i}\sigma(x)^{i}. \tag{47}\]

Taking derivatives over both sides of Eq. (46) yields

\[\sigma^{(k+1)}(x) =\sum_{i=1}^{k+1}iB_{k+1,i}\sigma(x)^{i-1}\sigma^{(1)}(x)\] \[=\sum_{i=1}^{k+1}iB_{k+1,i}\sigma(x)^{i-1}[\sigma(x)-\sigma(x)^{2}]\] \[=\sum_{i=1}^{k+1}iB_{k+1,i}\sigma(x)^{i}-\sum_{i=1}^{k+1}iB_{k+1,i} \sigma(x)^{i+1}\] \[=\sum_{i=1}^{k+1}iB_{k+1,i}\sigma(x)^{i}-\sum_{i=2}^{k+2}(i-1)B_{k+1,i-1}\sigma(x)^{i}\] \[=\sum_{i=1}^{k+2}[iB_{k+1,i}-(i-1)B_{k+1,i-1}]\sigma(x)^{i}\]

Comparing Eq. (47) and Eq. (48), we arrive at the following relationship

\[\left\{\begin{array}{ll}B_{1,1}=1,&i<j\\ B_{i,j}=jB_{i-1,j}-(j-1)B_{i-1,j-1}.&i\geq j\end{array}\right. \tag{49}\]
* **Tanh \(\mathbf{\sigma(x)}=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}\).** The first and second derivatives are \[\begin{split}\sigma^{(1)}(x)&=1-(\frac{e^{x}-e^{-x} }{e^{x}+e^{-x}})^{2}=1-\sigma(x)^{2},\\ \sigma^{(2)}(x)&=-2\sigma(x)\sigma^{(1)}(x)=-2\sigma( x)+2\sigma(x)^{3}.\end{split}\] (50)

Organize it into matrix form:

\[\left[\begin{array}{c}1\\ \sigma(x)\\ \sigma^{(1)}(x)\\ \sigma^{(2)}(x)\\ \vdots\\ \sigma^{(n)}(x)\end{array}\right]=\left[\begin{array}{ccccc}1&0&0&0&0\\ 0&1&0&0&0\\ 1&0&-1&0&0\\ 0&-2&0&2&0\\ \vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right]\left[\begin{array}{c}1 \\ \sigma(x)\\ \sigma(x)^{2}\\ \sigma(x)^{3}\\ \vdots\\ \sigma(x)^{n+1}\end{array}\right]. \tag{51}\]

This square matrix takes a lower triangular form, and we abbreviate it as \(\mathbf{C}\in\mathbb{R}^{n+2\times n+2}\). From above equation, the \(k\)-th and \((k+1)\)-th derivatives are

\[\sigma^{(k)}(x)=\sum_{i=1}^{k+2}C_{k+2,i}\sigma(x)^{i-1}, \tag{52}\]

\[\sigma^{(k+1)}(x)=\sum_{i=1}^{k+3}C_{k+3,i}\sigma(x)^{i-1}. \tag{53}\]

Taking derivatives over both sides of Eq. (52)

\[\begin{split}&\sigma^{(k+1)}(x)=\sum_{i=1}^{k+2}(i-1)C_{k+2,i} \sigma(x)^{i-2}\sigma^{(1)}(x)\\ &=\sum_{i=1}^{k+2}(i-1)C_{k+2,i}\sigma(x)^{i-2}[1-\sigma(x)^{2}] \\ &=\sum_{i=1}^{k+2}(i-1)C_{k+2,i}\sigma(x)^{i-2}-\sum_{i=1}^{k+2}(i- 1)C_{k+2,i}\sigma(x)^{i}\\ &=\sum_{i=0}^{k+1}iC_{k+2,i+1}\sigma(x)^{i-1}-\sum_{i=2}^{k+3}(i-2)C_{ k+2,i-1}\sigma(x)^{i-1}\\ &=\sum_{i=1}^{k+3}[iC_{k+2,i+1}-(i-2)C_{k+2,i-1}]\sigma(x)^{i-1} \end{split} \tag{54}\]

Comparing Eq. (53) and Eq. (54), we have

\[\left\{\begin{array}{ll}C_{1,1}=1,C_{2,1}=0,C_{2,2}=1,\\ C_{i,j}=0,&i<j\\ C_{i,j}=jC_{i-1,j+1}-(j-2)C_{i-1,j-1}.&i\geq j,i\geq 3\end{array}\right. \tag{55}\]

### _Pooling Layer_

The pooling layer divides the input into blocks and takes either the maximal or average value of each block as the output. To demonstrate the expansion of the pooling layer, we take the simple 2-D pooling layer shown in Fig. 2 as an example.

* **Max pooling** layer outputs the maximum entry of each block, i.e., \[\begin{split}\mathbf{y}_{1}^{(m+1)}&=max(\mathbf{y}_{1}^{(m )},\mathbf{y}_{2}^{(m)},\mathbf{y}_{5}^{(m)},\mathbf{y}_{6}^{(m)})=\mathbf{y} _{idx}^{(m)},\\ \mathbf{y}_{2}^{(m+1)}&=max(\mathbf{y}_{3}^{(m)},\mathbf{y}_{4}^{(m )},\mathbf{y}_{7}^{(m)},\mathbf{y}_{8}^{(m)})=\mathbf{y}_{idx_{2}}^{(m)},\\ \mathbf{y}_{3}^{(m+1)}&=max(\mathbf{y}_{9}^{(m)},\mathbf{y}_{10}^{(m )},\mathbf{y}_{13}^{(m)},\mathbf{y}_{14}^{(m)})=\mathbf{y}_{idx_{3}}^{(m)}, \\ \mathbf{y}_{4}^{(m+1)}&=max(\mathbf{y}_{11}^{(m)},\mathbf{y}_{12}^{(m )},\mathbf{y}_{15}^{(m)},\mathbf{y}_{16}^{(m)})=\mathbf{y}_{idx_{4}}^{(m)}, \end{split}\] (56) where \(idx_{1}\in\{1,2,5,6\},idx_{2}\in\{3,4,7,8\},idx_{3}\in\{9,10,13,14\},idx_{4} \in\{11,12,15,16\}\) are the subscripts of the maximum inputs in four blocks, respectively. Given \(\frac{\partial\mathbf{y}^{k}}{\partial\mathbf{y}_{i}^{(m+1)^{k}}}\), we can calculate the derivatives with respect to \(\mathbf{y}^{(m)}\) as \[\frac{\partial\mathbf{y}^{k}}{\partial\mathbf{y}_{j}^{(m)}}=\left\{\begin{array}[] {ll}\frac{\partial\mathbf{y}^{k}}{\partial\mathbf{y}_{i}^{(m+1)^{k}}},&if \ j=idx_{i}\\ 0.&else\end{array}\right.\] (57) Therefore, we only need to record the subscripts of the corresponding maximum inputs and assign \(\frac{\partial\mathbf{y}^{k}}{\partial\mathbf{y}_{i}^{(m+1)^{k}}}\) to \(\frac{\partial\mathbf{y}^{k}}{\partial\mathbf{y}_{i}^{(m)}}\). When the stride is less than the kernel size, one input may be related to multiple outputs and the derivative can be written as the following general formula \[\frac{\partial\mathbf{y}^{k}}{\partial\mathbf{y}_{j}^{(m)}}=\sum_{\begin{subarray} {c}i=1\\ j=idx_{i}\end{subarray}}^{o_{m+1}}\frac{\partial\mathbf{y}^{k}}{\partial\mathbf{y }_{i}^{(m+1)^{k}}}.\] (58)
* **Average pooling** layer takes the average over each block as the output, i.e., \[\begin{split}\mathbf{y}_{1}^{(m+1)}&=\frac{1}{4}(\mathbf{y }_{1}^{(m)}+\mathbf{y}_{2}^{(m)}+\mathbf{y}_{5}^{(m)}+\mathbf{y}_{6}^{(m)}), \\ \mathbf{y}_{2}^{(m+1)}&=\frac{1}{4}(\mathbf{y}_{3}^{(m)}+\mathbf{y }_{4}^{(m)}+\mathbf{y}_{7}^{(m)}+\mathbf{y}_{8}^{(m)}),\\ \mathbf{y}_{3}^{(m+1)}&=\frac{1}{4}(\mathbf{y}_{9}^{(m)}+\mathbf{y }_{10}^{(m)}+\mathbf{y}_{13}^{(m)}+\mathbf{y}_{14}^{(m)}),\\ \mathbf{y}_{4}^{(m+1)}&=\frac{1}{4}(\mathbf{y}_{11}^{(m)}+\mathbf{y }_{12}^{(m)}+\mathbf{y}_{15}^{(m)}+\mathbf{y}_{16}^{(m)}).\end{split}\] (59) As shown in Fig. (3), we can regard a average pooling layer as a special convolutional layer, with the input channel being the same as the output channel. After getting the equivalent convolutional layer of the average pooling layer, we can obtain the derivatives easily from the high-order derivative rule of the convolutional layer.

Fig. 2: **The illustration of a simple 2-D pooling layer downsampling \(\mathbf{y}^{(m+1)}\) to \(\mathbf{y}^{(m)}\).** Here the kernel size is \(2\times 2\), and the stride is \(2\times 2\).
## 4 Taylor Expansion of Neural Networks

As shown in Fig. 1, applying a forward propagation and back-propagation, we can get a mixed partial derivative vector \(\mathbf{v}_{0}^{*}\) on a reference input \(\mathbf{x}0\), which contains all the \(n\)-order derivatives \(\{\frac{\partial^{2}\mathbf{y}^{*}}{\partial\mathbf{x}_{i_{1}}\ldots\partial \mathbf{x}_{i_{k}}},i_{1},\ldots,i_{k}=1,\ldots,p\}\).

Denote the neural network as \(\mathbf{y}=f(\mathbf{x};\theta)\), with \(\theta\) being the network parameters, the \(n\)-order Taylor polynomial of \(f(\mathbf{x};\theta)\) can be expressed as

\[\begin{split} f(\mathbf{x})&=f(\mathbf{x}0;\theta )+\sum_{i=1}^{p}\frac{\frac{\partial f(\mathbf{x}\theta)}{\partial\mathbf{x}_ {i_{1}}\ldots\partial\mathbf{x}_{i_{k}}}|\mathbf{x}0}{1!}\Delta\mathbf{x}_{i} +\ldots\\ &+\sum_{i_{1},\ldots,i_{n}=1}^{p}\frac{\frac{\partial^{2}\mathbf{y }^{*}}{\partial\mathbf{x}_{i_{1}}\ldots\partial\mathbf{x}_{i_{n}}}|\mathbf{x }0}{n!}\Delta\mathbf{x}_{i_{1}}\ldots\Delta\mathbf{x}_{i_{n}},\end{split} \tag{60}\]

where \(\Delta\mathbf{x}=\mathbf{x}-\mathbf{x}0\), \(\frac{\partial^{k}f(\mathbf{x}\theta)}{\partial\mathbf{x}_{i_{1}}\ldots\partial \mathbf{x}_{i_{n}}}|\mathbf{x}0}|\mathbf{x}0\) is a \(k\)-order partial derivative on the reference input \(\mathbf{x}0\). When all the modules are infinitely differentiable, the network is equivalent to an infinite Taylor polynomial. If the high-order derivatives are much smaller than the low-order derivatives, we can ignore the high-order terms and Eq. (60) can provide an accurate and explicit explanation for the prediction.

### _Upper and Lower Bounds of Network and Taylor Polynomial_

For simplicity, we take the 1-D neural network as an example (i.e., \(\mathbf{x},\mathbf{x}0\in\mathbb{R}\)). Suppose \(f(\mathbf{x};\theta)\) has \(n\)-order continuous derivatives in the interval \([a,b]\) and \(\mathbf{x}0\in(a,b)\), for \(\forall\mathbf{x}\in[a,b]\), \(\exists\xi\in[\min(\mathbf{x},\mathbf{x}0),\max(\mathbf{x},\mathbf{x}0)]\), \(s.t.\)

\[f(\mathbf{x};\theta)=\sum_{k=0}^{n-1}\frac{f(\mathbf{x}0;\theta)^{(k)}}{k!} \Delta\mathbf{x}^{k}+\frac{f(\xi;\theta)^{(n)}}{n!}\Delta\mathbf{x}^{n}, \tag{61}\]

where \(f(\mathbf{x};\theta)^{(k)}\) is the \(k\)-order derivative with respect to the input \(\mathbf{x}\), and \(\frac{f(\xi;\theta)^{(n)}}{n!}\Delta\mathbf{x}^{n}\) is a \(n\)-order Lagrange remainder. After applying the proposed high-order Taylor expansion at \(\mathbf{x}0\), the \(f(\mathbf{x};\theta)\)'s \(n\)-order Taylor polynomial is derived as

\[f(\mathbf{x})=\sum_{k=0}^{n-1}\frac{f(\mathbf{x}0;\theta)^{(k)}}{k!}\Delta \mathbf{x}^{k}+\frac{f(\mathbf{x}0;\theta)^{(n)}}{n!}\Delta\mathbf{x}^{n}. \tag{62}\]

Setting

\[\begin{split} f_{1}(\mathbf{x})&=\sum_{k=0}^{n-1} \frac{f(\mathbf{x}0;\theta)^{(k)}}{k!}\Delta\mathbf{x}^{k}+\frac{\max_{ \mathbf{x}\in[a,b]}f(\mathbf{x};\theta)^{(n)}}{n!}\Delta\mathbf{x}^{n},\\ f_{2}(\mathbf{x})&=\sum_{k=0}^{n-1}\frac{f(\mathbf{ x}0;\theta)^{(k)}}{k!}\Delta\mathbf{x}^{k}+\frac{\min_{\mathbf{x}\in[a,b]}f( \mathbf{x};\theta)^{(n)}}{n!}\Delta\mathbf{x}^{n},\end{split} \tag{63}\]

we can provide the upper and lower boundaries of the network and the Taylor polynomial as

\[\begin{split} f_{d}(\mathbf{x})&\leq f(\mathbf{x}; \theta)\leq f_{u}(\mathbf{x}),\\ f_{d}(\mathbf{x})&\leq f(\mathbf{x})\leq f_{u}( \mathbf{x}),\end{split} \tag{64}\]

where

\[\begin{split} f_{u}(\mathbf{x})&=\max\left(f_{1}( \mathbf{x}),f_{2}(\mathbf{x})\right),\\ f_{d}(\mathbf{x})&=\min\left(f_{1}(\mathbf{x}),f_{2}( \mathbf{x})\right).\end{split} \tag{65}\]

Further, we can provide an upper boundary of the approximation error as

\[\begin{split}|f(\mathbf{x})&-f(\mathbf{x};\theta)| \leq|f_{u}(\mathbf{x})-f_{d}(\mathbf{x})|=|f_{1}(\mathbf{x})-f_{2}(\mathbf{x} )|\\ &=\frac{\max_{\mathbf{x}\in[a,b]}f(\mathbf{x};\theta)^{(n)}-\min_{ \mathbf{x}\in[a,b]}f(\mathbf{x};\theta)^{(n)}}{n!}|\Delta\mathbf{x}|^{n}.\end{split} \tag{66}\]

From the above equation, we can see that the approximation performance is closely related to three factors.

1. The range of \(f(\mathbf{x};\theta)^{(n)}\) in the interval \([a,b]\): \(r=\max_{\mathbf{x}\in[a,b]}f(\mathbf{x};\theta)^{(n)}-\min_{\mathbf{x}\in[a,b]}f (\mathbf{x};\theta)^{(n)}\). Specifically, \(r\) decreases with \(f(\mathbf{x};\theta)^{(n+1)}\), and a smaller \(r\) leads to a smaller approximation error.
2. The order of Taylor polynomial: \(n\). When \(n\) is large, the growth rate of \(n!\) far exceeds the growth speed of \(|\Delta\mathbf{x}|^{n}\), resulting in a decrease in \(|f(\mathbf{x})-f(\mathbf{x};\theta)|\).
3. The distance from \(\mathbf{x}\) to the reference point \(\mathbf{x}0\): \(|\Delta\mathbf{x}|\). The closer \(\mathbf{x}\) is to \(\mathbf{x}0\), the smaller the approximation error tends to be.

It should be noted that Eqns. (64)(66) are theoretical bounds for the error between neural networks and Taylor polynomial, and they hold only when the first \(n\) derivatives are accurate enough. However, in practice, due to the precision limitations of computer storage and computation, it cannot be guaranteed that the upper bounds for the approximation error can always be estimated in all cases.

### _Convergence Analysis of HOPE_

If \(\exists n\), \(s.t.\)\(\forall\mathbf{x}\in[a,b]\), \(|f(\mathbf{x};\theta)^{(n+1)}|\to 0\), we have

\[\max_{\mathbf{x}\in[a,b]}f(\mathbf{x};\theta)^{(n)}-\min_{\mathbf{x}\in[a,b]}f (\mathbf{x};\theta)^{(n)}\to 0.\]

Further from Eq. (66), the approximation error \(|f(\mathbf{x})-f(\mathbf{x};\theta)|\to 0\) and then the Taylor polynomial will converge to the target neural network, i.e., \(f(\mathbf{x})\to f(\mathbf{x};\theta)\). In this section, we will analyze the condition of \(|f(\mathbf{x};\theta)^{(n+1)}|\to 0\)

Fig. 3: **The illustration of a convolutional layer equivalent to the average pooling layer.** Here the number of input and output channels are both 3, and 3 kernels are involved. Convolution the input with kernel \(\#1\) acts as an average pooling applied specifically to the first input channel, and similar equivalence hold for the other two kernels.
From Eqns. (29)(31)(37), we can see that the \(k\)-order derivatives are related with \(\mathbf{W}_{i,j}^{k}\).

\[|f(\mathbf{x};\theta)^{(k)}|\propto|\mathbf{W}_{i,j}|^{k}. \tag{67}\]

When the elements in \(\mathbf{W}\) are concentrated near 0, high-order derivatives are more likely to approach 0. When the parameters are located far from 0, high-order derivatives may become increasingly larger, and thus the Taylor polynomial diverges, i.e., \(f(\mathbf{x})\nrightarrow f(\mathbf{x};\theta)\).

\[\begin{split}\lim_{\begin{subarray}{c}|\mathbf{W}_{i,j}| \to 0\\ k\rightarrow\infty\end{subarray}}|f(\mathbf{x};\theta)^{(k)}|=0,\\ \lim_{\begin{subarray}{c}|\mathbf{W}_{i,j}|>1\\ k\rightarrow\infty\end{subarray}}|f(\mathbf{x};\theta)^{(k)}|=+\infty.\end{split} \tag{68}\]

The above analysis tells that the parameter distribution of each layer has a great influence on the convergence of Taylor expansion. Therefore, we can refer to the above rules to design network structures or impose constraints on the network parameters during network training to achieve deep neural networks with high-order Taylor approximation, which facilitates leveraging the advantages of such explicit expansion. We will verify this conclusion in the Experiment section.

### _Time Complexity Analysis of HOPE_

As a fundamental building block and one of the most time-consuming operations in deep learning [15], back-propagation is implemented via automatic differentiation on the computational graph of the neural network in most deep learning frameworks, such as Autograd [37]. Here, we analyze and compare the time complexity of the computational-graph-based method and HOPE for a \(p\)-D neural network.

When calculating high-order derivatives, the length of the computational graph increases exponentially with the order of the derivative at base 2, because for each node one needs to accumulate the local derivatives along all the paths from the node to the input, resulting in an exponential increase in the number of nodes in the computational graph. For computational-graph-based method, mathematically, there are \(p^{k}\)\(k\)-order derivatives, and the length of their computational graphs is \(2^{k-1}\), so the time complexity of the computational graph is

\[T(n)=\sum_{k=1}^{n}p^{k}2^{k-1}\sim\mathcal{O}((2p)^{n}). \tag{69}\]

Differently, HOPE obtains all the derivatives at one time, with the main calculations lying in calculating the transformation matrix \(\mathbf{M}\) and conducting back-propagation. Since \(\mathbf{M}\) is a lower triangular matrix and the block matrices in \(k\)-th row need \(k\) operations, the complexity of calculating \(\mathbf{M}\) is \(T(n)=\sum_{k=1}^{n}k^{2}=\frac{n(n+1)(2n+1)}{6}\sim\mathcal{O}(n^{3})\). For linear layers, \(\mathbf{M}\) turns into a diagonal matrix and the complexity reduces to \(T(n)=\sum_{k=1}^{n}k=\frac{n(n+1)}{2}\sim\mathcal{O}(n^{2})\). For mixed partial derivatives in Eq. (21), \(\mathbf{M}^{*}\) is a diagonal matrix and the size of \(\mathbf{Q}_{k}\) is \(p^{k-1}\) times larger than \(\mathbf{W}\), and the complexity is about \(T(n)=\sum_{k=1}^{n}p^{k-1}=\frac{1-p^{2}}{1-p}\sim\mathcal{O}(p^{n})\). Therefore, the complexity of HOPE is

\[\mathcal{O}(n^{2})<T(n)<\mathcal{O}(p^{n}). \tag{70}\]

## 5 Experiments

In this section, we quantitatively demonstrate HOPE's significant advantages in terms of accuracy, speed, and memory consumption. We further explored the influence of the target network's parameter distribution on the convergence of its Taylor series, which verified the conclusion in Section 4.2. Besides, we also visualized the upper and lower bounds of a neural network and its Taylor polynomial at increasing order of polynomial terms. Finally, we conducted three experiments to show HOPE's applications, including function discovery, low-cost inference, and feature selection.

### _Approximation Accuracy_

Among all the computational-graph-based methods, Autograd [37] stands out as the most widely used and convenient approach. Therefore, we have selected Autograd as the benchmark for comparison in this section.

Fig. 4 shows the approximation accuracy of HOPE and Autograd on different neural networks. Specifically, we calculated the first 10-order derivatives with HOPE and Autograd separately, and get the Taylor polynomials with Eq. (60). In Fig. 4(a), we compared the output curves of the 1-D network, the Taylor polynomial of HOPE, and Autograd. When all the modules of the network are \(n\)-order differentiable, HOPE can perform high local approximation on this neural network, while Autograd suffers from large deviation as the input moves far away from the reference point, which indicates that HOPE can get the high-order derivatives more accurately. When the network contains modules like ReLU and Max Pooling, both HOPE and Autograd can only obtain the first-order derivative.

In Fig. 4(b), we drew the input-to-output mapping surfaces and the residues of a 2D-CNN. One can observe that HOPE's output is closer to the target surface, as shown in the residual of both methods, further demonstrating HOPE's superior accuracy in calculating high-order derivatives and capability to perform local approximations to neural networks.

In Fig. 4(c), we trained an MNIST [41] 01 classifier, and get all the first 10-order unmixed partial derivatives on a certain input image with HOPE, while Autograd is unable to decompose this model because the input dimension is too large. We varied the intensity of point \(\mathbf{x}_{18,9}\) while keeping the other input points unchanged, and plot the input-to-output mapping of the network and the Taylor polynomial, showing the local approximation ability of HOPE. The top \(10\%\) positive factors of the first 10-order derivatives are shown on the right. Since previous works on network interpretability usually analyze the feature contribution from the first-order derivatives of the adopted neural network, while the visualized derivatives indicate that higher-order derivatives also reflect the influence of input on output, we will explore how to utilize higher-order heat maps in our future work.

### _Running Efficiency_

We test the efficiency of the proposed approach HOPE and Autograd, in terms of running time and memory costs. The experiments are conducted on the Windows operating
system with a GPU of NVIDIA GeForce RTX 3080 Ti, 32 GB memory, and 20 CPU cores. The results are shown in Tab. (1).

We test on three MLPs with different input dimensions (\(1\sim 3\) ) but the same structure (10-layer MLP with a width of 1024), and an MNIST \(01\) classifier. We calculated all the mixed partial derivatives of the \(1\sim 3\) dimension networks, and for the MNIST network, we only calculate its unmixed partial derivatives. For a \(p\)-D network, there're \(n\) (\(p=1\)) or \(\frac{p^{n+1}-p}{p-1}\) (\(p>1\)) mixed partial derivatives, and for the MNIST network, there are \(784n\) unmixed partial derivatives.

HOPE takes less than 1s and 6% of the available memory in all cases, while the time consumption of Autograd is much longer and increases exponentially with \(n\), or even out of memory (OOM). The significant superiority largely validates the time and memory efficiency of HOPE.

### _Convergence Under Different Parameter Settings_

Tab. (2) shows the influence of the parameter distribution of the target neural network on the convergence of its Taylor series. We initialize the weights of an MLP (width 512, depth 5) to follow a uniform distribution \(\mathbf{W}_{i,j}^{(m+1)}\sim U(-\frac{w_{0}}{\sigma_{m}},\frac{w_{0}}{\sigma_ {m}})\). As the value of \(w_{0}\) decreases, the parameters tend to concentrate more closely around zero, and it is more likely that the high-order derivatives of the model become increasingly smaller according to the inference in Section 4.2. The data in this table is the absolute value of the \(n\)-order derivative divided by the first-order derivative (i.e., \(\frac{|\partial^{2}f}{\partial x^{n}}/\frac{\partial f}{\partial x}\)). We varied \(w_{0}\) from 0.01 to 100. When \(w_{0}=0.01\) and \(w_{0}=0.1\), the high-order derivatives are much smaller than the low-order derivatives. When \(w_{0}=1\), almost all of the derivatives are on the same order of magnitude. When \(w_{0}=10\) and \(w_{0}=100\), the high-order derivatives are far larger than the low-order derivatives, which means that we cannot ignore the higher-order terms and make local approximations to neural networks.

### _Upper and Lower Bounds of a Neural Network and its Taylor Polynomial_

Based on Eq. 66, one can calculate the theoretical error bound between a neural network and its Taylor polynomial. In Fig. 5, the first panel shows the maximum approximation error \(e_{1}\) in the interval [-6,6] at different orders, and the theoretical upper bound of error \(e_{2}\). A small \(e_{1}\) or \(e_{2}\) indicates that the model has small prediction errors at each point. One can see that the theoretical error \(e_{2}\) is always larger than \(e_{1}\), and when \(n>14\) the magnitude of \(e_{2}\) reduces to a small value, resulting in a high approximation accuracy (small \(e_{1}\)).

Fig. 4: **Approximation results of HOPE and Autograd.** (a) Approximation curves of 1-D networks with Sine (middle left), ReLU (bottom left), Average Pooling (middle right), and Max Pooling (bottom right). (b) Approximation surfaces (left) and approximation residuals (right) of a 2-D network. (c) Expansion results of an MNIST 01 classifier. High-order heat maps were calculated using HOPE (right), and a comparison (left) was made between the network and the Taylor polynomial with \(\mathbf{x}_{18,9}\) changed while keeping other input points unchanged.
Other panels in Fig. 5 show the evolution of the approximation curves with increasing orders. \(f_{1}(x)\) and \(f_{2}(x)\) are the theoretical bounds of the network and its Taylor polynomial, and their expressions can be found in Eq. 63. As the degree of the approximation increases, the bounds \(f_{1}(x)\) and \(f_{2}(x)\) gradually converge, resulting in a closer approximation between this neural network and its Taylor polynomial.

### _Applications_

**Interpretation of deep neural networks and function discovery.** Implicit neural functions [40, 42, 43, 44, 45] are of strong expression capability and can be used to describe some unknown systems from observations. Taking the 2-D function

\[y=\frac{\mathbf{x}_{1}^{2}+\mathbf{x}_{2}}{2},\ \mathbf{x}1,\mathbf{x}2\in[-1,1] \tag{71}\]

as an example, we uniformly sample in the range \([-1,1]^{2}\) and then use a 2-D MLP to fit its observations. Since the MLP is a "black-box", we expand it into 2-order Taylor polynomials on reference inputs (0.0, 0.0), (0.5, 0.5), and (-0.5, -0.5) separately, and achieve following explicit expressions

\[\mathbf{y} =-0.01+0.00\mathbf{x}_{1}+0.51\mathbf{x}_{2}+0.55\mathbf{x}_{1}^ {2}-0.00\mathbf{x}_{1}\mathbf{x}_{2}+0.03\mathbf{x}_{2}^{2}\] \[\approx 0.51\mathbf{x}_{2}+0.55\mathbf{x}_{1}^{2},\] \[\mathbf{y} =0.38+0.54(\mathbf{x}_{1}-0.5)+0.51(\mathbf{x}_{2}-0.5)+0.53( \mathbf{x}_{1}-0.5)^{2}\] \[+0.06(\mathbf{x}_{1}-0.5)(\mathbf{x}_{2}-0.5)+0.03(\mathbf{x}_{2} -0.5)^{2}\] \[=0.01-0.02\mathbf{x}_{1}+0.45\mathbf{x}_{2}+0.53\mathbf{x}_{1}^ {2}+0.06\mathbf{x}_{1}\mathbf{x}_{2}+0.03\mathbf{x}_{2}^{2}\] \[\approx 0.45\mathbf{x}_{2}+0.53\mathbf{x}_{1}^{2},\] \[\mathbf{y} =-0.13-0.51(\mathbf{x}_{1}+0.5)+0.49(\mathbf{x}_{2}+0.5)+0.50( \mathbf{x}_{1}+0.5)^{2}\] \[-0.07(\mathbf{x}_{1}+0.5)(\mathbf{x}_{2}+0.5)+0.02(\mathbf{x}_{2} +0.5)^{2}\] \[=-0.03-0.05\mathbf{x}_{1}+0.48\mathbf{x}_{2}+0.50\mathbf{x}_{1}^ {2}-0.07\mathbf{x}_{1}\mathbf{x}_{2}+0.02\mathbf{x}_{2}^{2}\] \[\approx 0.48\mathbf{x}_{2}+0.50\mathbf{x}_{1}^{2}. \tag{72}\]

The aforementioned equations provide local explanations for the "black-box" network. When all these local explanations align and reach a consistent conclusion, a global explanation can be obtained. The findings suggest that HOPE possesses the capability of local interpretation and also exhibits potential for global interpretation and function discovery. Furthermore, the results validate that the expanded polynomial can learn the latent function with the same fidelity as the trained neural network, within the whole interval \([-1,1]^{2}\). This also implies that HOPE can be employed to assess the quality of the model.

**Low-cost inference of deep neural networks.** To show the advantageous running efficiency after Taylor expansion, we test on a controller of a single-tank liquid system implemented with a neural network. We simulated the following liquid system, in which the opening of the water outlet valve \(v_{2}\) keeps constant, while the water output \(q_{2}\) is determined by the liquid level height \(h\), as illustrated in Fig. 6(a). Specifically, to achieve the desired liquid level height \(h_{s}\), the opening of the inlet valve \(v_{1}\) is manipulated to regulate the quantity of inlet water \(q_{1}\), and the system is a first-order differential system, with transfer function

\[G(s)=\frac{K}{Ts+1}. \tag{73}\]

Here \(K\)=1 is the system's gain, \(T\)=1 is the system's time constant, and \(s\) is the complex frequency domain variable.

We trained a neural network level controller \(\hat{q}_{1}=f(e;\theta)\), where the input is the liquid level difference \(e=h_{s}-h\), and the output is the estimated input flow rate \(\hat{q}_{1}\). The network

\begin{table}
\begin{tabular}{c||c|c|c|c|c|c|c} \hline  & \multicolumn{2}{c|}{\(\mathbf{x}\in\mathbb{R}\)} & \multicolumn{2}{c|}{\(\mathbf{x}\in\mathbb{R}^{2}\)} & \multicolumn{2}{c|}{\(\mathbf{x}\in\mathbb{R}^{3}\)} & \multicolumn{2}{c}{\(\mathbf{x}\in\mathbb{R}^{28.728}\)} \\ \cline{2-9}
**n** & HOPE & Autograd & HOPE & Autograd & HOPE & Autograd & HOPE & Autograd \\ \hline
1 & 0.0059 / 0.20 & 0.0328 / 0.10 & 0.0109 / 0.30 & 0.0617 / 0.10 & 0.0120 / 0.20 & 0.0462 / 0.10 & 0.0120 / 0.10 & 0.0339 / 0.80 \\
2 & 0.0159 / 0.40 & 0.0358 / 0.10 & 0.0160 / 0.30 & 0.0820 / 0.10 & 0.0162 / 0.30 & 0.0921 / 0.20 & 0.0431 / 0.10 & 3.3456 / 2.10 \\
3 & 0.0229 / 0.40 & 0.0438 / 0.10 & 0.0232 / 0.30 & 0.1278 / 0.10 & 0.0234 / 0.40 & 0.2558 / 0.20 & 0.0229 / 0.10 & 7.5919 / 4.80 \\
4 & 0.0428 / 0.50 & 0.0538 / 0.10 & 0.0448 / 0.60 & 0.3559 / 0.20 & 0.0468 / 0.50 & 1.2151 / 0.90 & 0.0342 / 0.10 & 18.457 / 10.8 \\
5 & 0.0649 / 0.60 & 0.0812 / 0.10 & 0.0669 / 0.60 & 1.1938 / 0.80 & 0.0711 / 0.60 & 6.9751 / 5.80 & 0.0399 / 0.10 & 101.02 / 25.7 \\
6 & 0.0857 / 0.80 & 0.1446 / 0.10 & 0.0867 / 0.80 & 5.0142 / 0.40 & 0.0974 / 0.80 & 59.799 / 43.7 & 0.0479 / 0.10 & 1129 / 56.4 \\
7 & 0.1086 / 0.80 & 0.3037 / 0.20 & 0.1196 / 0.90 & 31.890 / 21.2 & 0.1340 / 1.00 & OOM / 76+ & 0.0648 / 0.10 & OOM / 76+ \\
8 & 0.1387 / 0.90 & 0.7244 / 0.50 & 0.1405 / 0.90 & OOM / 76+ & 0.1839 / 1.40 & OOM / 76+ & 0.0731 / 0.20 & OOM / 76+ \\
9 & 0.1634 / 1.00 & 1.9066 / 1.50 & 0.1638 / 1.10 & OOM / 76+ & 0.3234 / 2.40 & OOM / 76+ & 0.0864 / 0.30 & OOM / 76+ \\
10 & 0.1942 / 1.20 & 5.4099 / 4.10 & 0.1982 / 1.50 & OOM / 76+ & 0.7079 / 5.10 & OOM / 76+ & 0.1000 / 0.30 & OOM / 76+ \\ \hline \end{tabular}
\end{table} TABLE I: **Comparison of time (s) / memory (%) between HOPE and Autograd.** Here \(n\) is the order of derivatives. When the program is not running, the system memory usage is recorded at 24%. If an Out-of-Memory (OOM) error is encountered, it indicates that the program is utilizing 76% or more of the available memory.

\begin{table}
\begin{tabular}{c||c|c|c|c|c} \hline \(n\) & \(w_{0}\)=0.01 & \(w_{0}\)=0.1 & \(w_{0}\)=1.0 & \(w_{0}\)=10 & \(w_{0}\)=100 \\ \hline
1 & 1.00e+00 & 1.00e+00 & 1.00e+00 & 1.00e+00 & 1.00e+00 \\
2 & 4.69e-03 & 6.39e-02 & 7.36e-02 & 3.14e+01 & 2.98e+03 \\
3 & 3.94e-05 & 6.99e-03 & 7.98e-01 & 3.83e+01 & 6.08e+03 \\
4 & 3.12e-07 & 5.18e-04 & 1.10e-01 & 2.28e+03 & 5.27e+07 \\
5 & 1.97e-09 & 5.30e-05 & 7.01e-01 & 3.49e+03 & 1.66e+08 \\
6 & 2.15e-11 & 4.42
structure is shown in Fig. 6(b). Here we set \(h_{s}=10\), and the label for training is designed as

\[f(e)=\left\{\begin{array}{ll}20,&6<e\leq 10\\ 15,&2<e\leq 6\\ 10,&-2<e\leq 2\\ 5,&-6<e\leq-2\\ 0.&-10<e\leq-6\end{array}\right. \tag{74}\]

We normalize both the inputs and labels to the range of -1 to 1, and the loss function is defined as

\[J(e;\theta)=\|f(\frac{e}{10};\theta)-(\frac{f(e)}{10}-1)\|_{2}. \tag{75}\]

Due to the continuity of this neural network, \(f(e;\theta)\) will be relatively smooth and not perfectly fit the function \(f(e)\).

We expand \(f(e;\theta)\) to 3-order Taylor polynomials, within the range [-1, 1], to replace the neural network perfectly. We conducted liquid-level control experiments on the system using both the original neural network and the 3-order polynomial, under different initial liquid levels and initial inflow rates. The results are shown in Fig. 6(c), from which we can see that the polynomial can exactly replace the neural network.

We conducted statistical analysis on the inference time of the deep neural network and its polynomial approximation by HOPE, across different input batch sizes. The results in

Fig. 5: **Upper and lower bounds of a network and its Taylor polynomials of different orders.** The first figure shows the errors between a neural network and the \(n\)-order Taylor polynomials in the interval \([-6,6]\). Here \(e_{1}=\max_{\mathbf{x}\in[-6,6]}|f(\mathbf{x})-f(\mathbf{x};\theta)|\), and \(e_{2}=\frac{\max f(\mathbf{x};\theta)^{(n)}-\min f(\mathbf{x};\theta)^{(n)}} {n!}6^{n}\). \(f_{1}(x)\) and \(f_{2}(x)\) are respectively the theoretical bound of the neural network and Taylor polynomial.

Fig. 6: **Illustration of the single-tank liquid level control system.** (a) A single-tank system. (b) The structure of the deep neural network implementing the controller, with the name of each module corresponds to the name of the class invoked in PyTorch. (c) Liquid level control curves of the neural-network controller and its polynomial expansion under different initial liquid levels and inflow rates.

\begin{table}
\begin{tabular}{c||c|c|c|c|c|c|c} \hline Batch size & 1 & 4 & 16 & 64 & 256 & 1024 & 4096 \\ \hline Network & 0.2056 & 0.2640 & 0.4163 & 2.5155 & 10.805 & 36.401 & 121.98 \\ HOPE & 0.0284 & 0.0284 & 0.0284 & 0.0284 & 0.0291 & 0.0294 & 0.0324 \\ \hline \end{tabular}
\end{table} TABLE III: Inference time (ms) under different input batches.
Tab. (3) indicate that the expanded polynomial has significantly shorter and more concentrated inference time. When the input batch grows to 4096, the polynomial inference time is approximately 3765 times faster than that of the neural network, while maintaining the same prediction accuracy. In terms of the model size, the file size of the deep neural network reaches 44,671 bytes, while the Taylor polynomial merely occupies 160 bytes, showcasing a remarkable reduction in storage space.

**Feature selection.** For a neural network taking multiple inputs, an explicit expansion can efficiently measure the quantitative contribution of each input element to the final output, and facilitate feature selection. To demonstrate this application, we trained an MNIST handwritten digit classifier, and then separate it into ten equivalent single-output classifiers for easier expansion, as illustrated in Fig. 7(a). Denoting the input image as \(\mathbf{x}\), prediction as \(\mathbf{y}\), and \(\mathbf{x}_{i,j}\)'s impact on \(\mathbf{y}\) as \(\Delta\mathbf{y}_{i,j}\in\mathbb{R}\), we have

\[\Delta\mathbf{y}_{i,j}\approx\sum_{k=1}^{n}\frac{\frac{\partial^{k}\mathbf{y} }{\partial\mathbf{x}_{i,j}^{k}}}{k!}\Delta\mathbf{x}_{i,j}^{k}, \tag{76}\]

which can be further converted into matrix form

\[\Delta\mathbf{y}\approx\sum_{k=1}^{n}\frac{\frac{\beta^{k}\mathbf{y}}{\partial \mathbf{x}^{k}}}{k!}\odot\Delta\mathbf{x}^{\circ k}. \tag{77}\]

Here \(\Delta\mathbf{x}\in\mathbb{R}^{28\times 28}\) represents the perturbation applied to the input image, \(\frac{\beta^{k}\mathbf{y}}{\partial\mathbf{x}^{k}}\in\mathbb{R}^{28\times 28}\) contains all the \(k\)-order unmixed partial derivatives, and \(\Delta\mathbf{y}\in\mathbb{R}^{28\times 28}\) is the heat map reflecting the impact of all input elements on the output.

We initialize \(\Delta\mathbf{x}\) as \(\mathbf{1}^{28\times 28}\) and expanded the ten single-output models to obtain 10 heat maps, and then applied perturbation analysis to get another 10 counterparts, as visualized in the upper and lower rows of Fig. 7(b). The results demonstrate that HOPE is capable of generating an equivalent heat map to the perturbation-based method. Besides, HOPE can infer significantly faster, taking only 0.002s to generate the heat map, whereas perturbation analysis requires 0.342s.

## 6 Conclusion

Aiming at providing a high precision polynomial interpretation of the "black-box" deep neural networks, we propose the network's high-order Taylor expansion, which is of high accuracy, low computational cost, good convergence and wide applications. Understanding the mechanism behind the deep neural networks is quite important, and we believe that neural networks will become more transparent with HOPE, accelerating the development and application of neural networks.

**Summary of the approach.** Specifically, we first derive the high-order derivative rule for a general composite function and then extend the rule to neural networks for fast and accurate calculation of its high-order derivatives. From all above derivatives, we can expand a black-box network into an explicit Taylor polynomial, providing a local explanation for the network's mapping from the input to the output. We also theoretically prove that a neural network is equivalent to its infinite Taylor polynomial if all of the modules are infinitely differential, and analyze the polynomial's convergence condition as well.

**Advantageous and applications.** HOPE has significant advantages in terms of accuracy, speed, and memory cost compared with computational-graph-based method. It works as a general expansion and thus of wide applicability for diverse deep neural networks, e.g., with different dimensions and layers.

The explicit Taylor expansion possesses the ability to conduct data-driven explicit model identification, facilitate the fast inference of a trained deep neural network, and select informative features contributing to the output.

**Limitations and future work.** It should be noted that HOPE can be used only for modules that are \(n\)-order differentiable. For networks including components such as ReLU, LeakReLU, or Max Pooling, both HOPE and Autograd can only obtain their first-order information. Further increasing the applicability is one of our ongoing work.

In the future, we would like to explore deeper the relationship between the convergence of the Taylor series and the parameter distribution, and apply it to the optimization and structure design of deep neural networks. We will also explore how to use high-order heat maps to determine the contribution of inputs to outputs more accurately. Moreover, HOPE can get the derivatives between any nodes of a neural network, which might inspire lightweight network design.

## Acknowledgments

This work is jointly funded by National Natural Science Foundation of China (Grant No. 61931012) and Beijing Municipal Natural Science Foundation (Grant No. Z200021).
